File indexing completed on 2025-12-16 09:18:07
0001 #include "caloHistGen.h"
0002
0003
0004 #include <calobase/RawCluster.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawClusterUtility.h>
0007 #include <calobase/RawTowerDefs.h>
0008
0009
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012 #include <calobase/TowerInfoDefs.h>
0013
0014
0015 #include <globalvertex/GlobalVertex.h>
0016 #include <globalvertex/GlobalVertexMap.h>
0017
0018
0019 #include <ffarawobjects/Gl1Packet.h>
0020
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/Fun4AllServer.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>
0027
0028
0029 #include <TFile.h>
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TH3.h>
0033 #include <TTree.h>
0034
0035
0036 #include <CLHEP/Vector/ThreeVector.h>
0037 #include <TLorentzVector.h>
0038 #include <cstdint>
0039 #include <iostream>
0040 #include <map>
0041 #include <utility>
0042
0043
0044 caloHistGen::caloHistGen(const std::string &name)
0045 : SubsysReco("CaloHistGen")
0046 , Outfile(name)
0047 {
0048 std::cout << "caloHistGen::caloHistGen(const std::string &name) Calling ctor" << std::endl;
0049 }
0050
0051
0052 int caloHistGen::Init(PHCompositeNode * )
0053 {
0054 delete out;
0055 out = new TFile(Outfile.c_str(), "RECREATE");
0056
0057 h_emcTowE = new TH3F("emcTowE", "tower eta, tower phi, tower energy", 96, -0.5, 95.5, 256, -0.5, 256.5, 1000, -10, 25);
0058 h_OHCalTowE = new TH3F("OHCalTowE", "outer hcal tower eta, tower phi, tower energy", 24, -0.5, 23.5, 64, -0.5, 63.5, 1000, -10, 25);
0059 h_IHCalTowE = new TH3F("IHCalTowE", "inner hcal tower eta, tower phi, tower energy", 24, -0.5, 23.5, 64, -0.5, 63.5, 1000, -10, 25);
0060
0061 h_emcTowChi2 = new TH3F("emcTowChi2", "tower eta, tower phi, tower Chi2", 96, -0.5, 95.5, 256, -0.5, 256.5, 1000, 0.5, 4e8);
0062 h_OHCalTowChi2 = new TH3F("OHCalTowChi2", "tower eta, tower phi, tower Chi2", 24, -0.5, 23.5, 64, -0.5, 63.5, 1000, 0.5, 4e8);
0063 h_IHCalTowChi2 = new TH3F("IHCalTowChi2", "tower eta, tower phi, tower Chi2", 24, -0.5, 23.5, 64, -0.5, 63.5, 1000, 0.5, 4e8);
0064
0065 h_emcTowTime = new TH3F("emcTowEnergy", "tower eta, tower phi, tower Time", 96, -0.5, 95.5, 256, -0.5, 256.5, 21, -10.5, 10.5);
0066 h_OHCalTowTime = new TH3F("OHCalTowEnergy", "tower eta, tower phi, tower Time", 24, -0.5, 23.5, 64, -0.5, 63.5, 21, -10.5, 10.5);
0067 h_IHCalTowTime = new TH3F("IHCalTowTime", "tower eta, tower phi, tower Time", 24, -0.5, 23.5, 64, -0.5, 63.5, 21, -10.5, 10.5);
0068
0069 h_emcTowPed = new TH3F("emcTowPed", "tower eta, tower phi, tower Ped", 96, -0.5, 95.5, 256, -0.5, 256.5, 100, 0, 5000);
0070 h_OHCalTowPed = new TH3F("OHCalTowPed", "tower eta, tower phi, tower Ped", 24, -0.5, 23.5, 64, -0.5, 63.5, 100, 0, 5000);
0071 h_IHCalTowPed = new TH3F("IHCalTowPed", "tower eta, tower phi, tower Ped", 24, -0.5, 23.5, 64, -0.5, 63.5, 100, 0, 5000);
0072
0073 h_clusInfo = new TH3F("clusInfo", "cluster eta, phi, energy", 140, -1.2, 1.2, 100, -1. * M_PI, M_PI, 500, -2, 25);
0074
0075 h_clusPt = new TH1F("clusPt", "cluster pT", 500, -2, 25);
0076 h_clusEcore = new TH1F("clusEcore", "cluster Ecore", 500, -2, 25);
0077
0078 h_clusChi2 = new TH1F("clusChi2_E", "cluster chi2", 100, 0, 100);
0079
0080 h_zVertex = new TH1F("zVertex", "z vertex from mbd", 200, -100, 100);
0081
0082 h_zdcNSE = new TH2F("zscNSChanE", "zdc N/S Energy", 2, -0.5, 1.5, 100, 0, 500);
0083
0084 h_zdcChanTime = new TH2F("zdcChanTime", "zdc timing per channel", 7, -0.5, 6.5, 21, -10.5, 10.5);
0085
0086 h_diPhotonEtaPhiE = new TH3F("h_diPhontonEtaPhiE", "diphoton spatial kinematics", 100, -1.1, 1.1, 100, -1 * M_PI, M_PI, 40, 0, 20);
0087
0088 h_diPhotonMassE = new TH2F("h_diPhontonMassE", "diphoton mass and energy", 200, 0, 1, 40, 0, 20);
0089
0090
0091 Fun4AllServer *se = Fun4AllServer::instance();
0092 se->Print("NODETREE");
0093
0094 std::cout << "caloHistGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0095 return Fun4AllReturnCodes::EVENT_OK;
0096 }
0097
0098
0099 int caloHistGen::process_event(PHCompositeNode *topNode)
0100 {
0101 float vertex = -9999;
0102 bool hasMBDvertex = true;
0103
0104 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0105 if (!vertexmap)
0106 {
0107 std::cout << "GlobalVertexMap node is missing" << std::endl;
0108 }
0109 if (vertexmap->empty())
0110 {
0111 std::cout << "caloHistGen::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0112
0113 hasMBDvertex = false;
0114 }
0115
0116 if (vertexmap && !vertexmap->empty())
0117 {
0118 GlobalVertex *vtx = vertexmap->begin()->second;
0119 if (vtx)
0120 {
0121 vertex = vtx->get_z();
0122 h_zVertex->Fill(vertex);
0123 }
0124 }
0125
0126
0127
0128 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, m_trigNode.c_str());
0129 if (!gl1PacketInfo)
0130 {
0131 std::cout << PHWHERE << "caloHistGen::process_event: " << m_trigNode << " node is missing. Output related to this node will be empty" << std::endl;
0132 }
0133
0134 if (gl1PacketInfo && checkTrig)
0135 {
0136 uint64_t triggervec = gl1PacketInfo->getScaledVector();
0137 for (int i = 0; i < 64; i++)
0138 {
0139 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0140
0141 triggervec = (triggervec >> 1U) & 0xffffffffU;
0142 if (!trig_decision && i == trigRequired[i])
0143 {
0144 return 0;
0145 }
0146 }
0147 }
0148
0149
0150 RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_clusterNode.c_str());
0151 if (!clusterContainer && storeClusters)
0152 {
0153 std::cout << PHWHERE << "caloHistGen::process_event: " << m_clusterNode << " node is missing. Output related to this node will be empty" << std::endl;
0154 return 0;
0155 }
0156
0157
0158 TowerInfoContainer *emcTowerContainer;
0159 emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0160 if (!emcTowerContainer && storeEMCal)
0161 {
0162 std::cout << PHWHERE << "caloHistGen::process_event: " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0163 }
0164
0165
0166 unsigned int tower_range = 0;
0167 if (storeEMCal && emcTowerContainer)
0168 {
0169 tower_range = emcTowerContainer->size();
0170 for (unsigned int iter = 0; iter < tower_range; iter++)
0171 {
0172 unsigned int towerkey = emcTowerContainer->encode_key(iter);
0173 unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0174 unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0175 double energy = emcTowerContainer->get_tower_at_channel(iter)->get_energy();
0176 int time = emcTowerContainer->get_tower_at_channel(iter)->get_time();
0177 float chi2 = emcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0178 float pedestal = emcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0179
0180 if (emcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0181 {
0182 h_emcTowE->Fill(ieta, iphi, energy);
0183 h_emcTowChi2->Fill(ieta, iphi, chi2);
0184 h_emcTowTime->Fill(ieta, iphi, time);
0185 h_emcTowPed->Fill(ieta, iphi, pedestal);
0186 }
0187 }
0188 }
0189
0190 if (storeClusters && storeEMCal)
0191 {
0192 RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0193 RawClusterContainer::ConstIterator clusterIter;
0194 for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0195 {
0196 RawCluster *recoCluster = clusterIter->second;
0197
0198 CLHEP::Hep3Vector hep_vertex(0, 0, 0);
0199 if (vertex != -9999)
0200 {
0201 hep_vertex.setZ(vertex);
0202 }
0203 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, hep_vertex);
0204 CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, hep_vertex);
0205
0206 float clusE = E_vec_cluster_Full.mag();
0207 float clusEcore = E_vec_cluster.mag();
0208 float clus_eta = E_vec_cluster.pseudoRapidity();
0209 float clus_phi = E_vec_cluster.phi();
0210 float clus_pt = E_vec_cluster.perp();
0211
0212 h_clusInfo->Fill(clus_eta, clus_phi, clusE);
0213 h_clusPt->Fill(clus_pt);
0214 h_clusEcore->Fill(clusEcore);
0215 }
0216 }
0217
0218
0219 if (doPi0Reco && storeEMCal)
0220 {
0221 RawClusterContainer::ConstRange clusters = clusterContainer->getClusters();
0222 RawClusterContainer::ConstIterator clusterIter1, clusterIter2;
0223 int clusterCounter1 = 0;
0224 for (clusterIter1 = clusters.first; clusterIter1 != clusters.second; clusterIter1++)
0225 {
0226 RawCluster *recoCluster1 = clusterIter1->second;
0227 clusterCounter1++;
0228 if (recoCluster1->get_chi2() > 4)
0229 {
0230 continue;
0231 }
0232
0233 CLHEP::Hep3Vector hep_vertex(0, 0, 0);
0234 if (vertex != -9999)
0235 {
0236 hep_vertex.setZ(vertex);
0237 }
0238 CLHEP::Hep3Vector E_vec_cluster1 = RawClusterUtility::GetECoreVec(*recoCluster1, hep_vertex);
0239
0240 int clusterCounter2 = 0;
0241 for (clusterIter2 = clusters.first; clusterIter2 != clusters.second; clusterIter2++)
0242 {
0243 if (clusterCounter2 <= clusterCounter1)
0244 {
0245 continue;
0246 }
0247 clusterCounter2++;
0248 RawCluster *recoCluster2 = clusterIter2->second;
0249 CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, hep_vertex);
0250 if (recoCluster2->get_chi2() > 4)
0251 {
0252 continue;
0253 }
0254 if (E_vec_cluster1.mag() < clusEMin || E_vec_cluster2.mag() < clusEMin)
0255 {
0256 continue;
0257 }
0258 if (std::abs(E_vec_cluster1.mag() - E_vec_cluster2.mag()) / (E_vec_cluster1.mag() + E_vec_cluster2.mag()) > maxAlpha)
0259 {
0260 continue;
0261 }
0262 TLorentzVector photon1, photon2, pi0;
0263 photon1.SetPtEtaPhiE(E_vec_cluster1.perp(), E_vec_cluster1.pseudoRapidity(), E_vec_cluster1.phi(), E_vec_cluster1.mag());
0264 photon2.SetPtEtaPhiE(E_vec_cluster2.perp(), E_vec_cluster2.pseudoRapidity(), E_vec_cluster2.phi(), E_vec_cluster2.mag());
0265 pi0 = photon1 + photon2;
0266 h_diPhotonMassE->Fill(pi0.M(), pi0.E());
0267 }
0268 }
0269 }
0270
0271 TowerInfoContainer *ohcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ohcTowerNode);
0272 TowerInfoContainer *ihcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ihcTowerNode.c_str());
0273
0274 if (!ohcTowerContainer || !ihcTowerContainer)
0275 {
0276 std::cout << PHWHERE << "caloHistGen::process_event: " << m_ohcTowerNode << " or " << m_ohcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0277 return Fun4AllReturnCodes::ABORTEVENT;
0278 }
0279
0280 if (storeHCals && ohcTowerContainer && ihcTowerContainer)
0281 {
0282 tower_range = ohcTowerContainer->size();
0283 for (unsigned int iter = 0; iter < tower_range; iter++)
0284 {
0285 unsigned int towerkey = ohcTowerContainer->encode_key(iter);
0286 unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0287 unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0288 int time = ohcTowerContainer->get_tower_at_channel(iter)->get_time();
0289 float chi2 = ohcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0290 double energy = ohcTowerContainer->get_tower_at_channel(iter)->get_energy();
0291 float pedestal = ohcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0292
0293 if (ohcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0294 {
0295 h_OHCalTowE->Fill(ieta, iphi, energy);
0296 h_OHCalTowChi2->Fill(ieta, iphi, chi2);
0297 h_OHCalTowTime->Fill(ieta, iphi, time);
0298 h_OHCalTowPed->Fill(ieta, iphi, pedestal);
0299 }
0300 }
0301
0302 tower_range = ihcTowerContainer->size();
0303 for (unsigned int iter = 0; iter < tower_range; iter++)
0304 {
0305 unsigned int towerkey = ihcTowerContainer->encode_key(iter);
0306 unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0307 unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0308 int time = ihcTowerContainer->get_tower_at_channel(iter)->get_time();
0309 float chi2 = ihcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0310 double energy = ihcTowerContainer->get_tower_at_channel(iter)->get_energy();
0311 float pedestal = ihcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0312
0313 if (ihcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0314 {
0315 h_IHCalTowE->Fill(ieta, iphi, energy);
0316 h_IHCalTowChi2->Fill(ieta, iphi, chi2);
0317 h_IHCalTowTime->Fill(ieta, iphi, time);
0318 h_IHCalTowPed->Fill(ieta, iphi, pedestal);
0319 }
0320 }
0321 }
0322
0323 TowerInfoContainer *zdcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_zdcTowerNode.c_str());
0324 if (!zdcTowerContainer)
0325 {
0326 std::cout << PHWHERE << "caloHistGen::process_event: " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0327 return Fun4AllReturnCodes::ABORTEVENT;
0328 }
0329
0330 if (storeZDC && zdcTowerContainer)
0331 {
0332 tower_range = zdcTowerContainer->size();
0333 float totalZDCSouth = 0;
0334 float totalZDCNorth = 0;
0335 for (unsigned int iter = 0; iter < tower_range; iter++)
0336 {
0337 int time = zdcTowerContainer->get_tower_at_channel(iter)->get_time();
0338 if (iter < 13 && iter % 2 == 0)
0339 {
0340 h_zdcChanTime->Fill((int) iter / 2, time);
0341 }
0342
0343 float energy = zdcTowerContainer->get_tower_at_channel(iter)->get_energy();
0344 if (iter == 0 || iter == 2 || iter == 4)
0345 {
0346 totalZDCSouth += energy;
0347 }
0348 if (iter == 8 || iter == 10 || iter == 12)
0349 {
0350 totalZDCNorth += energy;
0351 }
0352 }
0353 h_zdcNSE->Fill(0., totalZDCSouth);
0354 h_zdcNSE->Fill(1., totalZDCNorth);
0355 }
0356
0357 return Fun4AllReturnCodes::EVENT_OK;
0358 }
0359
0360
0361 int caloHistGen::ResetEvent(PHCompositeNode * )
0362 {
0363 return Fun4AllReturnCodes::EVENT_OK;
0364 }
0365
0366
0367 int caloHistGen::End(PHCompositeNode * )
0368 {
0369 std::cout << "caloHistGen::End(PHCompositeNode *topNode) Saving Histos" << std::endl;
0370
0371 out->cd();
0372
0373 h_emcTowE-> Write();
0374 h_OHCalTowE-> Write();
0375 h_IHCalTowE-> Write();
0376
0377 h_emcTowChi2-> Write();
0378 h_OHCalTowChi2-> Write();
0379 h_IHCalTowChi2-> Write();
0380
0381 h_emcTowTime-> Write();
0382 h_OHCalTowTime-> Write();
0383 h_IHCalTowTime -> Write();
0384
0385 h_emcTowPed-> Write();
0386 h_OHCalTowPed->Write();
0387 h_IHCalTowPed->Write();
0388
0389 h_clusInfo->Write();
0390 h_diPhotonEtaPhiE->Write();
0391 h_diPhotonMassE->Write();
0392
0393 h_clusPt->Write();
0394 h_clusEcore->Write();
0395
0396 h_clusChi2->Write();
0397
0398 h_zVertex->Write();
0399
0400 h_zdcNSE->Write();
0401
0402 h_zdcChanTime->Write();
0403
0404 out->Close();
0405 delete out;
0406
0407 return Fun4AllReturnCodes::EVENT_OK;
0408 }
