analysis/ZDCNeutronLocPol/ZDCNeutronLocPol.cc

0001 #include "ZDCNeutronLocPol.h"
0002
0003 #include <calobase/TowerInfoDefs.h>
0004 #include <caloreco/CaloWaveformFitting.h>
0005
0006 #include <calobase/TowerInfoDefs.h>
0007 #include <caloreco/CaloWaveformFitting.h>
0008 #include <ffarawobjects/CaloPacketContainerv1.h>
0009 #include <ffarawobjects/CaloPacketv1.h>
0010 #include <ffarawobjects/Gl1Packetv3.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <Event/packet.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>    // for PHIODataNode
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0019 #include <phool/PHObject.h>        // for PHObject
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022 #include <phool/recoConsts.h>
0023
0024 #include <TFile.h>
0025 #include <TH1.h>
0026 #include <TH2.h>
0027 #include <TTree.h>
0028 #include <phool/PHCompositeNode.h>
0029 #include <cmath>
0030 #include <fstream>
0031
0032 R__LOAD_LIBRARY(libuspin.so)
0033
0034 //____________________________________________________________________________..
0035 ZDCNeutronLocPol::ZDCNeutronLocPol(const std::string &name)
0036   : SubsysReco(name)
0037 {
0038   std::cout << "ZDCNeutronLocPol::ZDCNeutronLocPol(const std::string &name) Calling ctor" << std::endl;
0039 }
0040
0041 //____________________________________________________________________________..
0042 ZDCNeutronLocPol::~ZDCNeutronLocPol()
0043 {
0044   std::cout << "ZDCNeutronLocPol::~ZDCNeutronLocPol() Calling dtor" << std::endl;
0045 }
0046
0047 //____________________________________________________________________________..
0048 int ZDCNeutronLocPol::Init(PHCompositeNode * /*topNode*/)
0049 {
0050   std::cout << "ZDCNeutronLocPol::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0051
0052   if (!WaveformProcessingFast)
0053   {
0054     WaveformProcessingFast = new CaloWaveformFitting();
0055   }
0056
0057   smdHits = new TTree();
0058   smdHits = new TTree("smdHits", "smdHits");
0059   smdHits->SetDirectory(0);
0060   smdHits->Branch("bunchnumber", &bunchnumber, "bunchnumber/I");
0061   ///// Branches for in case that we might have mis-aligned GL1 and ZDC BCO in raw data.
0062   ///// Rarly happened at beginning of 24 pp runs, but just in case we store BCOs for debugging(if needed)
0063   // smdHits -> Branch("evtseq_gl1",  &evtseq_gl1, "evtseq_gl1/I");
0064   // smdHits -> Branch("evtseq_zdc",  &evtseq_zdc, "evtseq_zdc/I");
0065   // smdHits -> Branch("BCO_gl1",  &BCO_gl1, "BCO_gl1/l");
0066   // smdHits -> Branch("BCO_zdc",  &BCO_zdc, "BCO_zdc/l");
0067   smdHits->Branch("showerCutN", &showerCutN, "showerCutN/I");
0068   smdHits->Branch("showerCutS", &showerCutS, "showerCutS/I");
0069   smdHits->Branch("n_x", &n_x, "n_x/F");
0070   smdHits->Branch("n_y", &n_y, "n_y/F");
0071   smdHits->Branch("s_x", &s_x, "s_x/F");
0072   smdHits->Branch("s_y", &s_y, "s_y/F");
0073   smdHits->Branch("zdcN1_adc", &zdcN1_adc, "zdcN1_adc/F");
0074   smdHits->Branch("zdcN2_adc", &zdcN2_adc, "zdcN2_adc/F");
0075   smdHits->Branch("zdcN3_adc", &zdcN3_adc, "zdcN3_adc/F");
0076   smdHits->Branch("zdcS1_adc", &zdcS1_adc, "zdcS1_adc/F");
0077   smdHits->Branch("zdcS2_adc", &zdcS2_adc, "zdcS2_adc/F");
0078   smdHits->Branch("zdcS3_adc", &zdcS3_adc, "zdcS3_adc/F");
0079   smdHits->Branch("veto_NF", &veto_NF, "veto_NF/F");
0080   smdHits->Branch("veto_NB", &veto_NB, "veto_NB/F");
0081   smdHits->Branch("veto_SF", &veto_SF, "veto_SF/F");
0082   smdHits->Branch("veto_SB", &veto_SB, "veto_SB/F");
0083   smdHits->Branch("smdS1_adc", &smdS1_adc, "smdS1_adc/F");
0084   smdHits->Branch("smdS2_adc", &smdS2_adc, "smdS2_adc/F");
0085   smdHits->Branch("smdS6_adc", &smdS6_adc, "smdS6_adc/F");
0086   smdHits->Branch("smdS7_adc", &smdS7_adc, "smdS7_adc/F");
0087
0088   smdHits->Branch("smdN1_adc", &smdN1_adc, "smdN1_adc/F");
0089   smdHits->Branch("smdN2_adc", &smdN2_adc, "smdN2_adc/F");
0090   smdHits->Branch("smdN6_adc", &smdN6_adc, "smdN6_adc/F");
0091   smdHits->Branch("smdN7_adc", &smdN7_adc, "smdN7_adc/F");
0092
0093   smdHits->Branch("smdS1_v_adc", &smdS1_v_adc, "smdS1_v_adc/F");
0094   smdHits->Branch("smdS2_v_adc", &smdS2_v_adc, "smdS2_v_adc/F");
0095   smdHits->Branch("smdS7_v_adc", &smdS7_v_adc, "smdS7_v_adc/F");
0096   smdHits->Branch("smdS8_v_adc", &smdS8_v_adc, "smdS8_v_adc/F");
0097
0098   smdHits->Branch("smdN1_v_adc", &smdN1_v_adc, "smdN1_v_adc/F");
0099   smdHits->Branch("smdN2_v_adc", &smdN2_v_adc, "smdN2_v_adc/F");
0100   smdHits->Branch("smdN7_v_adc", &smdN7_v_adc, "smdN7_v_adc/F");
0101   smdHits->Branch("smdN8_v_adc", &smdN8_v_adc, "smdN8_v_adc/F");
0102
0103   // waveform
0104   h_waveformZDC = new TH2F("h_waveformZDC", "h_waveformZDC", 17, -0.5, 16.5, 512, -500, 20000);
0105   h_waveformSMD_North = new TH2F("h_waveformSMD_North", "h_waveformSMD_North", 17, -0.5, 16.5, 512, -500, 20000);
0106   h_waveformSMD_South = new TH2F("h_waveformSMD_South", "h_waveformSMD_South", 17, -0.5, 16.5, 512, -500, 20000);
0107   h_waveformVeto_North = new TH2F("h_waveformVeto_North", "h_waveformVeto_North", 17, -0.5, 16.5, 512, -500, 20000);
0108   h_waveformVeto_South = new TH2F("h_waveformVeto_South", "h_waveformVeto_South", 17, -0.5, 16.5, 512, -500, 20000);
0109   h_waveformAll = new TH2F("h_waveformAll", "h_waveformAll", 17, -0.5, 16.5, 512, -500, 20000);
0110
0111   for (int i = 0; i < 32; i++)
0112   {
0113     h_rawADC[i] = new TH1F(Form("rawADC%d", i), Form("rawADC%d", i), 128, 0, 4000);
0114     h_pedADC[i] = new TH1F(Form("pedADC%d", i), Form("pedADC%d", i), 128, 0, 4000);
0115     h_signalADC[i] = new TH1F(Form("signalDC%d", i), Form("signalADC%d", i), 128, 0, 4000);
0116   }
0117
0118   return Fun4AllReturnCodes::EVENT_OK;
0119 }
0120
0121 //____________________________________________________________________________..
0122 int ZDCNeutronLocPol::InitRun(PHCompositeNode * /*topNode*/)
0123 {
0124   std::cout << "ZDCNeutronLocPol::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0125
0126   evtcnt = 0;
0127
0128   return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130
0131 //____________________________________________________________________________..
0132 int ZDCNeutronLocPol::process_event(PHCompositeNode *topNode)
0133 {
0134   // std::cout << "ZDCNeutronLocPol::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0135   if (evtcnt % 10000 == 0)
0136   {
0137     std::cout << "Event: " << evtcnt << std::endl;
0138   }
0139
0140   p_gl1 = findNode::getClass<Gl1Packetv3>(topNode, 14001);
0141   p_zdc = findNode::getClass<CaloPacket>(topNode, 12001);
0142
0143   showerCutN = 0;
0144   showerCutS = 0;
0145
0146   if (p_gl1)
0147   {
0148     bunchnumber = p_gl1->getBunchNumber();
0149     if (!p_zdc)
0150     {
0151       return Fun4AllReturnCodes::EVENT_OK;
0152     }
0153     if (p_zdc)
0154     {
0155       int nThresholdNVer = 0;
0156       int nThresholdNHor = 0;
0157       int nThresholdSVer = 0;
0158       int nThresholdSHor = 0;
0159
0160       // evtseq_gl1 = p_gl1->getEvtSequence();
0161       // evtseq_zdc = p_zdc->getEvtSequence();
0162       // BCO_gl1 = p_gl1->getBCO();
0163       // BCO_zdc = p_zdc->getBCO();
0164
0165       // in this for loop we get: zdc_adc and smd_adc
0166       for (int c = 0; c < p_zdc->iValue(0, "CHANNELS"); c++)
0167       {
0168         std::vector<float> resultFast = anaWaveformFast(p_zdc, c);  // fast waveform fitting
0169         float signalFast = resultFast.at(0);
0170         float timingFast = resultFast.at(1);
0171         float pedFast = resultFast.at(2);
0172
0173         ////////// ZDC channels ///////
0174         if (c >= 0 && c < 16)
0175         {
0176           signalFast = (4 < timingFast && timingFast < 10) ? signalFast : 0;
0177           unsigned int towerkey = TowerInfoDefs::decode_zdc(c);
0178           int zdc_side = TowerInfoDefs::get_zdc_side(towerkey);
0179           if (zdc_side == 0)
0180           {
0181             if (c == 0)
0182             {
0183               zdcS1_adc = signalFast;
0184               continue;
0185             }
0186             if (c == 2)
0187             {
0188               zdcS2_adc = signalFast;
0189               continue;
0190             }
0191             if (c == 4)
0192             {
0193               zdcS3_adc = signalFast;
0194               continue;
0195             }
0196           }
0197           else if (zdc_side == 1)
0198           {
0199             if (c == 8)
0200             {
0201               zdcN1_adc = signalFast;
0202               continue;
0203             }
0204             if (c == 10)
0205             {
0206               zdcN2_adc = signalFast;
0207               continue;
0208             }
0209             if (c == 12)
0210             {
0211               zdcN3_adc = signalFast;
0212               continue;
0213             }
0214           }
0215         }
0216         /////////////end ZDC channels//////////////////
0217
0218         //////////// North Veto counters //////////////////
0219         if (c == 16)
0220         {
0221           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0222           veto_NF = signalFast;
0223           continue;
0224         }
0225         if (c == 17)
0226         {
0227           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0228           veto_NB = signalFast;
0229           continue;
0230         }
0231         //////////// end North veto counters /////////
0232
0233         ///////// North SMD ////////
0234         if (c >= 48 && c < 63)
0235         {
0236           signalFast = (9 < timingFast) ? signalFast : 0;
0237
0238           rawADC[c - 48] = signalFast + pedFast;
0239           pedADC[c - 48] = pedFast;
0240           signalADC[c - 48] = signalFast;
0241
0242           h_rawADC[c - 48]->Fill(signalFast + pedFast);
0243           h_pedADC[c - 48]->Fill(pedFast);
0244           h_signalADC[c - 48]->Fill(signalFast);
0245           if (c == 48)
0246           {
0247             smdN1_v_adc = signalFast;
0248           }
0249           if (c == 49)
0250           {
0251             smdN2_v_adc = signalFast;
0252           }
0253           if (c == 54)
0254           {
0255             smdN7_v_adc = signalFast;
0256           }
0257           if (c == 55)
0258           {
0259             smdN8_v_adc = signalFast;
0260           }
0261           if (c == 56)
0262           {
0263             smdN1_adc = signalFast;
0264           }
0265           if (c == 57)
0266           {
0267             smdN2_adc = signalFast;
0268           }
0269           if (c == 61)
0270           {
0271             smdN6_adc = signalFast;
0272           }
0273           if (c == 62)
0274           {
0275             smdN7_adc = signalFast;
0276           }
0277           if (c < 56 && signalFast > 5.)
0278           {
0279             nThresholdNHor++;
0280           }
0281           if (c >= 56 && signalFast > 5.)
0282           {
0283             nThresholdNVer++;
0284           }
0285           continue;
0286         }
0287         //////// end North SMD///////////////
0288
0289         //////////// South Veto counters //////////////////
0290         if (c == 80)
0291         {
0292           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0293           veto_SF = signalFast;
0294           continue;
0295         }
0296         if (c == 81)
0297         {
0298           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0299           veto_SB = signalFast;
0300           continue;
0301         }
0302         //////////// end South veto counters /////////
0303         //// South SMD ////////
0304         if (c >= 112 && c < 127)
0305         {
0306           signalFast = (7 < timingFast && timingFast < 12) ? signalFast : 0;
0307           rawADC[c - 96] = signalFast + pedFast;
0308           pedADC[c - 96] = pedFast;
0309           signalADC[c - 96] = signalFast;
0310
0311           h_rawADC[c - 96]->Fill(signalFast + pedFast);
0312           h_pedADC[c - 96]->Fill(pedFast);
0313           h_signalADC[c - 96]->Fill(signalFast);
0314           if (c == 112)
0315           {
0316             smdS1_v_adc = signalFast;
0317           }
0318           if (c == 113)
0319           {
0320             smdS2_v_adc = signalFast;
0321           }
0322           if (c == 118)
0323           {
0324             smdS7_v_adc = signalFast;
0325           }
0326           if (c == 119)
0327           {
0328             smdS8_v_adc = signalFast;
0329           }
0330           if (c == 120)
0331           {
0332             smdS1_adc = signalFast;
0333           }
0334           if (c == 121)
0335           {
0336             smdS2_adc = signalFast;
0337           }
0338           if (c == 125)
0339           {
0340             smdS6_adc = signalFast;
0341           }
0342           if (c == 126)
0343           {
0344             smdS7_adc = signalFast;
0345           }
0346           if (c < 120 && signalFast > 5.)
0347           {
0348             nThresholdSHor++;
0349           }
0350           if (c >= 120 && signalFast > 5.)
0351           {
0352             nThresholdSVer++;
0353           }
0354           continue;
0355         }
0356         /////// end South SMD //////////////
0357
0358       }  // end channel loop
0359       std::cout<<"Hello 3"<<std::endl;
0360       // shower cut (SMD noise threshold cut)
0361       if (nThresholdNHor > 1 && nThresholdNVer > 1)
0362       {
0363         showerCutN = 1;
0364       }
0365
0366       if (nThresholdSHor > 1 && nThresholdSVer > 1)
0367       {
0368         showerCutS = 1;
0369       }
0370
0371       ////////// Compute SMD position and fill hists
0372       CompSmdAdc();
0373       CompSmdPos();
0374       n_y = smd_pos[0];
0375       n_x = smd_pos[1];
0376       s_y = smd_pos[2];
0377
0378       s_x = smd_pos[3];
0379       smdHits->Fill();
0380     }  // end if p_zdc good
0381   }  // end if p_gl1 good
0382   else
0383   {
0384     std::cout<<"gl1 packet not found"<<std::endl;
0385   }
0386   evtcnt++;
0387
0388   return Fun4AllReturnCodes::EVENT_OK;
0389 }
0390
0391 //____________________________________________________________________________..
0392 int ZDCNeutronLocPol::End(PHCompositeNode * /*topNode*/)
0393 {
0394   std::cout << "ZDCNeutronLocPol::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0395   TFile *ofile = new TFile(outfile.c_str(), "RECREATE");
0396
0397   smdHits->Write();
0398
0399   // h_waveformZDC->Write();
0400   // h_waveformSMD_North->Write();
0401   // h_waveformSMD_South->Write();
0402   // h_waveformVeto_North->Write();
0403   // h_waveformVeto_South->Write();
0404   // h_waveformAll->Write();
0405   /*
0406   for (int i = 0; i < 32; i++)
0407   {
0408     h_rawADC[i]->Write();
0409     h_pedADC[i]->Write();
0410     h_signalADC[i]->Write();
0411   }
0412   */
0413   ofile->Write();
0414   ofile->Close();
0415   delete (smdHits);
0416
0417   return Fun4AllReturnCodes::EVENT_OK;
0418 }
0419
0420 std::vector<float> ZDCNeutronLocPol::anaWaveformFast(CaloPacket *p, const int channel)
0421 {
0422   std::vector<float> waveform;
0423   // Chris: preallocation = speed improvement
0424   waveform.reserve(p->iValue(0, "SAMPLES"));
0425   for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0426   {
0427     waveform.push_back(p->iValue(s, channel));
0428   }
0429   std::vector<std::vector<float>> multiple_wfs;
0430   multiple_wfs.push_back(waveform);
0431
0432   std::vector<std::vector<float>> fitresults_zdc;
0433   fitresults_zdc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
0434
0435   std::vector<float> result;
0436   result = fitresults_zdc.at(0);
0437   return result;
0438 }
0439
0440 void ZDCNeutronLocPol::CompSmdAdc()  // mulitplying by relative gains
0441 {
0442   for (int i = 0; i < 15; i++)  // last one is reserved for analog sum
0443   {
0444     // multiply SMD channels with their gain factors
0445     // to get the absolute ADC values in the same units
0446     // rgains come from CompSmdAdc()
0447     signalADC[i] = signalADC[i] * smd_north_rgain[i];
0448     signalADC[i + 16] = signalADC[i + 16] * smd_south_rgain[i];
0449   }
0450 }
0451
0452 void ZDCNeutronLocPol::CompSmdPos()  // computing position with weighted averages
0453 {
0454   float w_ave[4];  // 0 -> north hor (y); 1 -> noth vert (x); 2 -> south hor (y); 3 -> south vert. (x)
0455   float weights[4] = {0};
0456   memset(weights, 0, sizeof(weights));  // memset float works only for 0
0457   float w_sum[4];
0458   memset(w_sum, 0, sizeof(w_sum));
0459
0460   float smd_threshold = 5.;
0461
0462   // these constants convert the SMD channel number into real dimensions (cm's)
0463   const float hor_scale = 2.0 * 11.0 / 10.5 * sin(3.1415 / 4);  // from gsl_math.h
0464   const float ver_scale = 1.5 * 11.0 / 10.5;
0465   float hor_offset = (hor_scale * 8 / 2.0) * (7.0 / 8.0);
0466   float ver_offset = (ver_scale * 7 / 2.0) * (6.0 / 7.0);
0467
0468   for (int i = 0; i < 8; i++)
0469   {
0470     if (signalADC[i] > smd_threshold)
0471     {
0472       weights[0] += signalADC[i];            // summing weights
0473       w_sum[0] += (float) i * signalADC[i];  // summing for the average
0474     }
0475
0476     if (signalADC[i + 16] > smd_threshold)
0477     {
0478       weights[2] += signalADC[i + 16];
0479       w_sum[2] += ((float) i + 16.) * signalADC[i + 16];
0480     }
0481   }
0482   for (int i = 0; i < 7; i++)
0483   {
0484     if (signalADC[i + 8] > smd_threshold)
0485     {
0486       weights[1] += signalADC[i + 8];
0487       w_sum[1] += ((float) i + 8.) * signalADC[i + 8];
0488     }
0489
0490     if (signalADC[i + 24] > smd_threshold)
0491     {
0492       weights[3] += signalADC[i + 24];
0493       w_sum[3] += ((float) i + 24.) * signalADC[i + 24];
0494     }
0495   }
0496
0497   if (weights[0] > 0.0)
0498   {
0499     w_ave[0] = w_sum[0] / weights[0];  // average = sum / sumn of weights...
0500     smd_pos[0] = hor_scale * w_ave[0] - hor_offset;
0501   }
0502   else
0503   {
0504     smd_pos[0] = 0;
0505   }
0506
0507   if (weights[1] > 0.0)
0508   {
0509     w_ave[1] = w_sum[1] / weights[1];
0510     smd_pos[1] = ver_scale * (w_ave[1] - 8.0) - ver_offset;
0511   }
0512   else
0513   {
0514     smd_pos[1] = 0;
0515   }
0516
0517   if (weights[2] > 0.0)
0518   {
0519     w_ave[2] = w_sum[2] / weights[2];
0520     smd_pos[2] = hor_scale * (w_ave[2] - 16.0) - hor_offset;
0521   }
0522   else
0523   {
0524     smd_pos[2] = 0;
0525   }
0526
0527   if (weights[3] > 0.0)
0528   {
0529     w_ave[3] = w_sum[3] / weights[3];
0530     smd_pos[3] = ver_scale * (w_ave[3] - 24.0) - ver_offset;
0531   }
0532   else
0533   {
0534     smd_pos[3] = 0;
0535   }
0536 }
0537
0538 void ZDCNeutronLocPol::setFileName(const std::string &fname)
0539 {
0540   outfile = fname;
0541 }
0542
0543 void ZDCNeutronLocPol::setGainMatch(const std::string &gainfile)
0544 {
0545   // getting gains
0546   float col1;
0547   std::string col2;
0548    //std::string gainfile = gainname;
0549    std::ifstream gain_infile(gainfile);
0550
0551   if (!gain_infile)
0552   {
0553     std::cout << "No relative gain file: All relative gains set to 1." << std::endl;
0554     return;
0555   }
0556
0557   std::cout << "SMD Relative Gains: " << std::endl;
0558
0559   for (float &i : gain)
0560   {
0561     gain_infile >> col1 >> col2;
0562     i = col1;
0563   }
0564
0565   for (int i = 0; i < 16; i++)  // relative gains of SMD north channels
0566   {
0567     smd_north_rgain[i] = gain[i];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0568     std::cout << "N" << i << " " << gain[i] << std::endl;
0569   }
0570
0571   for (int i = 0; i < 16; i++)  // relative gains of SMD south channels
0572   {
0573     smd_south_rgain[i] = gain[i + 16];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0574     std::cout << "S" << i << " " << gain[i] << std::endl;
0575   }
0576
0577   gain_infile.close();
0578 }