CaloDataAnaRun24pp/src/caloHistGen.cc

0001 #include "caloHistGen.h"
0002
0003 // For clusters and geometry
0004 #include <calobase/RawCluster.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawClusterUtility.h>
0007 #include <calobase/RawTowerDefs.h>
0008
0009 // Tower stuff
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012 #include <calobase/TowerInfoDefs.h>
0013
0014 // for the vertex
0015 #include <globalvertex/GlobalVertex.h>
0016 #include <globalvertex/GlobalVertexMap.h>
0017
0018 // GL1 Information
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 // ROOT stuff
0029 #include <TFile.h>
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TH3.h>
0033 #include <TTree.h>
0034
0035 // for cluster vertex correction
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 * /*topNode*/)
0053 {
0054   delete out;  // make cppcheck happy (nullptrs can be deleted)
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, 255.5, 100, -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, 100, -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, 100, -10, 25);
0060
0061   h_emcTowChi2 = new TH3F("emcTowChi2", "tower eta, tower phi, tower Chi2", 96, -0.5, 95.5, 256, -0.5, 255.5, 100, 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, 100, 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, 100, 0.5, 4e8);
0064
0065   h_emcTowTime = new TH3F("emcTowEnergy", "tower eta, tower phi, tower Time", 96, -0.5, 95.5, 256, -0.5, 255.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, 255.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", 100, -2, 25);
0076   h_clusEcore = new TH1F("clusEcore", "cluster Ecore", 100, -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", 50, -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", 50, 0.02, 1, 40, 0, 20);
0089
0090   // so that the histos actually get written out
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   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0103   if (!vertexmap)
0104   {
0105     std::cout << "GlobalVertexMap node is missing" << std::endl;
0106   }
0107   if (vertexmap && !vertexmap->empty())
0108   {
0109     GlobalVertex *vtx = vertexmap->begin()->second;
0110     if (vtx)
0111     {
0112       vertex = vtx->get_z();
0113       h_zVertex->Fill(vertex);
0114     }
0115   }
0116
0117   Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, m_trigNode.c_str());
0118   if (!gl1PacketInfo && checkTrig)
0119   {
0120     std::cout << PHWHERE << "caloHistGen::process_event: " << m_trigNode << " node is missing. Output related to this node will be empty" << std::endl;
0121   }
0122
0123   if (gl1PacketInfo && checkTrig)
0124   {
0125     uint64_t triggervec = gl1PacketInfo->getScaledVector();
0126     for (int i = 0; i < 64; i++)
0127     {
0128       bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0129
0130       triggervec = (triggervec >> 1U) & 0xffffffffU;
0131       if (!trig_decision && trigRequired[i])
0132       {
0133         std::cout << "Trigger check failed, skipping event" << std::endl;
0134         return 0;
0135       }
0136     }
0137   }
0138
0139   // Information on clusters
0140   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_clusterNode.c_str());
0141   if (!clusterContainer && storeClusters)
0142   {
0143     std::cout << PHWHERE << "caloHistGen::process_event: " << m_clusterNode << " node is missing. Output related to this node will be empty" << std::endl;
0144     return 0;
0145   }
0146
0147   // tower information
0148   TowerInfoContainer *emcTowerContainer;
0149   emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0150   if (!emcTowerContainer && storeEMCal)
0151   {
0152     std::cout << PHWHERE << "caloHistGen::process_event: " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0153   }
0154
0155   // grab all the towers and fill their energies.
0156   unsigned int tower_range = 0;
0157   if (storeEMCal && emcTowerContainer)
0158   {
0159     tower_range = emcTowerContainer->size();
0160     for (unsigned int iter = 0; iter < tower_range; iter++)
0161     {
0162       unsigned int towerkey = emcTowerContainer->encode_key(iter);
0163       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0164       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0165       double energy = emcTowerContainer->get_tower_at_channel(iter)->get_energy();
0166       int time = emcTowerContainer->get_tower_at_channel(iter)->get_time();
0167       float chi2 = emcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0168       float pedestal = emcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0169
0170       if (emcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0171       {
0172         h_emcTowE->Fill(ieta, iphi, energy);
0173         h_emcTowChi2->Fill(ieta, iphi, chi2);
0174         h_emcTowTime->Fill(ieta, iphi, time);
0175         h_emcTowPed->Fill(ieta, iphi, pedestal);
0176       }
0177     }
0178   }
0179
0180   if (storeClusters && storeEMCal)
0181   {
0182     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0183     RawClusterContainer::ConstIterator clusterIter;
0184     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0185     {
0186       RawCluster *recoCluster = clusterIter->second;
0187
0188       CLHEP::Hep3Vector hep_vertex(0, 0, 0);
0189       if (vertex != -9999)
0190       {
0191         hep_vertex.setZ(vertex);
0192       }
0193       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, hep_vertex);
0194       CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, hep_vertex);
0195
0196       float clusE = E_vec_cluster_Full.mag();
0197       float clusEcore = E_vec_cluster.mag();
0198       float clus_eta = E_vec_cluster.pseudoRapidity();
0199       float clus_phi = E_vec_cluster.phi();
0200       float clus_pt = E_vec_cluster.perp();
0201
0202       h_clusInfo->Fill(clus_eta, clus_phi, clusE);
0203       h_clusPt->Fill(clus_pt);
0204       h_clusEcore->Fill(clusEcore);
0205     }
0206   }
0207
0208   // pi0 reconstruction
0209   float caloEnergy =  getTotalCaloEnergy(emcTowerContainer);
0210   bool doHIPi0Reco = true;
0211   if(isAuAu)
0212   {
0213     if(!(peripheralOnly &&  caloEnergy < caloFrac * emcaldownscale)) doHIPi0Reco = false;
0214   }
0215
0216   if (doPi0Reco && storeEMCal && doHIPi0Reco)
0217   {
0218     RawClusterContainer::ConstRange clusters = clusterContainer->getClusters();
0219     RawClusterContainer::ConstIterator clusterIter1, clusterIter2;
0220     int clusterCounter1 = 0;
0221     int clusterCounter2 = 0;
0222
0223     for (clusterIter1 = clusters.first; clusterIter1 != clusters.second; clusterIter1++)
0224     {
0225       RawCluster *recoCluster1 = clusterIter1->second;
0226       clusterCounter1++;
0227       if (recoCluster1->get_chi2() > 4)
0228       {
0229         continue;
0230       }
0231
0232       CLHEP::Hep3Vector hep_vertex(0, 0, 0);
0233       if (vertex != -9999)
0234       {
0235         hep_vertex.setZ(vertex);
0236       }
0237       CLHEP::Hep3Vector E_vec_cluster1 = RawClusterUtility::GetECoreVec(*recoCluster1, hep_vertex);
0238
0239       for (clusterIter2 = clusters.first; clusterIter2 != clusters.second; clusterIter2++)
0240       {
0241     clusterCounter2++;
0242         if (clusterCounter2 <= clusterCounter1)
0243         {
0244           continue;  // prevents double counting pairs
0245         }
0246         RawCluster *recoCluster2 = clusterIter2->second;
0247         CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, hep_vertex);
0248         if (recoCluster2->get_chi2() > 4)
0249         {
0250           continue;
0251         }
0252     float clusterLead = E_vec_cluster1.mag();
0253     float clusterSub = E_vec_cluster2.mag();
0254     if(E_vec_cluster2.mag() > clusterLead)
0255       {
0256         clusterLead = E_vec_cluster2.mag();
0257         clusterSub = E_vec_cluster1.mag();
0258       }
0259
0260         if (clusterLead < clus1EMin || clusterSub < clus2EMin)
0261         {
0262           continue;
0263         }
0264         if (std::fabs(E_vec_cluster1.mag() - E_vec_cluster2.mag()) / (E_vec_cluster1.mag() + E_vec_cluster2.mag()) > maxAlpha)
0265         {
0266           continue;
0267         }
0268         TLorentzVector photon1, photon2, pi0;
0269         photon1.SetPtEtaPhiE(E_vec_cluster1.perp(), E_vec_cluster1.pseudoRapidity(), E_vec_cluster1.phi(), E_vec_cluster1.mag());
0270         photon2.SetPtEtaPhiE(E_vec_cluster2.perp(), E_vec_cluster2.pseudoRapidity(), E_vec_cluster2.phi(), E_vec_cluster2.mag());
0271         pi0 = photon1 + photon2;
0272         h_diPhotonMassE->Fill(pi0.M(), pi0.E());
0273       }
0274     }
0275   }
0276
0277   // tower information
0278   TowerInfoContainer *ohcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ohcTowerNode);
0279   TowerInfoContainer *ihcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ihcTowerNode.c_str());
0280
0281   if (!ohcTowerContainer || !ihcTowerContainer)
0282   {
0283     std::cout << PHWHERE << "caloHistGen::process_event: " << m_ohcTowerNode << " or " << m_ohcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0284     return Fun4AllReturnCodes::ABORTEVENT;
0285   }
0286
0287   if (storeHCals && ohcTowerContainer && ihcTowerContainer)
0288   {
0289     tower_range = ohcTowerContainer->size();
0290     for (unsigned int iter = 0; iter < tower_range; iter++)
0291     {
0292       unsigned int towerkey = ohcTowerContainer->encode_key(iter);
0293       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0294       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0295       int time = ohcTowerContainer->get_tower_at_channel(iter)->get_time();
0296       float chi2 = ohcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0297       double energy = ohcTowerContainer->get_tower_at_channel(iter)->get_energy();
0298       float pedestal = ohcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0299
0300       if (ohcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0301       {
0302         h_OHCalTowE->Fill(ieta, iphi, energy);
0303         h_OHCalTowChi2->Fill(ieta, iphi, chi2);
0304         h_OHCalTowTime->Fill(ieta, iphi, time);
0305         h_OHCalTowPed->Fill(ieta, iphi, pedestal);
0306       }
0307     }
0308
0309     tower_range = ihcTowerContainer->size();
0310     for (unsigned int iter = 0; iter < tower_range; iter++)
0311     {
0312       unsigned int towerkey = ihcTowerContainer->encode_key(iter);
0313       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0314       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0315       int time = ihcTowerContainer->get_tower_at_channel(iter)->get_time();
0316       float chi2 = ihcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0317       double energy = ihcTowerContainer->get_tower_at_channel(iter)->get_energy();
0318       float pedestal = ihcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0319
0320       if (ihcTowerContainer->get_tower_at_channel(iter)->get_isGood())
0321       {
0322         h_IHCalTowE->Fill(ieta, iphi, energy);
0323         h_IHCalTowChi2->Fill(ieta, iphi, chi2);
0324         h_IHCalTowTime->Fill(ieta, iphi, time);
0325         h_IHCalTowPed->Fill(ieta, iphi, pedestal);
0326       }
0327     }
0328   }
0329
0330   TowerInfoContainer *zdcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_zdcTowerNode.c_str());
0331   if (!zdcTowerContainer)
0332   {
0333     std::cout << PHWHERE << "caloHistGen::process_event: " << m_zdcTowerNode.c_str() << " node is missing. Output related to this node will be empty" << std::endl;
0334     return Fun4AllReturnCodes::ABORTEVENT;
0335   }
0336
0337   if (storeZDC && zdcTowerContainer)
0338   {
0339     tower_range = zdcTowerContainer->size();
0340     float totalZDCSouth = 0;
0341     float totalZDCNorth = 0;
0342     for (unsigned int iter = 0; iter < tower_range; iter++)
0343     {
0344       int time = zdcTowerContainer->get_tower_at_channel(iter)->get_time();
0345       if (iter < 13 && iter % 2 == 0)
0346       {
0347         h_zdcChanTime->Fill((int) iter / 2, time);
0348       }
0349
0350       float energy = zdcTowerContainer->get_tower_at_channel(iter)->get_energy();
0351       if (iter == 0 || iter == 2 || iter == 4)
0352       {
0353         totalZDCSouth += energy;
0354       }
0355       if (iter == 8 || iter == 10 || iter == 12)
0356       {
0357         totalZDCNorth += energy;
0358       }
0359     }
0360     h_zdcNSE->Fill(0., totalZDCSouth);
0361     h_zdcNSE->Fill(1., totalZDCNorth);
0362   }
0363
0364   return Fun4AllReturnCodes::EVENT_OK;
0365 }
0366
0367 //____________________________________________________________________________..
0368 int caloHistGen::ResetEvent(PHCompositeNode * /*topNode*/)
0369 {
0370   return Fun4AllReturnCodes::EVENT_OK;
0371 }
0372
0373 //____________________________________________________________________________..
0374 int caloHistGen::End(PHCompositeNode * /*topNode*/)
0375 {
0376   std::cout << "caloHistGen::End(PHCompositeNode *topNode) Saving Histos" << std::endl;
0377
0378   out->cd();
0379
0380   h_emcTowE-> Write();
0381   h_OHCalTowE-> Write();
0382   h_IHCalTowE-> Write();
0383
0384   h_emcTowChi2-> Write();
0385   h_OHCalTowChi2-> Write();
0386   h_IHCalTowChi2-> Write();
0387
0388   h_emcTowTime-> Write();
0389   h_OHCalTowTime-> Write();
0390   h_IHCalTowTime -> Write();
0391
0392   h_emcTowPed-> Write();
0393   h_OHCalTowPed->Write();
0394   h_IHCalTowPed->Write();
0395
0396   h_clusInfo->Write();
0397   h_diPhotonEtaPhiE->Write();
0398   h_diPhotonMassE->Write();
0399
0400   h_clusPt->Write();
0401   h_clusEcore->Write();
0402
0403   h_clusChi2->Write();
0404
0405   h_zVertex->Write();
0406
0407   h_zdcNSE->Write();
0408
0409   h_zdcChanTime->Write();
0410
0411   out->Close();
0412   delete out;
0413
0414   return Fun4AllReturnCodes::EVENT_OK;
0415 }
0416
0417 float caloHistGen::getTotalCaloEnergy(TowerInfoContainer *towerContainer)
0418 {
0419
0420   float totalCaloEnergy = 0;
0421
0422   if(!towerContainer) return std::nanf("1");
0423
0424
0425   int tower_range = towerContainer->size();
0426   for(int iter = 0; iter < tower_range; iter++)
0427     {
0428       if(towerContainer->get_tower_at_channel(iter) -> get_isGood())
0429     {
0430       double energy = towerContainer->get_tower_at_channel(iter)->get_energy();
0431       totalCaloEnergy += energy;
0432     }
0433     }
0434
0435   return totalCaloEnergy;
0436 }