Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-17 09:21:16

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