QA/Calorimeters/CaloValid.cc

0001 #include "CaloValid.h"
0002
0003 // Calo includes
0004 #include <calobase/RawCluster.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawClusterUtility.h>
0007 #include <calobase/RawTowerGeomContainer.h>
0008 #include <calobase/TowerInfo.h>
0009 #include <calobase/TowerInfoContainer.h>
0010
0011 #include <mbd/MbdPmtContainer.h>
0012 #include <mbd/MbdPmtHit.h>
0013
0014 #include <globalvertex/GlobalVertex.h>
0015 #include <globalvertex/GlobalVertexMap.h>
0016
0017 #include <qautils/QAHistManagerDef.h>
0018
0019 #include <ffaobjects/EventHeader.h>
0020 #include <ffaobjects/RunHeader.h>  // for runnumber
0021
0022 #include <ffarawobjects/Gl1Packet.h>
0023
0024 #include <fun4all/Fun4AllHistoManager.h>
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026
0027 #include <phool/RunnumberRange.h>
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>  // for PHWHERE
0030
0031 #include <TH1.h>
0032 #include <TH2.h>
0033 #include <TH3.h>
0034 #include <TLorentzVector.h>
0035 #include <TProfile.h>
0036 #include <TProfile2D.h>
0037 #include <TSystem.h>
0038
0039 #include <cassert>
0040 #include <cmath>  // for log10, pow, sqrt, abs, M_PI
0041 #include <cstdint>
0042 #include <format>
0043 #include <iostream>  // for operator<<, endl, basic_...
0044 #include <limits>
0045 #include <map>  // for operator!=, _Rb_tree_con...
0046 #include <string>
0047 #include <utility>  // for pair
0048
0049 CaloValid::CaloValid(const std::string& name)
0050   : SubsysReco(name)
0051 {
0052 }
0053
0054 CaloValid::~CaloValid()
0055 {
0056   for (int i = 0; i < 128 * 192; i++)
0057   {
0058     delete h_cemc_channel_pedestal[i];
0059     delete h_cemc_channel_energy[i];
0060   }
0061   for (int i = 0; i < 32 * 48; i++)
0062   {
0063     delete h_ihcal_channel_pedestal[i];
0064     delete h_ihcal_channel_energy[i];
0065     delete h_ohcal_channel_pedestal[i];
0066     delete h_ohcal_channel_energy[i];
0067   }
0068   delete trigAna;
0069 }
0070
0071 int CaloValid::Init(PHCompositeNode* /*unused*/)
0072 {
0073   if (m_debug)
0074   {
0075     std::cout << "In CaloValid::Init" << std::endl;
0076   }
0077
0078   createHistos();
0079   trigAna = new TriggerAnalyzer();
0080   if (m_debug)
0081   {
0082     std::cout << "Leaving CaloValid::Init" << std::endl;
0083   }
0084   return Fun4AllReturnCodes::EVENT_OK;
0085 }
0086
0087 int CaloValid::InitRun(PHCompositeNode* topNode)
0088 {
0089   RunHeader* runhdr = findNode::getClass<RunHeader>(topNode, "RunHeader");
0090
0091   if (runhdr)
0092   {
0093     int runnumber = runhdr->get_RunNumber();
0094
0095     if (runnumber >= RunnumberRange::RUN2PP_FIRST && runnumber <= RunnumberRange::RUN2PP_LAST)
0096     {
0097       m_species = "pp";
0098       if (Verbosity() > 0)
0099       {
0100         std::cout << "This run is from Run-2 p+p.\n";
0101       }
0102     }
0103     else if (runnumber >= RunnumberRange::RUN2AUAU_FIRST && runnumber <= RunnumberRange::RUN2AUAU_LAST)
0104     {
0105       m_species = "AuAu";
0106       if (Verbosity() > 0)
0107       {
0108         std::cout << "This run is from Run-2 Au+Au.\n";
0109       }
0110     }
0111     else if (runnumber >= RunnumberRange::RUN3AUAU_FIRST && runnumber <= RunnumberRange::RUN3AUAU_LAST)
0112     {
0113       m_species = "AuAu";
0114       if (Verbosity() > 0)
0115       {
0116         std::cout << "This run is from Run-3 Au+Au.\n";
0117       }
0118     }
0119     else
0120     {
0121       if (Verbosity() > 0)
0122       {
0123         std::cout << "Run number is out of range. Check RunnumberRange.h .Using pp as default. \n";
0124       }
0125     }
0126   }
0127   else if (Verbosity() > 0)
0128   {
0129     std::cout << "No RunHeader node found. Using pp as default species.\n";
0130   }
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 int CaloValid::process_event(PHCompositeNode* topNode)
0136 {
0137   _eventcounter++;
0138   //  std::cout << "In process_event" << std::endl;
0139   process_towers(topNode);
0140
0141   return Fun4AllReturnCodes::EVENT_OK;
0142 }
0143
0144 int CaloValid::process_towers(PHCompositeNode* topNode)
0145 {
0146   //---------------------------Event header--------------------------------//
0147   EventHeader* eventheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0148   int event_number = 0;
0149   if (eventheader)
0150   {
0151     if (eventheader->isValid())
0152     {
0153       event_number = eventheader->get_EvtSequence();
0154     }
0155   }
0156   else
0157   {
0158     std::cout << "GlobalQA::process_event()  No event header" << std::endl;
0159   }
0160
0161   if (m_debug)
0162   {
0163     std::cout << _eventcounter << std::endl;
0164   }
0165   //  std::cout << "In process_towers" << std::endl;
0166   auto* hm = QAHistManagerDef::getHistoManager();
0167   assert(hm);
0168
0169   float totalcemc = 0.;
0170   float totalihcal = 0.;
0171   float totalohcal = 0.;
0172   float totalmbd = 0.;
0173
0174   float emcaldownscale;
0175   float ihcaldownscale;
0176   float ohcaldownscale;
0177   float mbddownscale;
0178   float adc_threshold;
0179   float emcal_hit_threshold;
0180   float ohcal_hit_threshold;
0181   float ihcal_hit_threshold;
0182
0183   if (m_species == "AuAu")
0184   {
0185     emcaldownscale = 1350000. / 800.;
0186     ihcaldownscale = 55000. / 300.;
0187     ohcaldownscale = 265000. / 600.;
0188     mbddownscale = 2800.0;
0189     adc_threshold = 15.;
0190
0191     emcal_hit_threshold = 0.5;  // GeV
0192     ohcal_hit_threshold = 0.5;
0193     ihcal_hit_threshold = 0.25;
0194   }
0195   else
0196   {
0197     emcaldownscale = 100000. / 800.;
0198     ihcaldownscale = 4000. / 300.;
0199     ohcaldownscale = 25000. / 600.;
0200     mbddownscale = 200.0;
0201     adc_threshold = 100.;
0202
0203     emcal_hit_threshold = 0.5;  // GeV
0204     ohcal_hit_threshold = 0.5;
0205     ihcal_hit_threshold = 0.25;
0206   }
0207
0208   //----------------------------------vertex------------------------------------------------------//
0209   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0210   if (!vertexmap)
0211   {
0212     std::cout << "CaloValid GlobalVertexMap node is missing" << std::endl;
0213   }
0214   float vtx_z = std::numeric_limits<float>::quiet_NaN();
0215
0216   if (vertexmap && !vertexmap->empty())
0217   {
0218     GlobalVertex* vtx = vertexmap->begin()->second;
0219     if (vtx)
0220     {
0221       vtx_z = vtx->get_z();
0222     }
0223     h_vtx_z_raw->Fill(vtx_z);
0224   }
0225
0226   //--------------------------- trigger and GL1-------------------------------//
0227   if (trigAna)
0228   {
0229     trigAna->decodeTriggers(topNode);
0230   }
0231   else
0232   {
0233     if (m_debug)
0234     {
0235       std::cout << "[ERROR] No TriggerAnalyzer pointer!\n";
0236     }
0237   }
0238
0239   std::vector<int> scaledActiveBits;
0240   scaledActiveBits.reserve(triggerIndices.size());
0241
0242   for (int bit : triggerIndices)
0243   {
0244     if (trigAna->didTriggerFire(bit))
0245     {
0246       scaledActiveBits.push_back(bit);
0247     }
0248   }
0249
0250   bool scaledBits[64] = {false};
0251   uint64_t raw[64] = {0};
0252   uint64_t live[64] = {0};
0253   // long long int scaled[64] = { 0 };
0254   Gl1Packet* gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, 14001);
0255   if (!gl1PacketInfo)
0256   {
0257     gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0258     if (!gl1PacketInfo)
0259     {
0260       std::cout << PHWHERE << "GlobalQA::process_event: GL1Packet node is missing"
0261                 << std::endl;
0262     }
0263   }
0264   uint64_t triggervec = 0;
0265   if (gl1PacketInfo)
0266   {
0267     triggervec = gl1PacketInfo->getScaledVector();
0268     for (int i = 0; i < 64; i++)
0269     {
0270       bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0271       scaledBits[i] = trig_decision;
0272
0273       raw[i] = gl1PacketInfo->lValue(i, 0);
0274       live[i] = gl1PacketInfo->lValue(i, 1);
0275       // scaled[i] = gl1PacketInfo->lValue(i, 2);
0276
0277       if (trig_decision)
0278       {
0279         h_triggerVec->Fill(i);
0280       }
0281       triggervec = (triggervec >> 1U) & 0xffffffffU;
0282     }
0283     // triggervec = gl1PacketInfo->getScaledVector(); commented out to get rid of never used warning using clang -tidy
0284   }
0285
0286   //---------------------------calibrated towers-------------------------------//
0287   {
0288     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0289     if (towers)
0290     {
0291       int size = towers->size();  // online towers should be the same!
0292       for (int channel = 0; channel < size; channel++)
0293       {
0294         TowerInfo* tower = towers->get_tower_at_channel(channel);
0295         float offlineenergy = tower->get_energy();
0296         unsigned int towerkey = towers->encode_key(channel);
0297         int ieta = towers->getTowerEtaBin(towerkey);
0298         int iphi = towers->getTowerPhiBin(towerkey);
0299         int _time = tower->get_time();
0300         h_cemc_e_chi2->Fill(offlineenergy, tower->get_chi2());
0301         float _timef = tower->get_time_float();
0302         h_emcaltime_cut->Fill(_time);
0303         bool isGood = tower->get_isGood();
0304         uint8_t status = tower->get_status();
0305         h_emcal_tower_e->Fill(offlineenergy);
0306         h_emcal_tower_e_wide_range->Fill(offlineenergy);
0307         if (tower->get_isSaturated())
0308         {
0309           h_emcal_tower_e_saturated->Fill(offlineenergy);
0310         }
0311         float pedestal = tower->get_pedestal();
0312         h_cemc_channel_pedestal[channel]->Fill(pedestal);
0313
0314         for (int is = 0; is < 8; is++)
0315         {
0316           if (status & 1U)  // clang-tidy mark 1 as unsigned
0317           {
0318             h_cemc_status->Fill(is);
0319           }
0320           status = status >> 1U;  // clang-tidy mark 1 as unsigned
0321         }
0322
0323         totalcemc += offlineenergy;
0324         h_emcaltime->Fill(_time);
0325         if (offlineenergy > emcal_hit_threshold)
0326         {
0327           h_cemc_etaphi_time->Fill(ieta, iphi, _timef);
0328           h_cemc_etaphi->Fill(ieta, iphi);
0329           if (isGood && (scaledBits[10] || scaledBits[11]))
0330           {
0331             h_cemc_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0332           }
0333           if (tower->get_isBadChi2())
0334           {
0335             h_cemc_etaphi_badChi2->Fill(ieta, iphi, 1);
0336           }
0337           else
0338           {
0339             h_cemc_etaphi_badChi2->Fill(ieta, iphi, 0);
0340           }
0341         }
0342         if (offlineenergy > 0.25)
0343         {
0344           h_cemc_etaphi_fracHit->Fill(ieta, iphi, 1);
0345         }
0346         else
0347         {
0348           h_cemc_etaphi_fracHit->Fill(ieta, iphi, 0);
0349         }
0350       }
0351     }
0352   }
0353   {
0354     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0355     if (towers)
0356     {
0357       int size = towers->size();  // online towers should be the same!
0358       for (int channel = 0; channel < size; channel++)
0359       {
0360         TowerInfo* tower = towers->get_tower_at_channel(channel);
0361         float offlineenergy = tower->get_energy();
0362         unsigned int towerkey = towers->encode_key(channel);
0363         int ieta = towers->getTowerEtaBin(towerkey);
0364         int iphi = towers->getTowerPhiBin(towerkey);
0365         int _time = tower->get_time();
0366         float _timef = tower->get_time_float();
0367         h_ihcaltime_cut->Fill(_time);
0368         h_ihcal_e_chi2->Fill(offlineenergy, tower->get_chi2());
0369         bool isGood = tower->get_isGood();
0370         h_ihcal_status->Fill(tower->get_status());
0371         float pedestal = tower->get_pedestal();
0372         h_ihcal_channel_pedestal[channel]->Fill(pedestal);
0373         h_ihcal_tower_e->Fill(offlineenergy);
0374         h_ihcal_tower_e_wide_range->Fill(offlineenergy);
0375         if (tower->get_isSaturated())
0376         {
0377           h_ihcal_tower_e_saturated->Fill(offlineenergy);
0378         }
0379
0380         uint8_t status = tower->get_status();
0381         for (int is = 0; is < 8; is++)
0382         {
0383           if (status & 1U)  // clang-tidy mark 1 as unsigned
0384           {
0385             h_ihcal_status->Fill(is);
0386           }
0387           status = status >> 1U;  // clang-tidy mark 1 as unsigned
0388         }
0389
0390         totalihcal += offlineenergy;
0391         h_ihcaltime->Fill(_time);
0392
0393         if (offlineenergy > ihcal_hit_threshold)
0394         {
0395           h_ihcal_etaphi->Fill(ieta, iphi);
0396           h_ihcal_etaphi_time->Fill(ieta, iphi, _timef);
0397           if (isGood && (scaledBits[10] || scaledBits[11]))
0398           {
0399             h_ihcal_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0400           }
0401           if (tower->get_isBadChi2())
0402           {
0403             h_ihcal_etaphi_badChi2->Fill(ieta, iphi, 1);
0404           }
0405           else
0406           {
0407             h_ihcal_etaphi_badChi2->Fill(ieta, iphi, 0);
0408           }
0409         }
0410       }
0411     }
0412   }
0413
0414   {
0415     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0416     if (towers)
0417     {
0418       int size = towers->size();  // online towers should be the same!
0419       for (int channel = 0; channel < size; channel++)
0420       {
0421         TowerInfo* tower = towers->get_tower_at_channel(channel);
0422         float offlineenergy = tower->get_energy();
0423         unsigned int towerkey = towers->encode_key(channel);
0424         int ieta = towers->getTowerEtaBin(towerkey);
0425         int iphi = towers->getTowerPhiBin(towerkey);
0426         int _time = tower->get_time();
0427         float _timef = tower->get_time_float();
0428         h_ohcaltime_cut->Fill(_time);
0429         h_ohcal_e_chi2->Fill(offlineenergy, tower->get_chi2());
0430         bool isGood = tower->get_isGood();
0431         h_ohcal_status->Fill(tower->get_status());
0432         float pedestal = tower->get_pedestal();
0433         h_ohcal_channel_pedestal[channel]->Fill(pedestal);
0434         h_ohcal_tower_e->Fill(offlineenergy);
0435         h_ohcal_tower_e_wide_range->Fill(offlineenergy);
0436         if (tower->get_isSaturated())
0437         {
0438           h_ohcal_tower_e_saturated->Fill(offlineenergy);
0439         }
0440
0441         uint8_t status = tower->get_status();
0442         for (int is = 0; is < 8; is++)
0443         {
0444           if (status & 1U)  // clang-tidy mark 1 as unsigned
0445           {
0446             h_ohcal_status->Fill(is);
0447           }
0448           status = status >> 1U;  // clang-tidy mark 1 as unsigned
0449         }
0450
0451         totalohcal += offlineenergy;
0452         h_ohcaltime->Fill(_time);
0453
0454         if (offlineenergy > ohcal_hit_threshold)
0455         {
0456           h_ohcal_etaphi->Fill(ieta, iphi);
0457           h_ohcal_etaphi_time->Fill(ieta, iphi, _timef);
0458           if (isGood && (scaledBits[10] || scaledBits[11]))
0459           {
0460             h_ohcal_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0461           }
0462           if (tower->get_isBadChi2())
0463           {
0464             h_ohcal_etaphi_badChi2->Fill(ieta, iphi, 1);
0465           }
0466           else
0467           {
0468             h_ohcal_etaphi_badChi2->Fill(ieta, iphi, 0);
0469           }
0470         }
0471       }
0472     }
0473   }
0474
0475   //-------------------------- raw tower ------------------------------//
0476
0477   {
0478     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_CEMC");
0479     if (towers)
0480     {
0481       int size = towers->size();  // online towers should be the same!
0482       for (int channel = 0; channel < size; channel++)
0483       {
0484         TowerInfo* tower = towers->get_tower_at_channel(channel);
0485         unsigned int towerkey = towers->encode_key(channel);
0486         int ieta = towers->getTowerEtaBin(towerkey);
0487         int iphi = towers->getTowerPhiBin(towerkey);
0488         float raw_time = tower->get_time_float();
0489         if (tower->get_isZS())
0490         {
0491           h_cemc_channel_energy[channel]->Fill(tower->get_energy());
0492         }
0493
0494         float raw_energy = tower->get_energy();
0495         if (raw_energy > adc_threshold)
0496         {
0497           h_cemc_etaphi_fracHitADC->Fill(ieta, iphi, 1);
0498           h_cemc_etaphi_time_raw->Fill(ieta, iphi, raw_time);
0499         }
0500         else
0501         {
0502           h_cemc_etaphi_fracHitADC->Fill(ieta, iphi, 0);
0503         }
0504       }
0505     }
0506   }
0507   {
0508     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALOUT");
0509     if (towers)
0510     {
0511       int size = towers->size();  // online towers should be the same!
0512       for (int channel = 0; channel < size; channel++)
0513       {
0514         TowerInfo* tower = towers->get_tower_at_channel(channel);
0515         unsigned int towerkey = towers->encode_key(channel);
0516         int ieta = towers->getTowerEtaBin(towerkey);
0517         int iphi = towers->getTowerPhiBin(towerkey);
0518         float raw_time = tower->get_time_float();
0519         if (tower->get_isZS())
0520         {
0521           h_ohcal_channel_energy[channel]->Fill(tower->get_energy());
0522         }
0523
0524         float raw_energy = tower->get_energy();
0525         if (raw_energy > adc_threshold)
0526         {
0527           h_ohcal_etaphi_time_raw->Fill(ieta, iphi, raw_time);
0528           h_ohcal_etaphi_fracHitADC->Fill(ieta, iphi, 1);
0529         }
0530         else
0531         {
0532           h_ohcal_etaphi_fracHitADC->Fill(ieta, iphi, 0);
0533         }
0534       }
0535     }
0536   }
0537   {
0538     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALIN");
0539     if (towers)
0540     {
0541       int size = towers->size();  // online towers should be the same!
0542       for (int channel = 0; channel < size; channel++)
0543       {
0544         TowerInfo* tower = towers->get_tower_at_channel(channel);
0545         unsigned int towerkey = towers->encode_key(channel);
0546         int ieta = towers->getTowerEtaBin(towerkey);
0547         float raw_time = tower->get_time_float();
0548         int iphi = towers->getTowerPhiBin(towerkey);
0549         if (tower->get_isZS())
0550         {
0551           h_ihcal_channel_energy[channel]->Fill(tower->get_energy());
0552         }
0553
0554         float raw_energy = tower->get_energy();
0555         if (raw_energy > adc_threshold)
0556         {
0557           h_ihcal_etaphi_time_raw->Fill(ieta, iphi, raw_time);
0558           h_ihcal_etaphi_fracHitADC->Fill(ieta, iphi, 1);
0559         }
0560         else
0561         {
0562           h_ihcal_etaphi_fracHitADC->Fill(ieta, iphi, 0);
0563         }
0564       }
0565     }
0566   }
0567
0568   //--------------------------- MBD ----------------------------------------//
0569   MbdPmtContainer* bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0570   if (!bbcpmts)
0571   {
0572     std::cout << "CaloVald::process_event: Could not find MbdPmtContainer," << std::endl;
0573     // return Fun4AllReturnCodes::ABORTEVENT;
0574   }
0575
0576   int hits = 0;
0577   if (bbcpmts)
0578   {
0579     int nPMTs = bbcpmts->get_npmt();
0580     for (int i = 0; i < nPMTs; i++)
0581     {
0582       MbdPmtHit* mbdpmt = bbcpmts->get_pmt(i);
0583       float pmtadc = mbdpmt->get_q();
0584       totalmbd += pmtadc;
0585       if (pmtadc > 0.4)
0586       {
0587         hits++;
0588       }
0589     }
0590   }
0591   h_mbd_hits->Fill(hits);
0592
0593   h_emcal_mbd_correlation->Fill(totalcemc / emcaldownscale, totalmbd / mbddownscale);
0594   h_ihcal_mbd_correlation->Fill(totalihcal / ihcaldownscale, totalmbd / mbddownscale);
0595   h_ohcal_mbd_correlation->Fill(totalohcal / ohcaldownscale, totalmbd / mbddownscale);
0596   h_emcal_hcal_correlation->Fill(totalcemc / emcaldownscale, totalohcal / ohcaldownscale);
0597
0598   //------------------------------ clusters & pi0 ------------------------------//
0599   RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0600   if (!clusterContainer)
0601   {
0602     std::cout << PHWHERE << "CaloValid::funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0603     return 0;
0604   }
0605
0606   //////////////////////////////////////////
0607   // geometry for hot tower/cluster masking
0608   std::string towergeomnodename = "TOWERGEOM_CEMC";
0609   RawTowerGeomContainer* m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0610   if (!m_geometry)
0611   {
0612     std::cout << Name() << "::"
0613               << "CreateNodeTree"
0614               << ": Could not find node " << towergeomnodename << std::endl;
0615     gSystem->Exit(1);
0616   }
0617
0618   // cuts
0619   float emcMinClusE1 = 1.3;  // 0.5;
0620   float emcMinClusE2 = 0.7;  // 0.5;
0621   float emcMaxClusE = 100;
0622   float maxAlpha = 0.6;
0623
0624   if (totalcemc < 0.2 * emcaldownscale)
0625   {
0626     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0627     RawClusterContainer::ConstIterator clusterIter;
0628     RawClusterContainer::ConstIterator clusterIter2;
0629
0630     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0631     {
0632       RawCluster* recoCluster = clusterIter->second;
0633
0634       CLHEP::Hep3Vector vertex(0, 0, 0);
0635       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0636
0637       float clusE = E_vec_cluster.mag();
0638       float clus_eta = E_vec_cluster.pseudoRapidity();
0639       float clus_phi = E_vec_cluster.phi();
0640       float clus_pt = E_vec_cluster.perp();
0641       float clus_chisq = recoCluster->get_chi2();
0642
0643       h_clusE->Fill(clusE);
0644
0645       if (clusE < emcMinClusE1 || clusE > emcMaxClusE)
0646       {
0647         continue;
0648       }
0649       if (clus_chisq > 4)
0650       {
0651         continue;
0652       }
0653
0654       h_etaphi_clus->Fill(clus_eta, clus_phi);
0655
0656       TLorentzVector photon1;
0657       photon1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clusE);
0658
0659       for (clusterIter2 = clusterEnd.first; clusterIter2 != clusterEnd.second; clusterIter2++)
0660       {
0661         if (clusterIter == clusterIter2)
0662         {
0663           continue;
0664         }
0665         RawCluster* recoCluster2 = clusterIter2->second;
0666
0667         CLHEP::Hep3Vector vertex2(0, 0, 0);
0668         CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, vertex2);
0669
0670         float clus2E = E_vec_cluster2.mag();
0671         float clus2_eta = E_vec_cluster2.pseudoRapidity();
0672         float clus2_phi = E_vec_cluster2.phi();
0673         float clus2_pt = E_vec_cluster2.perp();
0674         float clus2_chisq = recoCluster2->get_chi2();
0675
0676         if (clus2E < emcMinClusE2 || clus2E > emcMaxClusE)
0677         {
0678           continue;
0679         }
0680         if (clus2_chisq > 4)
0681         {
0682           continue;
0683         }
0684
0685         TLorentzVector photon2;
0686         photon2.SetPtEtaPhiE(clus2_pt, clus2_eta, clus2_phi, clus2E);
0687
0688         if (sqrt(pow(clusE - clus2E, 2)) / (clusE + clus2E) > maxAlpha)
0689         {
0690           continue;
0691         }
0692
0693         TLorentzVector pi0 = photon1 + photon2;
0694         float pi0Mass = pi0.M();
0695         unsigned int lt_eta = recoCluster->get_lead_tower().first;
0696         unsigned int lt_phi = recoCluster->get_lead_tower().second;
0697
0698         int IB_num = ((lt_eta / 8) * 32) + (lt_phi / 8);
0699
0700         for (int bit : scaledActiveBits)
0701         {
0702           if (std::find(triggerIndices.begin(), triggerIndices.end(), bit) == triggerIndices.end())
0703           {
0704             continue;
0705           }
0706           h_pi0_trigIB_mass->Fill(
0707               static_cast<double>(bit),
0708               static_cast<double>(IB_num),
0709               static_cast<double>(pi0Mass));
0710         }
0711         h_InvMass->Fill(pi0Mass);
0712       }
0713     }  // end cluster loop
0714   }
0715
0716   //----------------- Trigger / alignment ----------------------------//
0717   float leading_cluster_ecore = 0;
0718   float leading_cluster_eta = 0;
0719   float leading_cluster_phi = 0;
0720   int evtNum_overK = event_number / 1000;
0721
0722   if (clusterContainer)
0723   {
0724     RawClusterContainer::ConstRange clusterEnd =
0725         clusterContainer->getClusters();
0726     RawClusterContainer::ConstIterator clusterIter;
0727     RawClusterContainer::ConstIterator clusterIter2;
0728
0729     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second;
0730          clusterIter++)
0731     {
0732       RawCluster* recoCluster = clusterIter->second;
0733       if (recoCluster->get_chi2() > 2)
0734       {
0735         continue;
0736       }
0737
0738       CLHEP::Hep3Vector vertex(0, 0, 0);
0739       CLHEP::Hep3Vector E_vec_cluster =
0740           RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0741
0742       float clusE = E_vec_cluster.mag();
0743       float clusEta = E_vec_cluster.pseudoRapidity();
0744       float clusPhi = E_vec_cluster.phi();
0745       if (clusE > leading_cluster_ecore)
0746       {
0747         leading_cluster_ecore = clusE;
0748         leading_cluster_eta = clusEta;
0749         leading_cluster_phi = clusPhi;
0750       }
0751     }
0752     for (int i = 0; i < 64; i++)
0753     {
0754       if (scaledBits[i])
0755       {
0756         pr_ldClus_trig->Fill(i, leading_cluster_ecore);
0757         if (!(std::find(trigOfInterest.begin(), trigOfInterest.end(), i) != trigOfInterest.end()))
0758         {
0759           continue;
0760         }
0761         h_edist[i]->Fill(leading_cluster_eta, leading_cluster_phi);
0762         h_ldClus_trig[i]->Fill(leading_cluster_ecore);
0763         pr_evtNum_ldClus_trig[i]->Fill(evtNum_overK, leading_cluster_ecore);
0764         if (raw[i] > 0)
0765         {
0766           pr_rejection[i]->Fill(evtNum_overK,
0767                                 (float) raw[10] / (float) raw[i]);
0768           pr_livetime[i]->Fill(evtNum_overK,
0769                                (float) live[i] / (float) raw[i]);
0770         }
0771       }
0772     }
0773   }
0774
0775   return Fun4AllReturnCodes::EVENT_OK;
0776 }
0777
0778 int CaloValid::End(PHCompositeNode* topNode)
0779 {
0780   auto* hm = QAHistManagerDef::getHistoManager();
0781   assert(hm);
0782
0783   //------EmCal-----//
0784   {
0785     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0786     if (towers)
0787     {
0788       int size = towers->size();
0789
0790       auto* h_CaloValid_cemc_etaphi_pedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}cemc_etaphi_pedRMS", getHistoPrefix())));
0791       auto* h_CaloValid_cemc_etaphi_ZSpedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}cemc_etaphi_ZSpedRMS", getHistoPrefix())));
0792
0793       for (int channel = 0; channel < size; channel++)
0794       {
0795         unsigned int towerkey = towers->encode_key(channel);
0796         int ieta = towers->getTowerEtaBin(towerkey);
0797         int iphi = towers->getTowerPhiBin(towerkey);
0798         float ped_rms = h_cemc_channel_pedestal[channel]->GetRMS();
0799         h_CaloValid_cemc_etaphi_pedRMS->Fill(ieta, iphi, ped_rms);
0800         MirrorHistogram(h_cemc_channel_energy[channel]);
0801         double rmsZS = h_cemc_channel_energy[channel]->GetRMS();
0802         h_CaloValid_cemc_etaphi_ZSpedRMS->Fill(ieta, iphi, rmsZS);
0803       }
0804     }
0805   }
0806   //------IHCal------//
0807   {
0808     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0809     if (towers)
0810     {
0811       int size = towers->size();
0812
0813       auto* h_CaloValid_ihcal_etaphi_pedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}ihcal_etaphi_pedRMS", getHistoPrefix())));
0814       auto* h_CaloValid_ihcal_etaphi_ZSpedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}ihcal_etaphi_ZSpedRMS", getHistoPrefix())));
0815
0816       for (int channel = 0; channel < size; channel++)
0817       {
0818         unsigned int towerkey = towers->encode_key(channel);
0819         int ieta = towers->getTowerEtaBin(towerkey);
0820         int iphi = towers->getTowerPhiBin(towerkey);
0821         float ped_rms = h_ihcal_channel_pedestal[channel]->GetRMS();
0822         h_CaloValid_ihcal_etaphi_pedRMS->Fill(ieta, iphi, ped_rms);
0823         MirrorHistogram(h_ihcal_channel_energy[channel]);
0824         double rmsZS = h_ihcal_channel_energy[channel]->GetRMS();
0825         h_CaloValid_ihcal_etaphi_ZSpedRMS->Fill(ieta, iphi, rmsZS);
0826       }
0827     }
0828   }
0829   //------OHCal-----//
0830   {
0831     TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0832     if (towers)
0833     {
0834       int size = towers->size();
0835
0836       auto* h_CaloValid_ohcal_etaphi_pedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}ohcal_etaphi_pedRMS", getHistoPrefix())));
0837       auto* h_CaloValid_ohcal_etaphi_ZSpedRMS = dynamic_cast<TProfile2D*>(hm->getHisto(std::format("{}ohcal_etaphi_ZSpedRMS", getHistoPrefix())));
0838
0839       for (int channel = 0; channel < size; channel++)
0840       {
0841         unsigned int towerkey = towers->encode_key(channel);
0842         int ieta = towers->getTowerEtaBin(towerkey);
0843         int iphi = towers->getTowerPhiBin(towerkey);
0844         float ped_rms = h_ohcal_channel_pedestal[channel]->GetRMS();
0845         h_CaloValid_ohcal_etaphi_pedRMS->Fill(ieta, iphi, ped_rms);
0846         MirrorHistogram(h_ohcal_channel_energy[channel]);
0847         double rmsZS = h_ohcal_channel_energy[channel]->GetRMS();
0848         h_CaloValid_ohcal_etaphi_ZSpedRMS->Fill(ieta, iphi, rmsZS);
0849       }
0850     }
0851   }
0852   return Fun4AllReturnCodes::EVENT_OK;
0853 }
0854
0855 void CaloValid::MirrorHistogram(TH1* h)
0856 {
0857   int middleBin = h->GetXaxis()->FindBin(0.0);
0858
0859   for (int i = 1; i < middleBin; ++i)
0860   {
0861     int correspondingBin = middleBin + (middleBin - i);
0862     float negValue = h->GetBinContent(i);
0863     h->SetBinContent(correspondingBin, negValue);
0864   }
0865 }
0866
0867 TH2* CaloValid::LogYHist2D(const std::string& name, const std::string& title, int xbins_in, double xmin, double xmax, int ybins_in, double ymin, double ymax)
0868 {
0869   Double_t logymin = std::log10(ymin);
0870   Double_t logymax = std::log10(ymax);
0871   Double_t binwidth = (logymax - logymin) / ybins_in;
0872   Double_t* ybins = new Double_t[ybins_in + 2];  // allocate 1 extra bin to fix "malloc()L memory corruption crash
0873
0874   for (Int_t i = 0; i <= ybins_in + 1; i++)
0875   {
0876     ybins[i] = pow(10, logymin + (i * binwidth));
0877   }
0878
0879   TH2F* h = new TH2F(name.c_str(), title.c_str(), xbins_in, xmin, xmax, ybins_in, ybins);
0880   delete[] ybins;
0881   return h;
0882 }
0883 std::string CaloValid::getHistoPrefix() const { return std::string("h_") + Name() + std::string("_"); }
0884
0885 void CaloValid::createHistos()
0886 {
0887   auto* hm = QAHistManagerDef::getHistoManager();
0888   assert(hm);
0889
0890   // create and register your histos (all types) here
0891   h_emcal_mbd_correlation = new TH2F(std::format("{}emcal_mbd_correlation", getHistoPrefix()).c_str(), ";emcal;mbd", 100, 0, 1, 100, 0, 1);
0892   h_emcal_mbd_correlation->SetDirectory(nullptr);
0893   hm->registerHisto(h_emcal_mbd_correlation);
0894
0895   h_mbd_hits = new TH1F(std::format("{}mbd_hits", getHistoPrefix()).c_str(), "mb hits", 100, 0, 100);
0896   h_mbd_hits->SetDirectory(nullptr);
0897   hm->registerHisto(h_mbd_hits);
0898
0899   h_ohcal_mbd_correlation = new TH2F(std::format("{}ohcal_mbd_correlation", getHistoPrefix()).c_str(), ";ohcal;mbd", 100, 0, 1, 100, 0, 1);
0900   h_ohcal_mbd_correlation->SetDirectory(nullptr);
0901   hm->registerHisto(h_ohcal_mbd_correlation);
0902
0903   h_ihcal_mbd_correlation = new TH2F(std::format("{}ihcal_mbd_correlation", getHistoPrefix()).c_str(), ";ihcal;mbd", 100, 0, 1, 100, 0, 1);
0904   h_ihcal_mbd_correlation->SetDirectory(nullptr);
0905   hm->registerHisto(h_ihcal_mbd_correlation);
0906
0907   h_emcal_hcal_correlation = new TH2F(std::format("{}emcal_hcal_correlation", getHistoPrefix()).c_str(), ";emcal;hcal", 100, 0, 1, 100, 0, 1);
0908   h_emcal_hcal_correlation->SetDirectory(nullptr);
0909   hm->registerHisto(h_emcal_hcal_correlation);
0910
0911   h_cemc_etaphi = new TH2F(std::format("{}cemc_etaphi", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256);
0912   h_cemc_etaphi->SetDirectory(nullptr);
0913   hm->registerHisto(h_cemc_etaphi);
0914
0915   h_ihcal_etaphi = new TH2F(std::format("{}ihcal_etaphi", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64);
0916   h_ihcal_etaphi->SetDirectory(nullptr);
0917   hm->registerHisto(h_ihcal_etaphi);
0918
0919   h_ohcal_etaphi = new TH2F(std::format("{}ohcal_etaphi", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64);
0920   h_ohcal_etaphi->SetDirectory(nullptr);
0921   hm->registerHisto(h_ohcal_etaphi);
0922
0923   h_cemc_etaphi_wQA = new TH2F(std::format("{}cemc_etaphi_wQA", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256);
0924   h_cemc_etaphi_wQA->SetDirectory(nullptr);
0925   hm->registerHisto(h_cemc_etaphi_wQA);
0926
0927   h_ihcal_etaphi_wQA = new TH2F(std::format("{}ihcal_etaphi_wQA", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64);
0928   h_ihcal_etaphi_wQA->SetDirectory(nullptr);
0929   hm->registerHisto(h_ihcal_etaphi_wQA);
0930
0931   h_ohcal_etaphi_wQA = new TH2F(std::format("{}ohcal_etaphi_wQA", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64);
0932   h_ohcal_etaphi_wQA->SetDirectory(nullptr);
0933   hm->registerHisto(h_ohcal_etaphi_wQA);
0934
0935   h_ihcal_status = new TH1F(std::format("{}ihcal_status", getHistoPrefix()).c_str(), "", 256, 0, 256);
0936   h_ihcal_status->SetDirectory(nullptr);
0937   hm->registerHisto(h_ihcal_status);
0938
0939   h_ohcal_status = new TH1F(std::format("{}ohcal_status", getHistoPrefix()).c_str(), "", 256, 0, 256);
0940   h_ohcal_status->SetDirectory(nullptr);
0941   hm->registerHisto(h_ohcal_status);
0942
0943   h_cemc_status = new TH1F(std::format("{}cemc_status", getHistoPrefix()).c_str(), "", 256, 0, 256);
0944   h_cemc_status->SetDirectory(nullptr);
0945   hm->registerHisto(h_cemc_status);
0946
0947   h_cemc_e_chi2 = LogYHist2D(std::format("{}cemc_e_chi2", getHistoPrefix()), "", 270, -2, 25, 1000, 0.5, 4e8);
0948   h_cemc_e_chi2->SetDirectory(nullptr);
0949   hm->registerHisto(h_cemc_e_chi2);
0950
0951   h_ihcal_e_chi2 = LogYHist2D(std::format("{}ihcal_e_chi2", getHistoPrefix()), "", 270, -2, 25, 1000, 0.5, 4e8);
0952   h_ihcal_e_chi2->SetDirectory(nullptr);
0953   hm->registerHisto(h_ihcal_e_chi2);
0954
0955   h_ohcal_e_chi2 = LogYHist2D(std::format("{}ohcal_e_chi2", getHistoPrefix()), "", 270, -2, 25, 1000, 0.5, 4e8);
0956   h_ohcal_e_chi2->SetDirectory(nullptr);
0957   hm->registerHisto(h_ohcal_e_chi2);
0958
0959   h_cemc_etaphi_time = new TProfile2D(std::format("{}cemc_etaphi_time", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
0960   h_cemc_etaphi_time->SetDirectory(nullptr);
0961   hm->registerHisto(h_cemc_etaphi_time);
0962
0963   h_cemc_etaphi_time_raw = new TProfile2D(std::format("{}cemc_etaphi_time_raw", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
0964   h_cemc_etaphi_time_raw->SetDirectory(nullptr);
0965   hm->registerHisto(h_cemc_etaphi_time_raw);
0966
0967   h_ihcal_etaphi_time = new TProfile2D(std::format("{}ihcal_etaphi_time", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0968   h_ihcal_etaphi_time->SetDirectory(nullptr);
0969   hm->registerHisto(h_ihcal_etaphi_time);
0970
0971   h_ihcal_etaphi_time_raw = new TProfile2D(std::format("{}ihcal_etaphi_time_raw", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0972   h_ihcal_etaphi_time_raw->SetDirectory(nullptr);
0973   hm->registerHisto(h_ihcal_etaphi_time_raw);
0974
0975   h_ohcal_etaphi_time = new TProfile2D(std::format("{}ohcal_etaphi_time", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0976   h_ohcal_etaphi_time->SetDirectory(nullptr);
0977   hm->registerHisto(h_ohcal_etaphi_time);
0978
0979   h_ohcal_etaphi_time_raw = new TProfile2D(std::format("{}ohcal_etaphi_time_raw", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0980   h_ohcal_etaphi_time_raw->SetDirectory(nullptr);
0981   hm->registerHisto(h_ohcal_etaphi_time_raw);
0982
0983   h_cemc_etaphi_fracHitADC = new TProfile2D(std::format("{}cemc_etaphi_fracHitADC", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
0984   h_cemc_etaphi_fracHitADC->SetDirectory(nullptr);
0985   hm->registerHisto(h_cemc_etaphi_fracHitADC);
0986
0987   h_cemc_etaphi_fracHit = new TProfile2D(std::format("{}cemc_etaphi_fracHit", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
0988   h_cemc_etaphi_fracHit->SetDirectory(nullptr);
0989   hm->registerHisto(h_cemc_etaphi_fracHit);
0990
0991   h_ihcal_etaphi_fracHitADC = new TProfile2D(std::format("{}ihcal_etaphi_fracHitADC", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0992   h_ihcal_etaphi_fracHitADC->SetDirectory(nullptr);
0993   hm->registerHisto(h_ihcal_etaphi_fracHitADC);
0994
0995   h_ohcal_etaphi_fracHitADC = new TProfile2D(std::format("{}ohcal_etaphi_fracHitADC", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
0996   h_ohcal_etaphi_fracHitADC->SetDirectory(nullptr);
0997   hm->registerHisto(h_ohcal_etaphi_fracHitADC);
0998
0999   h_cemc_etaphi_pedRMS = new TProfile2D(std::format("{}cemc_etaphi_pedRMS", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, 0, 1000);
1000   h_cemc_etaphi_pedRMS->SetDirectory(nullptr);
1001   hm->registerHisto(h_cemc_etaphi_pedRMS);
1002
1003   h_ohcal_etaphi_pedRMS = new TProfile2D(std::format("{}ohcal_etaphi_pedRMS", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, 0, 1000);
1004   h_ohcal_etaphi_pedRMS->SetDirectory(nullptr);
1005   hm->registerHisto(h_ohcal_etaphi_pedRMS);
1006
1007   h_ihcal_etaphi_pedRMS = new TProfile2D(std::format("{}ihcal_etaphi_pedRMS", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, 0, 1000);
1008   h_ihcal_etaphi_pedRMS->SetDirectory(nullptr);
1009   hm->registerHisto(h_ihcal_etaphi_pedRMS);
1010
1011   h_cemc_etaphi_ZSpedRMS = new TProfile2D(std::format("{}cemc_etaphi_ZSpedRMS", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, 0, 1000);
1012   h_cemc_etaphi_ZSpedRMS->SetDirectory(nullptr);
1013   hm->registerHisto(h_cemc_etaphi_ZSpedRMS);
1014
1015   h_ohcal_etaphi_ZSpedRMS = new TProfile2D(std::format("{}ohcal_etaphi_ZSpedRMS", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, 0, 1000);
1016   h_ohcal_etaphi_ZSpedRMS->SetDirectory(nullptr);
1017   hm->registerHisto(h_ohcal_etaphi_ZSpedRMS);
1018
1019   h_ihcal_etaphi_ZSpedRMS = new TProfile2D(std::format("{}ihcal_etaphi_ZSpedRMS", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, 0, 1000);
1020   h_ihcal_etaphi_ZSpedRMS->SetDirectory(nullptr);
1021   hm->registerHisto(h_ihcal_etaphi_ZSpedRMS);
1022
1023   h_cemc_etaphi_badChi2 = new TProfile2D(std::format("{}cemc_etaphi_badChi2", getHistoPrefix()).c_str(), ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
1024   h_cemc_etaphi_badChi2->SetDirectory(nullptr);
1025   hm->registerHisto(h_cemc_etaphi_badChi2);
1026
1027   h_ihcal_etaphi_badChi2 = new TProfile2D(std::format("{}ihcal_etaphi_badChi2", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
1028   h_ihcal_etaphi_badChi2->SetDirectory(nullptr);
1029   hm->registerHisto(h_ihcal_etaphi_badChi2);
1030
1031   h_ohcal_etaphi_badChi2 = new TProfile2D(std::format("{}ohcal_etaphi_badChi2", getHistoPrefix()).c_str(), ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
1032   h_ohcal_etaphi_badChi2->SetDirectory(nullptr);
1033   hm->registerHisto(h_ohcal_etaphi_badChi2);
1034
1035   // 1D distributions
1036   h_InvMass = new TH1F(std::format("{}InvMass", getHistoPrefix()).c_str(), "Invariant Mass", 120, 0, 1.2);
1037   h_InvMass->SetDirectory(nullptr);
1038   hm->registerHisto(h_InvMass);
1039
1040   // for (int channel = 0; channel < 128*192; channel++) {
1041   h_channel_pedestal_0 = new TH1F(std::format("{}channel_pedestal_0", getHistoPrefix()).c_str(), "Test Pedestal", 1000, -500., 500.);
1042   h_channel_pedestal_0->SetDirectory(nullptr);
1043   hm->registerHisto(h_channel_pedestal_0);
1044
1045   // vertex distributions
1046   h_vtx_z_raw = new TH1D(std::format("{}vtx_z_raw", getHistoPrefix()).c_str(), "hvtx_z_raw", 201, -100.5, 100.5);
1047   h_vtx_z_raw->SetDirectory(nullptr);
1048   hm->registerHisto(h_vtx_z_raw);
1049
1050   h_vtx_z_cut = new TH1D(std::format("{}vtx_z_cut", getHistoPrefix()).c_str(), "hvtx_z_cut", 201, -100.5, 100.5);
1051   h_vtx_z_cut->SetDirectory(nullptr);
1052   hm->registerHisto(h_vtx_z_cut);
1053
1054   // raw timing information
1055   h_emcaltime_cut = new TH1D(std::format("{}emcaltime_cut", getHistoPrefix()).c_str(), "hemcaltime_cut", 50, -17.5, 32.5);
1056   h_emcaltime_cut->SetDirectory(nullptr);
1057   hm->registerHisto(h_emcaltime_cut);
1058
1059   h_ihcaltime_cut = new TH1D(std::format("{}ihcaltime_cut", getHistoPrefix()).c_str(), "hihcaltime_cut", 50, -17.5, 32.5);
1060   h_ihcaltime_cut->SetDirectory(nullptr);
1061   hm->registerHisto(h_ihcaltime_cut);
1062
1063   h_ohcaltime_cut = new TH1D(std::format("{}ohcaltime_cut", getHistoPrefix()).c_str(), "hohcaltime_cut", 50, -17.5, 32.5);
1064   h_ohcaltime_cut->SetDirectory(nullptr);
1065   hm->registerHisto(h_ohcaltime_cut);
1066
1067   // extracted timing information
1068   h_emcaltime = new TH1D(std::format("{}emcaltime", getHistoPrefix()).c_str(), "hemcaltime", 50, -17.5, 32.5);
1069   h_emcaltime->SetDirectory(nullptr);
1070   hm->registerHisto(h_emcaltime);
1071
1072   h_ihcaltime = new TH1D(std::format("{}ihcaltime", getHistoPrefix()).c_str(), "hihcaltime", 50, -17.5, 32.5);
1073   h_ihcaltime->SetDirectory(nullptr);
1074   hm->registerHisto(h_ihcaltime);
1075
1076   h_ohcaltime = new TH1D(std::format("{}ohcaltime", getHistoPrefix()).c_str(), "hohcaltime", 50, -17.5, 32.5);
1077   h_ohcaltime->SetDirectory(nullptr);
1078   hm->registerHisto(h_ohcaltime);
1079
1080   h_emcal_tower_e = new TH1F(std::format("{}emcal_tower_e", getHistoPrefix()).c_str(), "emcal_tower_e", 5000, -0.1, 1);
1081   h_emcal_tower_e->SetDirectory(nullptr);
1082   hm->registerHisto(h_emcal_tower_e);
1083
1084   h_emcal_tower_e_wide_range = new TH1F(std::format("{}emcal_tower_e_wide_range", getHistoPrefix()).c_str(), "emcal_tower_e_wide_range", 1000, -10., 40.);
1085   h_emcal_tower_e_wide_range->SetDirectory(nullptr);
1086   hm->registerHisto(h_emcal_tower_e_wide_range);
1087
1088   h_emcal_tower_e_saturated = new TH1F(std::format("{}emcal_tower_e_saturated", getHistoPrefix()).c_str(), "emcal_tower_e_saturated", 1000, -10., 40.);
1089   h_emcal_tower_e_saturated->SetDirectory(nullptr);
1090   hm->registerHisto(h_emcal_tower_e_saturated);
1091
1092   h_ihcal_tower_e = new TH1F(std::format("{}ihcal_tower_e", getHistoPrefix()).c_str(), "ihcal_tower_e", 5000, -0.1, 1);
1093   h_ihcal_tower_e->SetDirectory(nullptr);
1094   hm->registerHisto(h_ihcal_tower_e);
1095
1096   h_ihcal_tower_e_wide_range = new TH1F(std::format("{}ihcal_tower_e_wide_range", getHistoPrefix()).c_str(), "ihcal_tower_e_wide_range", 1000, -10., 40.);
1097   h_ihcal_tower_e_wide_range->SetDirectory(nullptr);
1098   hm->registerHisto(h_ihcal_tower_e_wide_range);
1099
1100   h_ihcal_tower_e_saturated = new TH1F(std::format("{}ihcal_tower_e_saturated", getHistoPrefix()).c_str(), "ihcal_tower_e_saturated", 1000, -10., 40.);
1101   h_ihcal_tower_e_saturated->SetDirectory(nullptr);
1102   hm->registerHisto(h_ihcal_tower_e_saturated);
1103
1104   h_ohcal_tower_e = new TH1F(std::format("{}ohcal_tower_e", getHistoPrefix()).c_str(), "ohcal_tower_e", 5000, -0.1, 1);
1105   h_ohcal_tower_e->SetDirectory(nullptr);
1106   hm->registerHisto(h_ohcal_tower_e);
1107
1108   h_ohcal_tower_e_wide_range = new TH1F(std::format("{}ohcal_tower_e_wide_range", getHistoPrefix()).c_str(), "ohcal_tower_e_wide_range", 1000, -10., 40.);
1109   h_ohcal_tower_e_wide_range->SetDirectory(nullptr);
1110   hm->registerHisto(h_ohcal_tower_e_wide_range);
1111
1112   h_ohcal_tower_e_saturated = new TH1F(std::format("{}ohcal_tower_e_saturated", getHistoPrefix()).c_str(), "ohcal_tower_e_saturated", 1000, -10., 40.);
1113   h_ohcal_tower_e_saturated->SetDirectory(nullptr);
1114   hm->registerHisto(h_ohcal_tower_e_saturated);
1115
1116   // cluster QA
1117   h_etaphi_clus = new TH2F(std::format("{}etaphi_clus", getHistoPrefix()).c_str(), "", 140, -1.2, 1.2, 64, -1 * M_PI, M_PI);
1118   h_etaphi_clus->SetDirectory(nullptr);
1119   hm->registerHisto(h_etaphi_clus);
1120
1121   h_clusE = new TH1F(std::format("{}clusE", getHistoPrefix()).c_str(), "", 100, 0, 10);
1122   h_clusE->SetDirectory(nullptr);
1123   hm->registerHisto(h_clusE);
1124
1125   h_triggerVec = new TH1F(std::format("{}triggerVec", getHistoPrefix()).c_str(), "", 64, 0, 64);
1126   h_triggerVec->SetDirectory(nullptr);
1127   hm->registerHisto(h_triggerVec);
1128
1129   //---------EMCal--------//
1130   {
1131     int size = 128 * 192;
1132     for (int channel = 0; channel < size; channel++)
1133     {
1134       std::string hname = std::format("h_cemc_channel_pedestal_{}", channel);
1135       h_cemc_channel_pedestal[channel] = new TH1F(hname.c_str(), hname.c_str(), 2000, -0.5, 2000.5);
1136       h_cemc_channel_pedestal[channel]->SetDirectory(nullptr);
1137
1138       std::string hnameE = std::format("h_cemc_channel_energy_{}", channel);
1139       h_cemc_channel_energy[channel] = new TH1F(hnameE.c_str(), hnameE.c_str(), 1000, -50, 50);
1140       h_cemc_channel_energy[channel]->SetDirectory(nullptr);
1141     }
1142   }
1143   //--------OHCal--------//
1144   {
1145     int size = 32 * 48;
1146     for (int channel = 0; channel < size; channel++)
1147     {
1148       std::string hname = std::format("h_ohcal_channel_pedestal_{}", channel);
1149       h_ohcal_channel_pedestal[channel] = new TH1F(hname.c_str(), hname.c_str(), 2000, -0.5, 2000.5);
1150       h_ohcal_channel_pedestal[channel]->SetDirectory(nullptr);
1151
1152       std::string hnameE = std::format("h_ohcal_channel_energy_{}", channel);
1153       h_ohcal_channel_energy[channel] = new TH1F(hnameE.c_str(), hnameE.c_str(), 1000, -50, 50);
1154       h_ohcal_channel_energy[channel]->SetDirectory(nullptr);
1155     }
1156   }
1157   //--------IHCal-------//
1158   {
1159     int size = 32 * 48;
1160     for (int channel = 0; channel < size; channel++)
1161     {
1162       std::string hname = std::format("h_ihcal_channel_pedestal_{}", channel);
1163       h_ihcal_channel_pedestal[channel] = new TH1F(hname.c_str(), hname.c_str(), 2000, -0.5, 2000.5);
1164       h_ihcal_channel_pedestal[channel]->SetDirectory(nullptr);
1165
1166       std::string hnameE = std::format("h_ihcal_channel_energy_{}", channel);
1167       h_ihcal_channel_energy[channel] = new TH1F(hnameE.c_str(), hnameE.c_str(), 1000, -50, 50);
1168       h_ihcal_channel_energy[channel]->SetDirectory(nullptr);
1169     }
1170   }
1171   h_pi0_trigIB_mass = new TH3F(
1172       "h_pi0_trigIB_mass",
1173       ";Trigger Bit; iIB (eta*32 + phi); #pi^{0} Mass (GeV/c^{2})",
1174       64, -0.5, 63.5,    // triggers
1175       384, -0.5, 383.5,  // IB index
1176       120, 0.0, 1.2      // mass range
1177   );
1178   h_pi0_trigIB_mass->SetDirectory(nullptr);
1179   hm->registerHisto(h_pi0_trigIB_mass);
1180
1181   // Trigger QA
1182   // h_triggerVec = new TH1F("h_CaloValid_triggerVec", "", 64, 0, 64); commented out due to memory allocation issue w/ redef from line 1009
1183   pr_ldClus_trig =
1184       new TProfile("pr_CaloValid_ldClus_trig", "", 64, 0, 64, 0, 10);
1185   for (int i = 0; i < 64; i++)
1186   {
1187     if (!(std::find(trigOfInterest.begin(), trigOfInterest.end(), i) != trigOfInterest.end()))
1188     {
1189       continue;
1190     }
1191     h_edist[i] = new TH2F(
1192         std::format("h_CaloValid_edist_trig{}", i).c_str(), "",
1193         64, -1.2, 1.2, 128, -3.1415, 3.1415);
1194     h_ldClus_trig[i] = new TH1F(
1195         std::format("h_CaloValid_ldClus_trig{}", i).c_str(), "",
1196         18, 1, 10);
1197     pr_evtNum_ldClus_trig[i] = new TProfile(
1198         std::format("pr_CaloValid_evtNum_ldClus_trig{}", i)
1199             .c_str(),
1200         "", 100000, 0, 100000, 0, 10);
1201     pr_rejection[i] = new TProfile(
1202         std::format("pr_CaloValid_rejection_trig{}", i).c_str(),
1203         "", 100000, 0, 100000, 0, 50000);
1204     pr_livetime[i] = new TProfile(
1205         std::format("pr_CaloValid_livetime_trig{}", i).c_str(),
1206         "", 100000, 0, 100000, 0, 10);
1207
1208     hm->registerHisto(h_edist[i]);
1209     hm->registerHisto(h_ldClus_trig[i]);
1210     hm->registerHisto(pr_evtNum_ldClus_trig[i]);
1211     hm->registerHisto(pr_rejection[i]);
1212     hm->registerHisto(pr_livetime[i]);
1213   }
1214   hm->registerHisto(h_triggerVec);
1215   hm->registerHisto(pr_ldClus_trig);
1216 }