Calorimeters/caloTreeGen/caloTreeGen.cc

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