Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:40 

0001 #include "caloTreeGen.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 <TTree.h>
0032 
0033 // for cluster vertex correction
0034 #include <CLHEP/Vector/ThreeVector.h>
0035 
0036 #include <cstdint>
0037 #include <iostream>
0038 #include <map>
0039 #include <utility>
0040 
0041 //____________________________________________________________________________..
0042 caloTreeGen::caloTreeGen(const std::string &name)
0043   : SubsysReco("CaloTreeGen")
0044   , T(nullptr)
0045   , Outfile(name)
0046 {
0047   std::cout << "caloTreeGen::caloTreeGen(const std::string &name) Calling ctor" << std::endl;
0048 }
0049 
0050 //____________________________________________________________________________..
0051 int caloTreeGen::Init(PHCompositeNode * /*topNode*/)
0052 {
0053   delete out;  // make cppcheck happy (nullptrs can be deleted)
0054   out = new TFile(Outfile.c_str(), "RECREATE");
0055 
0056   T = new TTree("T", "T");
0057 
0058   // Electromagnetic Calorimeter
0059   if (storeEMCal)
0060   {
0061     T->Branch("emcTowE", &m_emcTowE);
0062     T->Branch("emcTowiEta", &m_emciEta);
0063     T->Branch("emcTowiPhi", &m_emciPhi);
0064     T->Branch("emcTime", &m_emcTime);
0065     T->Branch("emcChi2", &m_emcChi2);
0066     T->Branch("emcPed", &m_emcPed);
0067 
0068     // EMCal Cluster information
0069     if (storeEMCal && storeClusters)
0070     {
0071       T->Branch("clusterE", &m_clusterE);
0072       T->Branch("clusterPhi", &m_clusterPhi);
0073       T->Branch("clusterEta", &m_clusterEta);
0074       T->Branch("clusterPt", &m_clusterPt);
0075       T->Branch("clusterChi2", &m_clusterChi);
0076       T->Branch("clusterNtow", &m_clusterNtow);
0077       T->Branch("clusterTowMaxE", &m_clusterTowMaxE);
0078       T->Branch("clusterECore", &m_clusterECore);
0079 
0080       // Information for towers within clusters
0081       // Enabled by setting "DoFineClusters" in the macro
0082       if (storeEMCal && storeClusters && storeClusterDetails)
0083       {
0084         T->Branch("clusTowPhi", "vector<vector<int> >", &m_clusTowPhi);
0085         T->Branch("clusTowEta", "vector<vector<int> >", &m_clusTowEta);
0086         T->Branch("clusTowE", "vector<vector<float> >", &m_clusTowE);
0087       }
0088     }
0089   }
0090   // Outer Hadronic Calorimeter
0091   if (storeHCals)
0092   {
0093     T->Branch("ohcTowE", &m_ohcTowE);
0094     T->Branch("ohcTowiEta", &m_ohciTowEta);
0095     T->Branch("ohcTowiPhi", &m_ohciTowPhi);
0096     T->Branch("ohcTime", &m_ohcTime);
0097     T->Branch("ohcChi2", &m_ohcChi2);
0098     T->Branch("ohcPed", &m_ohcPed);
0099 
0100     // Inner Hadronic Calorimeter
0101     T->Branch("ihcTowE", &m_ihcTowE);
0102     T->Branch("ihcTowiEta", &m_ihciTowEta);
0103     T->Branch("ihcTowiPhi", &m_ihciTowPhi);
0104     T->Branch("ihcTime", &m_ihcTime);
0105     T->Branch("ihcChi2", &m_ihcChi2);
0106     T->Branch("ihcPed", &m_ihcPed);
0107   }
0108   // ZDC information
0109   if (storeZDC)
0110   {
0111     T->Branch("zdcTowE", &m_zdcTowE);
0112     T->Branch("zdcTowside", &m_zdcSide);
0113 
0114     // SMD information
0115     T->Branch("smdE", &m_smdE);
0116     T->Branch("smdSide", &m_smdSide);
0117   }
0118   // Total
0119   T->Branch("totalCaloEEMCal", &totalCaloEEMCal);
0120   T->Branch("totalCaloEOHCal", &totalCaloEOHCal);
0121   T->Branch("totalCaloEIHCal", &totalCaloEIHCal);
0122   T->Branch("totalCaloEZDC", &totalCaloEZDC);
0123   T->Branch("zvertex", &m_vertex);
0124 
0125   // GL1 Information
0126   if (storeTrig)
0127   {
0128     T->Branch("triggerVector", &m_triggerVector);
0129   }
0130 
0131   zVertex = new TH1F("zVertex", "zVertex", 200, -100, 100);
0132 
0133   // so that the histos actually get written out
0134   Fun4AllServer *se = Fun4AllServer::instance();
0135   se->Print("NODETREE");
0136 
0137   std::cout << "caloTreeGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0138   return Fun4AllReturnCodes::EVENT_OK;
0139 }
0140 
0141 //____________________________________________________________________________..
0142 int caloTreeGen::process_event(PHCompositeNode *topNode)
0143 {
0144   
0145   if(eventNum >= 100)
0146     {
0147       return Fun4AllReturnCodes::EVENT_OK;
0148     }
0149 
0150   m_vertex = -9999;
0151   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0152   if (!vertexmap)
0153   {
0154     std::cout << "GlobalVertexMap node is missing" << std::endl;
0155   }
0156   if (vertexmap && !vertexmap->empty())
0157   {
0158     GlobalVertex *vtx = vertexmap->begin()->second;
0159     if (vtx)
0160     {
0161       m_vertex = vtx->get_z();
0162       zVertex->Fill(m_vertex);
0163     }
0164   }
0165 
0166   // Information on clusters
0167   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_clusterNode.c_str());
0168   if (!clusterContainer && storeClusters)
0169   {
0170     std::cout << PHWHERE << "caloTreeGen::process_event: Cluster " << m_clusterNode << " node is missing. Output related to this node will be empty" << std::endl;
0171     return 0;
0172   }
0173 
0174   // tower information
0175   TowerInfoContainer *emcTowerContainer;
0176   emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0177   if (!emcTowerContainer && storeEMCal)
0178   {
0179     std::cout << PHWHERE << "caloTreeGen::process_event: EMCal " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0180   }
0181 
0182   // grab all the towers and fill their energies.
0183   unsigned int tower_range = 0;
0184   if (storeEMCal && emcTowerContainer)
0185   {
0186     tower_range = emcTowerContainer->size();
0187     totalCaloEEMCal = 0;
0188     for (unsigned int iter = 0; iter < tower_range; iter++)
0189     {
0190       unsigned int towerkey = emcTowerContainer->encode_key(iter);
0191       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0192       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0193       double energy = emcTowerContainer->get_tower_at_channel(iter)->get_energy();
0194       if (energy < emcalThresh)
0195       {
0196         continue;
0197       }
0198       int time = emcTowerContainer->get_tower_at_channel(iter)->get_time();
0199       float chi2 = emcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0200       float pedestal = emcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0201       totalCaloEEMCal += energy;
0202       m_emciPhi.push_back(iphi);
0203       m_emciEta.push_back(ieta);
0204       m_emcTowE.push_back(energy);
0205       m_emcTime.push_back(time);
0206       m_emcChi2.push_back(chi2);
0207       m_emcPed.push_back(pedestal);
0208     }
0209   }
0210 
0211   if (storeClusters && storeEMCal)
0212   {
0213     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0214     RawClusterContainer::ConstIterator clusterIter;
0215     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0216     {
0217       RawCluster *recoCluster = clusterIter->second;
0218 
0219       CLHEP::Hep3Vector vertex(0, 0, 0);
0220       if (m_vertex != -9999)
0221       {
0222         vertex.setZ(m_vertex);
0223       }
0224       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0225       CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, vertex);
0226 
0227       float clusE = E_vec_cluster_Full.mag();
0228       if (clusE < clusterThresh)
0229       {
0230         continue;
0231       }
0232 
0233       float clusEcore = E_vec_cluster.mag();
0234       float clus_eta = E_vec_cluster.pseudoRapidity();
0235       float clus_phi = E_vec_cluster.phi();
0236       float clus_pt = E_vec_cluster.perp();
0237       float clus_chi = recoCluster->get_chi2();
0238       float nTowers = recoCluster->getNTowers();
0239       float maxTowerEnergy = getMaxTowerE(recoCluster);
0240 
0241       m_clusterE.push_back(clusE);
0242       m_clusterECore.push_back(clusEcore);
0243       m_clusterPhi.push_back(clus_phi);
0244       m_clusterEta.push_back(clus_eta);
0245       m_clusterPt.push_back(clus_pt);
0246       m_clusterChi.push_back(clus_chi);
0247       m_clusterNtow.push_back(nTowers);
0248       m_clusterTowMaxE.push_back(maxTowerEnergy);
0249       if (storeClusterDetails)
0250       {
0251         m_clusTowPhi.push_back(returnClusterTowPhi(recoCluster));
0252         m_clusTowEta.push_back(returnClusterTowEta(recoCluster));
0253         m_clusTowE.push_back(returnClusterTowE(recoCluster));
0254       }
0255     }
0256   }
0257 
0258   // tower information
0259   TowerInfoContainer *ohcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ohcTowerNode);
0260   TowerInfoContainer *ihcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ihcTowerNode.c_str());
0261 
0262   if ((!ohcTowerContainer || !ihcTowerContainer) && storeHCals)
0263   {
0264     std::cout << PHWHERE << "caloTreeGen::process_event: OHCal " << m_ohcTowerNode << " or IHCal " << m_ihcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0265     return Fun4AllReturnCodes::ABORTEVENT;
0266   }
0267 
0268   if (storeHCals && ohcTowerContainer && ihcTowerContainer)
0269   {
0270     tower_range = ohcTowerContainer->size();
0271     totalCaloEOHCal = 0;
0272     for (unsigned int iter = 0; iter < tower_range; iter++)
0273     {
0274       unsigned int towerkey = ohcTowerContainer->encode_key(iter);
0275       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0276       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0277       int time = ohcTowerContainer->get_tower_at_channel(iter)->get_time();
0278       float chi2 = ohcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0279       double energy = ohcTowerContainer->get_tower_at_channel(iter)->get_energy();
0280       if (energy < ohcalThresh)
0281       {
0282         continue;
0283       }
0284       float pedestal = ohcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0285 
0286       totalCaloEOHCal += energy;
0287       m_ohcTowE.push_back(energy);
0288       m_ohciTowEta.push_back(ieta);
0289       m_ohciTowPhi.push_back(iphi);
0290       m_ohcTime.push_back(time);
0291       m_ohcChi2.push_back(chi2);
0292       m_ohcPed.push_back(pedestal);
0293     }
0294 
0295     tower_range = ihcTowerContainer->size();
0296     totalCaloEIHCal = 0;
0297     for (unsigned int iter = 0; iter < tower_range; iter++)
0298     {
0299       unsigned int towerkey = ihcTowerContainer->encode_key(iter);
0300       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0301       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0302       int time = ihcTowerContainer->get_tower_at_channel(iter)->get_time();
0303       float chi2 = ihcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0304       double energy = ihcTowerContainer->get_tower_at_channel(iter)->get_energy();
0305       if (energy < ihcalThresh)
0306       {
0307         continue;
0308       }
0309       float pedestal = ihcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0310 
0311       totalCaloEIHCal += energy;
0312       m_ihcTowE.push_back(energy);
0313       m_ihciTowEta.push_back(ieta);
0314       m_ihciTowPhi.push_back(iphi);
0315       m_ihcTime.push_back(time);
0316       m_ihcChi2.push_back(chi2);
0317       m_ihcPed.push_back(pedestal);
0318     }
0319   }
0320 
0321   TowerInfoContainer *zdcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_zdcTowerNode.c_str());
0322   if (!zdcTowerContainer && storeZDC)
0323   {
0324     std::cout << PHWHERE << "caloTreeGen::process_event: ZDC " << m_zdcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0325   }
0326 
0327   if (storeZDC && zdcTowerContainer)
0328   {
0329     tower_range = zdcTowerContainer->size();
0330     totalCaloEZDC = 0;
0331     for (unsigned int iter = 0; iter < tower_range; iter++)
0332     {
0333       if (iter < 16)
0334       {
0335         float energy = zdcTowerContainer->get_tower_at_channel(iter)->get_energy();
0336         m_zdcTowE.push_back(energy);
0337         unsigned int towerkey = zdcTowerContainer->encode_key(iter);
0338         unsigned int side = TowerInfoDefs::get_zdc_side(towerkey);
0339         m_zdcSide.push_back(side);
0340         totalCaloEZDC += energy;
0341       }
0342       if (iter > 15 && iter < 48)
0343       {
0344         // smd north stuff
0345         float energy = zdcTowerContainer->get_tower_at_channel(iter)->get_energy();
0346         m_smdE.push_back(energy);
0347         unsigned int towerkey = zdcTowerContainer->encode_key(iter);
0348         unsigned int side = TowerInfoDefs::get_zdc_side(towerkey);
0349         m_smdSide.push_back(side);
0350       }
0351     }
0352   }
0353 
0354   Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, m_trigNode.c_str());
0355   if (!gl1PacketInfo && storeTrig)
0356   {
0357     std::cout << PHWHERE << "caloTreeGen::process_event: Gl1 " << m_trigNode << " node is missing. Output related to this node will be empty" << std::endl;
0358   }
0359 
0360   if (storeTrig && gl1PacketInfo)
0361   {
0362     uint64_t triggervec = gl1PacketInfo->getScaledVector();
0363     for (int i = 0; i < 64; i++)
0364     {
0365       bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0366       m_triggerVector.push_back(trig_decision);
0367       triggervec = (triggervec >> 1U) & 0xffffffffU;
0368     }
0369   }
0370   T->Fill();
0371   eventNum++;
0372 
0373   return Fun4AllReturnCodes::EVENT_OK;
0374 }
0375 
0376 //____________________________________________________________________________..
0377 int caloTreeGen::ResetEvent(PHCompositeNode * /*topNode*/)
0378 {
0379   m_clusterE.clear();
0380   m_clusterPhi.clear();
0381   m_clusterEta.clear();
0382   m_clusterPt.clear();
0383   m_clusterChi.clear();
0384   m_clusterTowMaxE.clear();
0385   m_clusterNtow.clear();
0386   m_clusterECore.clear();
0387 
0388   m_emcTowE.clear();
0389   m_emciEta.clear();
0390   m_emciPhi.clear();
0391   m_emcTime.clear();
0392   m_emcChi2.clear();
0393   m_emcPed.clear();
0394 
0395   m_ihcTowE.clear();
0396   m_ihciTowEta.clear();
0397   m_ihciTowPhi.clear();
0398   m_ihcTime.clear();
0399   m_ihcChi2.clear();
0400   m_ihcPed.clear();
0401 
0402   m_ohcTowE.clear();
0403   m_ohciTowEta.clear();
0404   m_ohciTowPhi.clear();
0405   m_ohcTime.clear();
0406   m_ohcChi2.clear();
0407   m_ohcPed.clear();
0408 
0409   m_clusTowPhi.clear();
0410   m_clusTowEta.clear();
0411   m_clusTowE.clear();
0412 
0413   return Fun4AllReturnCodes::EVENT_OK;
0414 }
0415 
0416 //____________________________________________________________________________..
0417 int caloTreeGen::End(PHCompositeNode * /*topNode*/)
0418 {
0419   std::cout << "caloTreeGen::End(PHCompositeNode *topNode) Saving TTree" << std::endl;
0420 
0421   out->cd();
0422   T->Write();
0423   zVertex->Write();
0424   out->Close();
0425   delete out;
0426 
0427   return Fun4AllReturnCodes::EVENT_OK;
0428 }
0429 
0430 //____________________________________________________________________________..
0431 float caloTreeGen::getMaxTowerE(RawCluster *cluster)
0432 {
0433   RawCluster::TowerConstRange towers = cluster->get_towers();
0434   RawCluster::TowerConstIterator toweriter;
0435 
0436   float maxEnergy = 0;
0437   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0438   {
0439     float towE = toweriter->second;
0440 
0441     if (towE > maxEnergy)
0442     {
0443       maxEnergy = towE;
0444     }
0445   }
0446   return maxEnergy;
0447 }
0448 //____________________________________________________________________________..
0449 std::vector<int> caloTreeGen::returnClusterTowEta(RawCluster *cluster)
0450 {
0451   RawCluster::TowerConstRange towers = cluster->get_towers();
0452   RawCluster::TowerConstIterator toweriter;
0453 
0454   std::vector<int> towerIDsEta;
0455   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0456   {
0457     towerIDsEta.push_back(RawTowerDefs::decode_index1(toweriter->first));
0458   }
0459 
0460   return towerIDsEta;
0461 }
0462 //____________________________________________________________________________..
0463 std::vector<int> caloTreeGen::returnClusterTowPhi(RawCluster *cluster)
0464 {
0465   RawCluster::TowerConstRange towers = cluster->get_towers();
0466   RawCluster::TowerConstIterator toweriter;
0467 
0468   std::vector<int> towerIDsPhi;
0469   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0470   {
0471     towerIDsPhi.push_back(RawTowerDefs::decode_index2(toweriter->first));
0472   }
0473   return towerIDsPhi;
0474 }
0475 //____________________________________________________________________________..
0476 std::vector<float> caloTreeGen::returnClusterTowE(RawCluster *cluster)
0477 {
0478   RawCluster::TowerConstRange towers = cluster->get_towers();
0479   RawCluster::TowerConstIterator toweriter;
0480 
0481   std::vector<float> towerE;
0482   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0483   {
0484     towerE.push_back(toweriter->second);
0485   }
0486 
0487   return towerE;
0488 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!