Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:20:33

0001 // use #include "" only for your local include and put
0002 // those in the first line(s) before any #include <>
0003 // otherwise you are asking for weird behavior
0004 // (more info - check the difference in include path search when using "" versus <>)
0005 
0006 #include "ZdcMon.h"
0007 
0008 #include <onlmon/OnlMon.h>  // for OnlMon
0009 #include <onlmon/OnlMonDB.h>
0010 #include <onlmon/OnlMonServer.h>
0011 
0012 #include <Event/msg_profile.h>
0013 #include <calobase/TowerInfoDefs.h>
0014 #include <caloreco/CaloWaveformFitting.h>
0015 
0016 #include <Event/Event.h>
0017 #include <Event/EventTypes.h>
0018 #include <Event/msg_profile.h>
0019 
0020 #include <TH1.h>
0021 #include <TH2.h>
0022 #include <TProfile.h>
0023 #include <TRandom.h>
0024 
0025 #include <cmath>
0026 #include <cstdio>  // for printf
0027 #include <fstream>
0028 #include <iostream>
0029 #include <sstream>
0030 #include <string>  // for allocator, string, char_traits
0031 
0032 enum
0033 {
0034   TRGMESSAGE = 1,
0035   FILLMESSAGE = 2
0036 };
0037 
0038 ZdcMon::ZdcMon(const std::string &name)
0039   : OnlMon(name)
0040 {
0041   // leave ctor fairly empty, its hard to debug if code crashes already
0042   // during a new ZdcMon()
0043   return;
0044 }
0045 
0046 ZdcMon::~ZdcMon()
0047 {
0048   // you can delete NULL pointers it results in a NOOP (No Operation)
0049   return;
0050 }
0051 
0052 int ZdcMon::Init()
0053 {
0054   const float MAX_ENERGY1 = 12000.;
0055   const float MAX_ENERGY2 = 30000.;
0056   const float MAX_WFAMP = 20000.;
0057   const float MIN_ENERGY1 = 0.;
0058   const float MIN_ENERGY2 = 0.;
0059   const int BIN_NUMBER1 = 250;
0060   const int BIN_NUMBER2 = 550;
0061   const int SMD_ADC_BIN = 360;
0062   const float MAX_SMD_ADC = 18000.;
0063   const int BIN_WF = 1000;
0064     
0065   //  gRandom->SetSeed(rand());
0066   // read our calibrations from ZdcMonData.dat
0067   const char *zdccalib = getenv("ZDCCALIB");
0068   if (!zdccalib)
0069   {
0070     std::cout << "ZDCCALIB environment variable not set" << std::endl;
0071     exit(1);
0072   }
0073   std::string fullfile = std::string(zdccalib) + "/" + "ZdcMonData.dat";
0074   std::ifstream calib(fullfile);
0075   calib.close();
0076 
0077   //getting gains
0078   float col1, col2, col3;
0079   std::string gainfile = std::string(zdccalib) + "/" + "/ZdcCalib.pmtgain";
0080   std::ifstream gain_infile(gainfile);
0081   std::ostringstream msg(gainfile);
0082 
0083   if (!gain_infile)
0084   {
0085     msg << gainfile << " could not be opened.";
0086     OnlMonServer *se = OnlMonServer::instance();
0087     se->send_message(this, MSG_SOURCE_ZDC, MSG_SEV_FATAL, msg.str(), 2);
0088     exit(1);
0089   }
0090 
0091   for (float &i : gain)
0092   {
0093     gain_infile >> col1 >> col2 >> col3;
0094     i = col1;
0095   }
0096 
0097   for (int i = 0; i < 16; i++)  // relative gains of SMD north channels
0098   {
0099     smd_south_rgain[i] = gain[i];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0100   }
0101 
0102   for (int i = 0; i < 16; i++)  // relative gains of SMD north channels
0103   {
0104     smd_north_rgain[i] = gain[i + 16];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0105   }
0106 
0107   gain_infile.close();
0108 
0109   // use printf for stuff which should go the screen but not into the message
0110   // system (all couts are redirected)
0111   printf("doing the Init\n");
0112 
0113   // Create hitograms
0114   zdc_adc_north = new TH1F("zdc_adc_north", "ZDC ADC north", BIN_NUMBER2, MIN_ENERGY2, MAX_ENERGY2);
0115   zdc_adc_south = new TH1F("zdc_adc_south", "ZDC ADC south", BIN_NUMBER2, MIN_ENERGY2, MAX_ENERGY2);
0116 
0117   zdc_N1 = new TH1F("zdc_N1", "ZDC1 ADC north", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0118   zdc_N2 = new TH1F("zdc_N2", "ZDC2 ADC north", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0119   zdc_N3 = new TH1F("zdc_N3", "ZDC3 ADC north", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0120   zdc_S1 = new TH1F("zdc_S1", "ZDC1 ADC south", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0121   zdc_S2 = new TH1F("zdc_S2", "ZDC2 ADC south", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0122   zdc_S3 = new TH1F("zdc_S3", "ZDC3 ADC south", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0123 
0124   veto_NF = new TH1F("veto_NF", "veto north front", SMD_ADC_BIN, MIN_ENERGY2, MAX_SMD_ADC);
0125   veto_NB = new TH1F("veto_NB", "veto north back", SMD_ADC_BIN, MIN_ENERGY2, MAX_SMD_ADC);
0126   veto_SF = new TH1F("veto_SF", "veto south front", SMD_ADC_BIN, MIN_ENERGY2, MAX_SMD_ADC);
0127   veto_SB = new TH1F("veto_SB", "veto south back", SMD_ADC_BIN, MIN_ENERGY2, MAX_SMD_ADC);
0128 
0129   //waveform
0130     
0131    h_waveform_timez = new TH1F("h_waveform_timez", "", 16, 0.5, 16 + 0.5);
0132    h_waveform_timess = new TH1F("h_waveform_timess", "", 16, 0.5, 16 + 0.5);
0133    h_waveform_timesn = new TH1F("h_waveform_timesn", "", 16, 0.5, 16 + 0.5);
0134    h_waveform_timevs = new TH1F("h_waveform_timevs", "", 16, 0.5, 16 + 0.5);
0135    h_waveform_timevn = new TH1F("h_waveform_timevn", "", 16, 0.5, 16 + 0.5);
0136 
0137   h_waveformZDC = new TH2F("h_waveformZDC", "h_waveformZDC", 31, 0.5, 31 + 0.5, BIN_WF, 0, MAX_WFAMP);
0138   h_waveformSMD_North = new TH2F("h_waveformSMD_North", "h_waveformSMD_North", 31, 0.5, 31 + 0.5, BIN_WF, 0, MAX_WFAMP);
0139   h_waveformSMD_South = new TH2F("h_waveformSMD_South", "h_waveformSMD_South", 31, 0.5, 31 + 0.5, BIN_WF, 0, MAX_WFAMP);
0140   h_waveformVeto_North = new TH2F("h_waveformVeto_North", "h_waveformVeto_North", 31, 0.5, 31 + 0.5, BIN_WF, 0, MAX_WFAMP);
0141   h_waveformVeto_South = new TH2F("h_waveformVeto_South", "h_waveformVeto_South", 31, 0.5, 31 + 0.5, BIN_WF, 0, MAX_WFAMP);
0142       
0143   // SMD
0144   // Individual SMD_ADC Values
0145   // Horizontal (expert plot)
0146   for (int i = 0; i < 8; i++)
0147   {
0148     smd_adc_n_hor_ind[i] = new TH1I(Form("smd_adc_n_hor_ind%d", i), Form("smd_adc_n_hor_ind%d", i), SMD_ADC_BIN, 0, MAX_SMD_ADC);
0149     smd_adc_s_hor_ind[i] = new TH1I(Form("smd_adc_s_hor_ind%d", i), Form("smd_adc_s_hor_ind%d", i), SMD_ADC_BIN, 0, MAX_SMD_ADC);
0150   }
0151   // Vertical (expert plot)
0152   for (int i = 0; i < 7; i++)
0153   {
0154     smd_adc_n_ver_ind[i] = new TH1I(Form("smd_adc_n_ver_ind%d", i), Form("smd_adc_n_ver_ind%d", i), SMD_ADC_BIN, 0, MAX_SMD_ADC);
0155     smd_adc_s_ver_ind[i] = new TH1I(Form("smd_adc_s_ver_ind%d", i), Form("smd_adc_s_ver_ind%d", i), SMD_ADC_BIN, 0, MAX_SMD_ADC);
0156   }
0157 
0158   // SMD Hit Multiplicity
0159   smd_north_hor_hits = new TH1F("smd_north_hor_hits", "smd_north_hor_hits", 9, -0.5, 8.5);
0160   smd_north_ver_hits = new TH1F("smd_north_ver_hits", "smd_north_ver_hits", 8, -0.5, 7.5);
0161   smd_south_hor_hits = new TH1F("smd_south_hor_hits", "smd_south_hor_hits", 9, -0.5, 8.5);
0162   smd_south_ver_hits = new TH1F("smd_south_ver_hits", "smd_south_ver_hits", 8, -0.5, 7.5);
0163 
0164   // north smd
0165   smd_hor_north = new TH1F("smd_hor_north", "Beam centroid distribution, SMD North y", 296, -5.92, 5.92);
0166   smd_ver_north = new TH1F("smd_ver_north", "Beam centroid distribution, SMD North x", 220, -5.5, 5.5);
0167   
0168   smd_sum_hor_north = new TH1F("smd_sum_hor_north", "SMD North y", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0169   smd_sum_ver_north = new TH1F("smd_sum_ver_north", "SMD North x", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0170 
0171   smd_hor_north_small = new TH1F("smd_hor_north_small", "Beam centroid distribution, SMD North y, zdc <= 200", 296, -5.92, 5.92);
0172   smd_ver_north_small = new TH1F("smd_ver_north_small", "Beam centroid distribution, SMD North x, zdc <= 200", 220, -5.5, 5.5);
0173   smd_hor_north_good = new TH1F("smd_hor_north_good", "Beam centroid distribution, SMD North y, zdc > 200", 296, -5.92, 5.92);
0174   smd_ver_north_good = new TH1F("smd_ver_north_good", "Beam centroid distribution, SMD North x, zdc > 200", 220, -5.5, 5.5);
0175 
0176   // south smd
0177   smd_hor_south = new TH1F("smd_hor_south", "Beam centroid distribution, SMD South y", 296, -5.92, 5.92);
0178   smd_ver_south = new TH1F("smd_ver_south", "Beam centroid distribution, SMD South x", 220, -5.5, 5.5);
0179 
0180   smd_sum_hor_south = new TH1F("smd_sum_hor_south", "SMD South y", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0181   smd_sum_ver_south = new TH1F("smd_sum_ver_south", "SMD South x", BIN_NUMBER1, MIN_ENERGY1, MAX_ENERGY1);
0182   
0183   smd_hor_south_good = new TH1F("smd_hor_south_good", "Beam centroid distribution, SMD South y, zdc1 > 65 zdc2>20 and veto<200", 296, -5.92, 5.92);
0184   smd_ver_south_good = new TH1F("smd_ver_south_good", "Beam centroid distribution, SMD South x, zdc1 > 65 zdc2>20 and veto<200", 220, -5.5, 5.5);
0185 
0186 
0187   // smd values (expert plot)
0188   smd_value = new TH2F("smd_value", "SMD channel# vs value", 1024, 0, 4096, 33, -0.5, 32.5);
0189   smd_value_good = new TH2F("smd_value_good", "SMD channel# vs value, zdc > 200", 1024, 0, 4096, 33, -0.5, 32.5);
0190   smd_value_small = new TH2F("smd_value_small", "SMD channel# vs value, zdc <= 200", 1024, 0, 4096, 33, -0.5, 32.5);
0191   smd_xy_north = new TH2F("smd_xy_north", "SMD hit position north", 110, -5.5, 5.5, 119, -5.92, 5.92);
0192   smd_xy_south = new TH2F("smd_xy_south", "SMD hit position south", 110, -5.5, 5.5, 119, -5.92, 5.92);
0193 
0194   OnlMonServer *se = OnlMonServer::instance();
0195 
0196   //register histograms with server otherwise client won't get them
0197   //zdc
0198   se->registerHisto(this, zdc_adc_north);
0199   se->registerHisto(this, zdc_adc_south);
0200   se->registerHisto(this, zdc_N1);
0201   se->registerHisto(this, zdc_N2);
0202   se->registerHisto(this, zdc_N3);
0203   se->registerHisto(this, zdc_S1);
0204   se->registerHisto(this, zdc_S2);
0205   se->registerHisto(this, zdc_S3);
0206   se->registerHisto(this, h_waveformZDC);
0207   se->registerHisto(this, h_waveformSMD_North);
0208   se->registerHisto(this, h_waveformSMD_South);
0209   se->registerHisto(this, h_waveformVeto_North);
0210   se->registerHisto(this, h_waveformVeto_South);
0211     
0212   se->registerHisto(this, h_waveform_timez);
0213   se->registerHisto(this, h_waveform_timess);
0214   se->registerHisto(this, h_waveform_timesn);
0215   se->registerHisto(this, h_waveform_timevs);
0216   se->registerHisto(this, h_waveform_timevn);
0217  
0218     
0219   //veto
0220   se->registerHisto(this, veto_NF);
0221   se->registerHisto(this, veto_NB);
0222   se->registerHisto(this, veto_SF);
0223   se->registerHisto(this, veto_SB);
0224  
0225 
0226   // SMD
0227   // Individual smd_adc channel histos
0228 
0229   for (int i = 0; i < 8; i++)
0230   {
0231     se->registerHisto(this, smd_adc_n_hor_ind[i]);
0232     se->registerHisto(this, smd_adc_s_hor_ind[i]);
0233   }
0234   for (int i = 0; i < 7; i++)
0235   {
0236     se->registerHisto(this, smd_adc_n_ver_ind[i]);
0237     se->registerHisto(this, smd_adc_s_ver_ind[i]);
0238   }
0239 
0240   // SMD Hit Multiplicity
0241   se->registerHisto(this, smd_north_hor_hits);
0242   se->registerHisto(this, smd_north_ver_hits);
0243   se->registerHisto(this, smd_south_hor_hits);
0244   se->registerHisto(this, smd_south_ver_hits);
0245   
0246 
0247   // north SMD
0248   se->registerHisto(this, smd_hor_north);
0249   se->registerHisto(this, smd_ver_north);
0250   se->registerHisto(this, smd_sum_hor_north);
0251   se->registerHisto(this, smd_sum_ver_north);
0252   se->registerHisto(this, smd_hor_north_small);
0253   se->registerHisto(this, smd_ver_north_small);
0254   se->registerHisto(this, smd_hor_north_good);
0255   se->registerHisto(this, smd_ver_north_good);
0256   // south SMD
0257   se->registerHisto(this, smd_hor_south);
0258   se->registerHisto(this, smd_ver_south);
0259   se->registerHisto(this, smd_sum_hor_south);
0260   se->registerHisto(this, smd_sum_ver_south);
0261   se->registerHisto(this, smd_hor_south_good);
0262   se->registerHisto(this, smd_ver_south_good);
0263   // SMD values
0264   se->registerHisto(this, smd_value);
0265   se->registerHisto(this, smd_value_good);
0266   se->registerHisto(this, smd_value_small);
0267   se->registerHisto(this, smd_xy_north);
0268   se->registerHisto(this, smd_xy_south);
0269 
0270   WaveformProcessingFast = new CaloWaveformFitting();
0271 
0272   Reset();
0273 
0274   return 0;
0275 }
0276 
0277 int ZdcMon::BeginRun(const int /* runno */)
0278 {
0279   // if you need to read calibrations on a run by run basis
0280   // this is the place to do it
0281 
0282   return 0;
0283 }
0284 
0285 std::vector<float> ZdcMon::anaWaveformFast(Packet *p, const int channel)
0286 {
0287   std::vector<float> waveform;
0288   // Chris: preallocation = speed improvement
0289   waveform.reserve(p->iValue(0, "SAMPLES"));
0290   for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0291   {
0292     waveform.push_back(p->iValue(s, channel));
0293   }
0294   std::vector<std::vector<float>> multiple_wfs;
0295   multiple_wfs.push_back(waveform);
0296 
0297   std::vector<std::vector<float>> fitresults_zdc;
0298   fitresults_zdc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
0299 
0300   std::vector<float> result;
0301   result = fitresults_zdc.at(0);
0302   return result;
0303 }
0304 
0305 int ZdcMon::process_event(Event *e /* evt */)
0306 {
0307   evtcnt++;
0308     
0309   const float waveform_hit_threshold = 100 * 5.;
0310 
0311   int packet = 12001;
0312 
0313   float totalzdcsouthsignal = 0.;
0314   float totalzdcnorthsignal = 0.;
0315   float smd_south_xsum = 0.;
0316   float smd_south_ysum = 0.;
0317   float smd_north_xsum = 0.;
0318   float smd_north_ysum = 0.;
0319 
0320   float veto_cut = 200.0;
0321   float ZDC1cut = 65.0;
0322   float ZDC2cut = 20.0;
0323 
0324   std::vector <float> z;
0325   z.clear();
0326     
0327   std::vector <float> sm;
0328   sm.clear();
0329     
0330   std::vector <float> tz;
0331   tz.clear();
0332     
0333   std::vector <float> tsmd;
0334   tsmd.clear();
0335     
0336   std::vector<float> resultFast;
0337   resultFast.clear();
0338 
0339   float zdctimelow  = 5.0;
0340   float zdctimehigh  = 9.0;
0341     
0342   float smdNtimelow  = 9.0;
0343   float smdNtimehigh  = 14.0;
0344 
0345   float smdStimelow  = 6.0;
0346   float smdStimehigh  = 12.0;
0347     
0348   float vetoStimelow  = 6.0;
0349   float vetoStimehigh  = 12.0;
0350     
0351   float vetoNtimelow  = 5.0;
0352   float vetoNtimehigh  = 9.0;
0353  
0354   Packet *p = e->getPacket(packet);
0355   if (p)
0356   {
0357         
0358     for (int c = 0; c < p->iValue(0, "CHANNELS"); c++)
0359     {
0360       resultFast = anaWaveformFast(p, c);  // fast waveform fitting
0361       float signalFast = resultFast.at(0);
0362       float time = resultFast.at(1);
0363       float pedestal = resultFast.at(2);
0364       float signal = signalFast;
0365         
0366         for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0367         {
0368             if (c < 16)
0369             {
0370                 if (signal > waveform_hit_threshold) h_waveformZDC->Fill(s, p->iValue(s, c) - pedestal);
0371             }
0372             
0373             if (c > 15 && c < 18)
0374             {
0375                 if (signal > waveform_hit_threshold) h_waveformVeto_North->Fill(s, p->iValue(s, c) - pedestal);
0376             }
0377             
0378             if (c > 47 && c < 64)
0379             {
0380                 if (signal > waveform_hit_threshold) h_waveformSMD_North->Fill(s, p->iValue(s, c) - pedestal);
0381             }
0382             
0383             if (c > 79 && c < 82)
0384             {
0385                 if (signal > waveform_hit_threshold) h_waveformVeto_South->Fill(s, p->iValue(s, c) - pedestal);
0386             }
0387             
0388             if (c > 111)
0389             {
0390                 if (signal > waveform_hit_threshold) h_waveformSMD_South->Fill(s, p->iValue(s, c) - pedestal);
0391             }
0392         }
0393 
0394      // fill zdc raw signal first
0395      if(c < 16)
0396      {
0397          if(signal > waveform_hit_threshold)
0398          {
0399            h_waveform_timez->Fill(time);
0400          }
0401  
0402        z.push_back(signal);
0403        tz.push_back(time);
0404 
0405       if((c == 0) || (c == 2) || (c == 4))
0406       {
0407          totalzdcsouthsignal += signal;
0408           
0409          if(c == 0) zdc_S1->Fill(signal);
0410          if(c == 2) zdc_S2->Fill(signal);
0411          if(c == 4) zdc_S3->Fill(signal);
0412       }
0413 
0414       else if((c==8) || (c ==10) || (c == 12))
0415       {
0416          totalzdcnorthsignal += signal;
0417           
0418          if(c == 8) zdc_N1->Fill(signal);
0419          if(c == 10) zdc_N2->Fill(signal);
0420          if(c == 12) zdc_N3->Fill(signal);
0421       }
0422      }
0423 
0424     // veto N
0425     else if(c > 15 && c < 18)
0426     {
0427         if(signal > waveform_hit_threshold)
0428         {
0429            h_waveform_timevn->Fill(time);
0430         }
0431         
0432         if(((time >= vetoNtimelow) && (time <= vetoNtimehigh)))
0433         {
0434 
0435             if(c == 16)
0436             {
0437                 v[0] = signal; if(signal > 0.) veto_NF->Fill(signal);
0438             }
0439             
0440             if (c == 17)
0441             {
0442                 v[1] = signal; if(signal > 0.) veto_NB->Fill(signal);
0443             }
0444         }
0445         
0446     }
0447         
0448      //smd N
0449      else if(c > 47 && c < 64)
0450      {
0451         sm.push_back(signal);
0452         tsmd.push_back(time);
0453          
0454          if(signal > waveform_hit_threshold)
0455          {
0456             h_waveform_timesn->Fill(time);
0457          }
0458      }
0459         
0460     // veto S
0461     else if(c > 79 && c < 82)
0462     {
0463         if(((time >= vetoStimelow) && (time <= vetoStimehigh)))
0464         {
0465 
0466             if(c == 80)
0467             {
0468                 v[2] = signal; if(signal > 0.) veto_SF->Fill(signal);
0469             }
0470                 
0471             if (c == 81)
0472             {
0473                 v[3] = signal; if(signal > 0.) veto_SB->Fill(signal);
0474             }
0475         }
0476         
0477         if(signal > waveform_hit_threshold)
0478         {
0479             h_waveform_timevs->Fill(time);
0480         }
0481         
0482             
0483     }
0484         
0485           
0486     //smd S
0487      else if(c > 111)
0488      {
0489            sm.push_back(signal);
0490            tsmd.push_back(time);
0491          
0492          if(signal > waveform_hit_threshold)
0493          {
0494              h_waveform_timess->Fill(time);
0495          }
0496      }
0497 
0498 
0499   }// channel loop end
0500 
0501       
0502    zdc_adc_south->Fill(totalzdcsouthsignal);
0503    zdc_adc_north->Fill(totalzdcnorthsignal);
0504 
0505         
0506   int zsize = z.size();
0507   int ssize = sm.size();
0508       
0509    if(zsize != 16)
0510     {
0511       std::cout<< "zdc channel mapping error" << std::endl;
0512       if(tz.size() != 16)
0513       exit(1);
0514     }
0515 
0516    if(ssize != 32)
0517     {
0518       std::cout<< "smd channel mapping error" << std::endl;
0519       if(tsmd.size() != 32)
0520       exit(1);
0521     }
0522 
0523 
0524      for (int i = 0; i < zsize; i++) //zdc
0525      {
0526          if ((tz[i] >= zdctimelow) && (tz[i] <= zdctimehigh))
0527          {
0528              zdc_adc[i] = z[i];
0529          }
0530          else
0531          {
0532             zdc_adc[i] = 0.0;
0533          }
0534      }//zdc
0535 
0536      for (int j = 0; j < ssize; j++) //smd
0537       {
0538           if(j < 16) //smd north [0,15] --> [48,63]
0539           {
0540             if((tsmd[j] >= smdNtimelow) && (tsmd[j] <= smdNtimehigh))
0541               {
0542                   
0543                   smd_adc[j] = sm[j];
0544                   if(j<=7) smd_north_ysum += sm[j];
0545                   if(j >= 8 && j<=14) smd_north_xsum += sm[j]; //skip sum ch, 15->63
0546               }
0547               else
0548               {
0549                  smd_adc[j] = 0.0;
0550               }
0551            }
0552         
0553           if (j >= 16 && j <= 31) //smd south [16,31] --> [112,127]
0554           {
0555              if((tsmd[j] >= smdStimelow) && (tsmd[j] <= smdStimehigh ))
0556               {
0557                   smd_adc[j] = sm[j];
0558                   if(j >= 16 && j<=23) smd_south_ysum += sm[j];
0559                   if(j >= 24 && j<=30) smd_south_xsum += sm[j]; //skip sum ch, 31->127
0560               }
0561               else
0562               {
0563                   smd_adc[j] = 0.0;
0564               }
0565           }
0566       }//smd
0567 
0568 
0569 
0570   // call the functions
0571     CompSmdAdc();
0572     CompSmdPos();
0573 
0574     // BOOLEANS, INTs AND OTHER DEFINITIONS
0575 
0576     bool fill_hor_south = false;
0577     bool fill_ver_south = false;
0578 
0579     bool fill_hor_north = false;
0580     bool fill_ver_north = false;
0581 
0582     int s_ver = 0;
0583     int s_hor = 0;
0584 
0585     int n_ver = 0;
0586     int n_hor = 0;
0587 
0588     float zdc_adc_threshold = 200.;
0589     float smd_adc_threshold = 20.;
0590    
0591     // counters
0592     int smd_n_h_counter = 0;
0593     int smd_n_v_counter = 0;
0594     int smd_s_h_counter = 0;
0595     int smd_s_v_counter = 0;
0596     
0597       
0598     float hor_cut = 8.0;
0599     float ver_cut = 8.0;
0600     float ped_cut = 0.0; //suppress zeroes
0601 
0602     for (int i = 0; i < 8; i++)
0603     {
0604 
0605      if (smd_adc[i] > hor_cut)
0606       {
0607         n_hor++;
0608       }
0609       
0610 
0611      if (smd_adc[i] > ped_cut) smd_adc_n_hor_ind[i]->Fill(smd_adc[i]);
0612 
0613      smd_value->Fill(smd_adc[i], float(i));
0614    
0615      if (zdc_adc[8] > veto_cut)
0616       {
0617          smd_value_good->Fill(smd_adc[i], float(i));
0618       }
0619       else
0620       {
0621          smd_value_small->Fill(smd_adc[i], float(i));
0622       }
0623 
0624       if ((smd_adc[i] > smd_adc_threshold) && (zdc_adc[0] > zdc_adc_threshold))
0625       {
0626         smd_n_h_counter++;
0627       }
0628 
0629       if (smd_adc[i] > smd_adc_threshold)
0630       {
0631         smd_n_h_counter++;
0632       }
0633 
0634       if (smd_adc[i + 16] > hor_cut)
0635       {
0636         s_hor++;
0637       }
0638     
0639       if (smd_adc[i + 16] > ped_cut) smd_adc_s_hor_ind[i]->Fill(smd_adc[i + 16]);
0640       
0641       smd_value->Fill(smd_adc[i + 16], float(i) + 16);
0642      
0643       if (zdc_adc[0] > veto_cut)
0644       {
0645         smd_value_good->Fill(smd_adc[i + 16], float(i) + 16);
0646       }
0647       else
0648       {
0649         smd_value_small->Fill(smd_adc[i + 16], float(i) + 16);
0650       }
0651 
0652 
0653       if (smd_adc[i + 16] > smd_adc_threshold)
0654       {
0655         smd_s_h_counter++;
0656       }
0657         
0658     }
0659 
0660     for (int i = 0; i < 7; i++)
0661     {
0662       if (smd_adc[i + 8] > ver_cut)
0663       {
0664         n_ver++;
0665       }
0666      
0667   
0668       if(smd_adc[i + 8] > ped_cut) smd_adc_n_ver_ind[i]->Fill(smd_adc[i + 8]);
0669 
0670       smd_value->Fill(smd_adc[i + 8], float(i) + 8);
0671    
0672       if (zdc_adc[8] > veto_cut)
0673       {
0674         smd_value_good->Fill(smd_adc[i + 8], float(i) + 8);
0675       }
0676       else
0677       {
0678         smd_value_small->Fill(smd_adc[i + 8], float(i) + 8);
0679       }
0680 
0681       if (smd_adc[i + 8] > smd_adc_threshold)
0682       {
0683         smd_n_v_counter++;
0684       }
0685       //****************************
0686 
0687       //****smd south vertical individual channels****
0688       if (smd_adc[i + 24] > ver_cut)
0689       {
0690         s_ver++;
0691       }
0692       
0693       if (smd_adc[i + 24] > ped_cut)   smd_adc_s_ver_ind[i]->Fill(smd_adc[i + 24]);
0694 
0695       smd_value->Fill(smd_adc[i + 24], float(i) + 24);
0696     
0697       if (zdc_adc[0] > veto_cut)
0698       {
0699         smd_value_good->Fill(smd_adc[i + 24], float(i) + 24);
0700       }
0701       else
0702       {
0703         smd_value_small->Fill(smd_adc[i + 24], float(i) + 24);
0704       }
0705 
0706       if (smd_adc[i + 24] > smd_adc_threshold)
0707       {
0708         smd_s_v_counter++;
0709       }
0710       //****************************
0711 
0712       // Fill out the SMD counters with doubles instead of integers.
0713       double nh = smd_n_h_counter + 0.0;
0714       smd_north_hor_hits->Fill(nh);
0715       double nv = smd_n_v_counter + 0.0;
0716       smd_north_ver_hits->Fill(nv);
0717       double sh = smd_s_h_counter + 0.0;
0718       smd_south_hor_hits->Fill(sh);
0719       double sv = smd_s_v_counter + 0.0;
0720       smd_south_ver_hits->Fill(sv);
0721     }
0722 
0723     
0724     bool fired_smd_hor_n = (n_hor > 1);
0725     bool fired_smd_ver_n = (n_ver > 1);
0726 
0727     bool fired_smd_hor_s = (s_hor > 1);
0728     bool fired_smd_ver_s = (s_ver > 1);
0729 
0730    //compute, if smd is overloaded
0731     bool smd_ovld_north = false;
0732     bool smd_ovld_south = false;
0733 
0734     // FILLING OUT THE HISTOGRAMS
0735     
0736 
0737     if (fired_smd_hor_s && fired_smd_ver_s && !smd_ovld_south)
0738     {
0739       fill_hor_south = true;
0740       fill_ver_south = true;
0741       smd_hor_south->Fill(smd_pos[2]);
0742       smd_ver_south->Fill(smd_pos[3]);
0743       smd_sum_ver_south->Fill(smd_south_xsum);
0744       smd_sum_hor_south->Fill(smd_south_ysum);
0745     }
0746 
0747     if(fill_hor_south && fill_ver_south && (zdc_adc[0] > ZDC1cut) && (zdc_adc[2] > ZDC2cut) && (v[2] < veto_cut) && (v[3] < veto_cut))
0748     {
0749       smd_xy_south->Fill(smd_pos[3], smd_pos[2]);
0750       smd_hor_south_good->Fill(smd_pos[2]);
0751       smd_ver_south_good->Fill(smd_pos[3]);
0752     }
0753 
0754    
0755     if (fired_smd_hor_n && fired_smd_ver_n && !smd_ovld_north)
0756     {
0757       fill_hor_north = true;
0758       fill_ver_north = true;
0759       smd_hor_north->Fill(smd_pos[0]);
0760       smd_ver_north->Fill(smd_pos[1]);
0761       smd_sum_ver_north->Fill(smd_north_xsum);
0762       smd_sum_hor_north->Fill(smd_north_ysum);
0763       
0764     }
0765 
0766     if(fill_hor_north && fill_ver_north && (zdc_adc[8] > ZDC1cut) && (zdc_adc[10] > ZDC2cut) && (v[0] < veto_cut) && (v[1] < veto_cut))
0767     {
0768       smd_xy_north->Fill(smd_pos[1], smd_pos[0]);
0769       smd_hor_north_good->Fill(smd_pos[0]);
0770       smd_ver_north_good->Fill(smd_pos[1]);
0771     }
0772 
0773    
0774   }  // if packet good
0775 
0776   z.clear();
0777   sm.clear();
0778   tz.clear();
0779   tsmd.clear();
0780     
0781   delete p;
0782   return 0;
0783 }
0784 
0785 int ZdcMon::Reset()
0786 {
0787   // reset our internal counters
0788   evtcnt = 0;
0789   idummy = 0;
0790   return 0;
0791 }
0792 
0793 void ZdcMon::CompSmdAdc()  // mulitplying by relative gains
0794 {
0795   for (int i = 0; i < 15; i++)  // last one is reserved for analogue sum
0796   {
0797     // multiply SMD channels with their gain factors
0798     // to get the absolute ADC values in the same units
0799     //rgains come from CompSmdAdc()
0800     smd_adc[i] = smd_adc[i] * smd_north_rgain[i];            // sout -> north for PHENIX -> sPHENIX
0801     smd_adc[i + 16] = smd_adc[i + 16] * smd_south_rgain[i];  // north -> south for PHENIX-> sPHENIX
0802   }
0803 }
0804 
0805 void ZdcMon::CompSmdPos()  //computing position with weighted averages
0806 {
0807   float w_ave[4];  // 0 -> north hor; 1 -> noth vert; 2 -> south hor; 3 -> south vert.
0808   float weights[4] = {0};
0809   memset(weights, 0, sizeof(weights));  // memset float works only for 0
0810   float w_sum[4];
0811   memset(w_sum, 0, sizeof(w_sum));
0812 
0813   // these constants convert the SMD channel number into real dimensions (cm's)
0814   const float hor_scale = 2.0 * 11.0 / 10.5 * sin(PI / 4);  // from gsl_math.h
0815   const float ver_scale = 1.5 * 11.0 / 10.5;
0816   float hor_offset = (hor_scale * 8 / 2.0) * (7.0 / 8.0);
0817   float ver_offset = (ver_scale * 7 / 2.0) * (6.0 / 7.0);
0818 
0819   for (int i = 0; i < 8; i++)
0820   {
0821     weights[0] += smd_adc[i];  //summing weights
0822     weights[2] += smd_adc[i + 16];
0823     w_sum[0] += (float) i * smd_adc[i];  //summing for the average
0824     w_sum[2] += ((float) i + 16.) * smd_adc[i + 16];
0825   }
0826   for (int i = 0; i < 7; i++)
0827   {
0828     weights[1] += smd_adc[i + 8];
0829     weights[3] += smd_adc[i + 24];
0830     w_sum[1] += ((float) i + 8.) * smd_adc[i + 8];
0831     w_sum[3] += ((float) i + 24.) * smd_adc[i + 24];
0832   }
0833 
0834   if (weights[0] > 0.0)
0835   {
0836     w_ave[0] = w_sum[0] / weights[0];  //average = sum / sumn of weights...
0837     smd_pos[0] = hor_scale * w_ave[0] - hor_offset;
0838   }
0839   else
0840   {
0841     smd_pos[0] = 0;
0842   }
0843   if (weights[1] > 0.0)
0844   {
0845     w_ave[1] = w_sum[1] / weights[1];
0846     smd_pos[1] = ver_scale * (w_ave[1] - 8.0) - ver_offset;
0847   }
0848   else
0849   {
0850     smd_pos[1] = 0;
0851   }
0852 
0853   if (weights[2] > 0.0)
0854   {
0855     w_ave[2] = w_sum[2] / weights[2];
0856     smd_pos[2] = hor_scale * (w_ave[2] - 16.0) - hor_offset;
0857   }
0858   else
0859   {
0860     smd_pos[2] = 0;
0861   }
0862 
0863   if (weights[3] > 0.0)
0864   {
0865     w_ave[3] = w_sum[3] / weights[3];
0866     smd_pos[3] = ver_scale * (w_ave[3] - 24.0) - ver_offset;
0867   }
0868   else
0869   {
0870     smd_pos[3] = 0;
0871   }
0872 }
0873 
0874 
0875