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