Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:20:01 

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/getClass.h>
0025 #include <phool/phool.h>
0026 
0027 // ROOT stuff
0028 #include <TFile.h>
0029 #include <TH1.h>
0030 #include <TTree.h>
0031 
0032 // for cluster vertex correction
0033 #include <CLHEP/Vector/ThreeVector.h>
0034 
0035 #include <algorithm>
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   , Outfile(name)
0045 {
0046   std::cout << "caloTreeGen::caloTreeGen(const std::string &name) Calling ctor" << std::endl;
0047 }
0048 
0049 //____________________________________________________________________________..
0050 int caloTreeGen::Init(PHCompositeNode * /*topNode*/)
0051 {
0052   delete out;  // make cppcheck happy (nullptrs can be deleted)
0053   out = new TFile(Outfile.c_str(), "RECREATE");
0054 
0055   T = new TTree("T", "T");
0056 
0057   // Electromagnetic Calorimeter
0058   if (storeEMCal)
0059   {
0060     T->Branch("emcTowE", &m_emcTowE);
0061     T->Branch("emcTowiEta", &m_emciEta);
0062     T->Branch("emcTowiPhi", &m_emciPhi);
0063     T->Branch("emcTime", &m_emcTime);
0064     T->Branch("emcChi2", &m_emcChi2);
0065     T->Branch("emcPed", &m_emcPed);
0066 
0067     // EMCal Cluster information
0068     if (storeEMCal && storeClusters)
0069     {
0070       T->Branch("clusterE", &m_clusterE);
0071       T->Branch("clusterPhi", &m_clusterPhi);
0072       T->Branch("clusterEta", &m_clusterEta);
0073       T->Branch("clusterPt", &m_clusterPt);
0074       T->Branch("clusterChi2", &m_clusterChi);
0075       T->Branch("clusterNtow", &m_clusterNtow);
0076       T->Branch("clusterTowMaxE", &m_clusterTowMaxE);
0077       T->Branch("clusterECore", &m_clusterECore);
0078 
0079       // Information for towers within clusters
0080       // Enabled by setting "DoFineClusters" in the macro
0081       if (storeEMCal && storeClusters && storeClusterDetails)
0082       {
0083         T->Branch("clusTowPhi", "vector<vector<int> >", &m_clusTowPhi);
0084         T->Branch("clusTowEta", "vector<vector<int> >", &m_clusTowEta);
0085         T->Branch("clusTowE", "vector<vector<float> >", &m_clusTowE);
0086       }
0087     }
0088   }
0089   // Outer Hadronic Calorimeter
0090   if (storeHCals)
0091   {
0092     T->Branch("ohcTowE", &m_ohcTowE);
0093     T->Branch("ohcTowiEta", &m_ohciTowEta);
0094     T->Branch("ohcTowiPhi", &m_ohciTowPhi);
0095     T->Branch("ohcTime", &m_ohcTime);
0096     T->Branch("ohcChi2", &m_ohcChi2);
0097     T->Branch("ohcPed", &m_ohcPed);
0098 
0099     // Inner Hadronic Calorimeter
0100     T->Branch("ihcTowE", &m_ihcTowE);
0101     T->Branch("ihcTowiEta", &m_ihciTowEta);
0102     T->Branch("ihcTowiPhi", &m_ihciTowPhi);
0103     T->Branch("ihcTime", &m_ihcTime);
0104     T->Branch("ihcChi2", &m_ihcChi2);
0105     T->Branch("ihcPed", &m_ihcPed);
0106   }
0107   // ZDC information
0108   if (storeZDC)
0109   {
0110     T->Branch("zdcTowE", &m_zdcTowE);
0111     T->Branch("zdcTowside", &m_zdcSide);
0112 
0113     // SMD information
0114     T->Branch("smdE", &m_smdE);
0115     T->Branch("smdSide", &m_smdSide);
0116   }
0117   // Total
0118   T->Branch("totalCaloEEMCal", &totalCaloEEMCal);
0119   T->Branch("totalCaloEOHCal", &totalCaloEOHCal);
0120   T->Branch("totalCaloEIHCal", &totalCaloEIHCal);
0121   T->Branch("totalCaloEZDC", &totalCaloEZDC);
0122   T->Branch("zvertex", &m_vertex);
0123 
0124   // GL1 Information
0125   if (storeTrig)
0126   {
0127     T->Branch("triggerVector", &m_triggerVector);
0128   }
0129 
0130   zVertex = new TH1F("zVertex", "zVertex", 200, -100, 100);
0131 
0132   // so that the histos actually get written out
0133   Fun4AllServer *se = Fun4AllServer::instance();
0134   se->Print("NODETREE");
0135 
0136   std::cout << "caloTreeGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0137   return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139 
0140 //____________________________________________________________________________..
0141 int caloTreeGen::process_event(PHCompositeNode *topNode)
0142 {
0143   
0144   if(eventNum >= 100)
0145     {
0146       return Fun4AllReturnCodes::EVENT_OK;
0147     }
0148 
0149   m_vertex = -9999;
0150   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0151   if (!vertexmap)
0152   {
0153     std::cout << "GlobalVertexMap node is missing" << std::endl;
0154   }
0155   if (vertexmap && !vertexmap->empty())
0156   {
0157     GlobalVertex *vtx = vertexmap->begin()->second;
0158     if (vtx)
0159     {
0160       m_vertex = vtx->get_z();
0161       zVertex->Fill(m_vertex);
0162     }
0163   }
0164 
0165   // Information on clusters
0166   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_clusterNode);
0167   if (!clusterContainer && storeClusters)
0168   {
0169     std::cout << PHWHERE << "caloTreeGen::process_event: " << m_clusterNode << " node is missing. Output related to this node will be empty" << std::endl;
0170     return 0;
0171   }
0172 
0173   // tower information
0174   TowerInfoContainer *emcTowerContainer;
0175   emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode);
0176   if (!emcTowerContainer && storeEMCal)
0177   {
0178     std::cout << PHWHERE << "caloTreeGen::process_event: " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0179   }
0180 
0181   // grab all the towers and fill their energies.
0182   unsigned int tower_range = 0;
0183   if (storeEMCal && emcTowerContainer)
0184   {
0185     tower_range = emcTowerContainer->size();
0186     totalCaloEEMCal = 0;
0187     for (unsigned int iter = 0; iter < tower_range; iter++)
0188     {
0189       unsigned int towerkey = emcTowerContainer->encode_key(iter);
0190       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0191       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0192       double energy = emcTowerContainer->get_tower_at_channel(iter)->get_energy();
0193       if (energy < emcalThresh)
0194       {
0195         continue;
0196       }
0197       int time = emcTowerContainer->get_tower_at_channel(iter)->get_time();
0198       float chi2 = emcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0199       float pedestal = emcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0200       totalCaloEEMCal += energy;
0201       m_emciPhi.push_back(iphi);
0202       m_emciEta.push_back(ieta);
0203       m_emcTowE.push_back(energy);
0204       m_emcTime.push_back(time);
0205       m_emcChi2.push_back(chi2);
0206       m_emcPed.push_back(pedestal);
0207     }
0208   }
0209 
0210   if (storeClusters && storeEMCal)
0211   {
0212     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0213     RawClusterContainer::ConstIterator clusterIter;
0214     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0215     {
0216       RawCluster *recoCluster = clusterIter->second;
0217 
0218       CLHEP::Hep3Vector vertex(0, 0, 0);
0219       if (m_vertex != -9999)
0220       {
0221         vertex.setZ(m_vertex);
0222       }
0223       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0224       CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, vertex);
0225 
0226       float clusE = E_vec_cluster_Full.mag();
0227       if (clusE < clusterThresh)
0228       {
0229         continue;
0230       }
0231 
0232       float clusEcore = E_vec_cluster.mag();
0233       float clus_eta = E_vec_cluster.pseudoRapidity();
0234       float clus_phi = E_vec_cluster.phi();
0235       float clus_pt = E_vec_cluster.perp();
0236       float clus_chi = recoCluster->get_chi2();
0237       float nTowers = recoCluster->getNTowers();
0238       float maxTowerEnergy = getMaxTowerE(recoCluster);
0239 
0240       m_clusterE.push_back(clusE);
0241       m_clusterECore.push_back(clusEcore);
0242       m_clusterPhi.push_back(clus_phi);
0243       m_clusterEta.push_back(clus_eta);
0244       m_clusterPt.push_back(clus_pt);
0245       m_clusterChi.push_back(clus_chi);
0246       m_clusterNtow.push_back(nTowers);
0247       m_clusterTowMaxE.push_back(maxTowerEnergy);
0248       if (storeClusterDetails)
0249       {
0250         m_clusTowPhi.push_back(returnClusterTowPhi(recoCluster));
0251         m_clusTowEta.push_back(returnClusterTowEta(recoCluster));
0252         m_clusTowE.push_back(returnClusterTowE(recoCluster));
0253       }
0254     }
0255   }
0256 
0257   // tower information
0258   TowerInfoContainer *ohcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ohcTowerNode);
0259   TowerInfoContainer *ihcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_ihcTowerNode);
0260 
0261   if (!ohcTowerContainer || !ihcTowerContainer)
0262   {
0263     std::cout << PHWHERE << "caloTreeGen::process_event: " << m_ohcTowerNode << " or " << m_ohcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0264     return Fun4AllReturnCodes::ABORTEVENT;
0265   }
0266 
0267   if (storeHCals && ohcTowerContainer && ihcTowerContainer)
0268   {
0269     tower_range = ohcTowerContainer->size();
0270     totalCaloEOHCal = 0;
0271     for (unsigned int iter = 0; iter < tower_range; iter++)
0272     {
0273       unsigned int towerkey = ohcTowerContainer->encode_key(iter);
0274       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0275       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0276       int time = ohcTowerContainer->get_tower_at_channel(iter)->get_time();
0277       float chi2 = ohcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0278       double energy = ohcTowerContainer->get_tower_at_channel(iter)->get_energy();
0279       if (energy < ohcalThresh)
0280       {
0281         continue;
0282       }
0283       float pedestal = ohcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0284 
0285       totalCaloEOHCal += energy;
0286       m_ohcTowE.push_back(energy);
0287       m_ohciTowEta.push_back(ieta);
0288       m_ohciTowPhi.push_back(iphi);
0289       m_ohcTime.push_back(time);
0290       m_ohcChi2.push_back(chi2);
0291       m_ohcPed.push_back(pedestal);
0292     }
0293 
0294     tower_range = ihcTowerContainer->size();
0295     totalCaloEIHCal = 0;
0296     for (unsigned int iter = 0; iter < tower_range; iter++)
0297     {
0298       unsigned int towerkey = ihcTowerContainer->encode_key(iter);
0299       unsigned int ieta = TowerInfoDefs::getCaloTowerEtaBin(towerkey);
0300       unsigned int iphi = TowerInfoDefs::getCaloTowerPhiBin(towerkey);
0301       int time = ihcTowerContainer->get_tower_at_channel(iter)->get_time();
0302       float chi2 = ihcTowerContainer->get_tower_at_channel(iter)->get_chi2();
0303       double energy = ihcTowerContainer->get_tower_at_channel(iter)->get_energy();
0304       if (energy < ihcalThresh)
0305       {
0306         continue;
0307       }
0308       float pedestal = ihcTowerContainer->get_tower_at_channel(iter)->get_pedestal();
0309 
0310       totalCaloEIHCal += energy;
0311       m_ihcTowE.push_back(energy);
0312       m_ihciTowEta.push_back(ieta);
0313       m_ihciTowPhi.push_back(iphi);
0314       m_ihcTime.push_back(time);
0315       m_ihcChi2.push_back(chi2);
0316       m_ihcPed.push_back(pedestal);
0317     }
0318   }
0319 
0320   TowerInfoContainer *zdcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, m_zdcTowerNode);
0321   if (!zdcTowerContainer)
0322   {
0323     std::cout << PHWHERE << "caloTreeGen::process_event: " << m_emcTowerNode << " node is missing. Output related to this node will be empty" << std::endl;
0324     return Fun4AllReturnCodes::ABORTEVENT;
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);
0355   if (!gl1PacketInfo && storeTrig)
0356   {
0357     std::cout << PHWHERE << "caloTreeGen::process_event: " << 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     maxEnergy = std::max(towE, maxEnergy);
0442   }
0443   return maxEnergy;
0444 }
0445 //____________________________________________________________________________..
0446 std::vector<int> caloTreeGen::returnClusterTowEta(RawCluster *cluster)
0447 {
0448   RawCluster::TowerConstRange towers = cluster->get_towers();
0449   RawCluster::TowerConstIterator toweriter;
0450 
0451   std::vector<int> towerIDsEta;
0452   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0453   {
0454     towerIDsEta.push_back(RawTowerDefs::decode_index1(toweriter->first));
0455   }
0456 
0457   return towerIDsEta;
0458 }
0459 //____________________________________________________________________________..
0460 std::vector<int> caloTreeGen::returnClusterTowPhi(RawCluster *cluster)
0461 {
0462   RawCluster::TowerConstRange towers = cluster->get_towers();
0463   RawCluster::TowerConstIterator toweriter;
0464 
0465   std::vector<int> towerIDsPhi;
0466   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0467   {
0468     towerIDsPhi.push_back(RawTowerDefs::decode_index2(toweriter->first));
0469   }
0470   return towerIDsPhi;
0471 }
0472 //____________________________________________________________________________..
0473 std::vector<float> caloTreeGen::returnClusterTowE(RawCluster *cluster)
0474 {
0475   RawCluster::TowerConstRange towers = cluster->get_towers();
0476   RawCluster::TowerConstIterator toweriter;
0477 
0478   std::vector<float> towerE;
0479   for (toweriter = towers.first; toweriter != towers.second; toweriter++)
0480   {
0481     towerE.push_back(toweriter->second);
0482   }
0483 
0484   return towerE;
0485 }
[ 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!