emcalTreeMaker/src/caloTreeGen.cc

0001 #include "caloTreeGen.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006
0007 //Fun4All
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <fun4all/Fun4AllServer.h>
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>
0014 #include <ffaobjects/EventHeader.h>
0015
0016 //ROOT stuff
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019 #include <TH3F.h>
0020 #include <TFile.h>
0021 #include <TLorentzVector.h>
0022 #include <TTree.h>
0023
0024 //For clusters and geometry
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterContainer.h>
0027 #include <calobase/RawClusterUtility.h>
0028 #include <calobase/RawTowerGeomContainer.h>
0029
0030 //Tower stuff
0031 #include <calobase/TowerInfoContainer.h>
0032 #include <calobase/TowerInfo.h>
0033
0034 //thingy
0035 #include <CLHEP/Geometry/Point3D.h>
0036
0037
0038 //____________________________________________________________________________..
0039 caloTreeGen::caloTreeGen(const std::string &name):
0040 SubsysReco("caloTreeGen")
0041   ,T(nullptr)
0042   ,Outfile(name)
0043   ,doClusters(1)
0044 ,totalCaloE(0)
0045 ,doFineCluster(0)
0046 {
0047   std::cout << "caloTreeGen::caloTreeGen(const std::string &name) Calling ctor" << std::endl;
0048 }
0049
0050 //____________________________________________________________________________..
0051 caloTreeGen::~caloTreeGen()
0052 {
0053   std::cout << "caloTreeGen::~caloTreeGen() Calling dtor" << std::endl;
0054 }
0055
0056 //____________________________________________________________________________..
0057 int caloTreeGen::Init(PHCompositeNode *topNode)
0058 {
0059
0060   out = new TFile(Outfile.c_str(),"RECREATE");
0061
0062
0063   T = new TTree("T","T");
0064
0065   //emc
0066   T -> Branch("emcTowE",&m_emcTowE);
0067   T -> Branch("emcTowEta",&m_emcTowEta);
0068   T -> Branch("emcTowPhi",&m_emcTowPhi);
0069   T -> Branch("emcTiming",&m_emcTiming);
0070   T -> Branch("emciEta",&m_emciEta);
0071   T -> Branch("emciPhi",&m_emciPhi);
0072
0073   T -> Branch("clusterEFull",&m_clusterE);
0074   T -> Branch("clusterPhi",&m_clusterPhi);
0075   T -> Branch("clusterEta", &m_clusterEta);
0076   T -> Branch("clustrPt", &m_clusterPt);
0077   T -> Branch("clusterChi", &m_clusterChi);
0078   T -> Branch("clusterNtow",&m_clusterNtow);
0079   T -> Branch("clusterTowMax",&m_clusterTowMax);
0080   T -> Branch("clusterECore",&m_clusterECore);
0081
0082   T -> Branch("totalCaloE",&totalCaloE);
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  //so that the histos actually get written out
0089   Fun4AllServer *se = Fun4AllServer::instance();
0090   se -> Print("NODETREE");
0091
0092   std::cout << "caloTreeGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0093   return Fun4AllReturnCodes::EVENT_OK;
0094 }
0095
0096 //____________________________________________________________________________..
0097 int caloTreeGen::InitRun(PHCompositeNode *topNode)
0098 {
0099   std::cout << "caloTreeGen::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0100   return Fun4AllReturnCodes::EVENT_OK;
0101 }
0102
0103 //____________________________________________________________________________..
0104 int caloTreeGen::process_event(PHCompositeNode *topNode)
0105 {
0106
0107    //Information on clusters
0108   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_CEMC");
0109   if(!clusterContainer && doClusters)
0110     {
0111       std::cout << PHWHERE << "caloTreeGen::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0112       return 0;
0113     }
0114
0115   //tower information
0116   TowerInfoContainer *emcTowerContainer;
0117   emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC");
0118   if(!emcTowerContainer)
0119     {
0120       std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERS_CEMC"  << std::endl;
0121       return Fun4AllReturnCodes::ABORTEVENT;
0122     }
0123
0124
0125   //Tower geometry node for location information
0126   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0127   if (!towergeom && doClusters)
0128     {
0129       std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERGEOM_CEMC"  << std::endl;
0130       return Fun4AllReturnCodes::ABORTEVENT;
0131     }
0132
0133   float calib = 1.;
0134
0135   //grab all the towers and fill their energies.
0136   unsigned int tower_range = emcTowerContainer->size();
0137   totalCaloE = 0;
0138   for(unsigned int iter = 0; iter < tower_range; iter++)
0139     {
0140       unsigned int towerkey = emcTowerContainer->encode_key(iter);
0141       unsigned int ieta = getCaloTowerEtaBin(towerkey);
0142       unsigned int iphi = getCaloTowerPhiBin(towerkey);
0143       m_emciEta.push_back(ieta);
0144       m_emciPhi.push_back(iphi);
0145
0146       double energy = emcTowerContainer -> get_tower_at_channel(iter) -> get_energy()/calib;
0147       totalCaloE += energy;
0148       double time = emcTowerContainer -> get_tower_at_channel(iter) -> get_time();
0149
0150       if(doClusters)
0151     {
0152       double phi = towergeom -> get_phicenter(iphi);
0153       double eta = towergeom -> get_etacenter(ieta);
0154       m_emcTowPhi.push_back(phi);
0155       m_emcTowEta.push_back(eta);
0156     }
0157       m_emcTowE.push_back(energy);
0158       m_emcTiming.push_back(time);
0159
0160     }
0161
0162   if(doClusters)
0163     {
0164       RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
0165       RawClusterContainer::ConstIterator clusterIter;
0166       for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0167     {
0168       RawCluster *recoCluster = clusterIter -> second;
0169
0170       CLHEP::Hep3Vector vertex(0,0,0);
0171       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0172       CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*recoCluster, vertex);
0173
0174       float clusE = E_vec_cluster.mag();
0175       float clus_eta = E_vec_cluster.pseudoRapidity();
0176       float clus_phi = E_vec_cluster.phi();
0177       float clus_pt = E_vec_cluster.perp();
0178       float clus_chi = recoCluster -> get_chi2();
0179       float nTowers = recoCluster ->getNTowers();
0180       float maxTowerEnergy = getMaxTowerE(recoCluster,emcTowerContainer);
0181
0182       m_clusterE.push_back(E_vec_cluster_Full.mag());
0183       m_clusterECore.push_back(clusE);
0184       m_clusterPhi.push_back(clus_phi);
0185       m_clusterEta.push_back(clus_eta);
0186       m_clusterPt.push_back(clus_pt);
0187       m_clusterChi.push_back(clus_chi);
0188       m_clusterNtow.push_back(nTowers);
0189       m_clusterTowMax.push_back(maxTowerEnergy);
0190       if(doFineCluster)
0191         {
0192           m_clusTowPhi.push_back(returnClusterTowPhi(recoCluster,emcTowerContainer));
0193           m_clusTowEta.push_back(returnClusterTowEta(recoCluster,emcTowerContainer));
0194           m_clusTowE.push_back(returnClusterTowE(recoCluster,emcTowerContainer));
0195         }
0196     }
0197     }
0198
0199   T -> Fill();
0200
0201   return Fun4AllReturnCodes::EVENT_OK;
0202 }
0203
0204 //____________________________________________________________________________..
0205 int caloTreeGen::ResetEvent(PHCompositeNode *topNode)
0206 {
0207   m_clusterE.clear();
0208   m_clusterPhi.clear();
0209   m_clusterEta.clear();
0210   m_clusterPt.clear();
0211   m_clusterChi.clear();
0212   m_clusterTowMax.clear();
0213   m_clusterNtow.clear();
0214   m_clusterECore.clear();
0215
0216   m_emcTowPhi.clear();
0217   m_emcTowEta.clear();
0218   m_emcTowE.clear();
0219   m_emcTiming.clear();
0220   m_emciEta.clear();
0221   m_emciPhi.clear();
0222
0223   m_clusTowPhi.clear();
0224   m_clusTowEta.clear();
0225   m_clusTowE.clear();
0226
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229
0230 //____________________________________________________________________________..
0231 int caloTreeGen::EndRun(const int runnumber)
0232 {
0233   std::cout << "caloTreeGen::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0234   return Fun4AllReturnCodes::EVENT_OK;
0235 }
0236
0237 //____________________________________________________________________________..
0238 int caloTreeGen::End(PHCompositeNode *topNode)
0239 {
0240   std::cout << "caloTreeGen::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0241
0242   out -> cd();
0243   T -> Write();
0244   out -> Close();
0245   delete out;
0246
0247   //hm -> dumpHistos(Outfile.c_str(), "UPDATE");
0248   return Fun4AllReturnCodes::EVENT_OK;
0249 }
0250
0251 //____________________________________________________________________________..
0252 int caloTreeGen::Reset(PHCompositeNode *topNode)
0253 {
0254  std::cout << "caloTreeGen::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0255   return Fun4AllReturnCodes::EVENT_OK;
0256 }
0257
0258 //____________________________________________________________________________..
0259 void caloTreeGen::Print(const std::string &what) const
0260 {
0261   std::cout << "caloTreeGen::Print(const std::string &what) const Printing info for " << what << std::endl;
0262 }
0263 // convert from calorimeter key to phi bin
0264 //____________________________________________________________________________..
0265 unsigned int caloTreeGen::getCaloTowerPhiBin(const unsigned int key)
0266 {
0267   unsigned int etabin = key >> 16U;
0268   unsigned int phibin = key - (etabin << 16U);
0269   return phibin;
0270 }
0271
0272 //____________________________________________________________________________..
0273 // convert from calorimeter key to eta bin
0274 unsigned int caloTreeGen::getCaloTowerEtaBin(const unsigned int key)
0275 {
0276   unsigned int etabin = key >> 16U;
0277   return etabin;
0278 }
0279 //____________________________________________________________________________..
0280 float caloTreeGen::getMaxTowerE(RawCluster *cluster, TowerInfoContainer *towerContainer)
0281 {
0282   RawCluster::TowerConstRange towers = cluster -> get_towers();
0283   RawCluster::TowerConstIterator toweriter;
0284
0285   float maxEnergy = 0;
0286   for(toweriter = towers.first; toweriter != towers.second; toweriter++)
0287     {
0288       float towE = toweriter -> second;
0289
0290       if( towE > maxEnergy)  maxEnergy = towE;
0291     }
0292   return maxEnergy;
0293 }
0294 //____________________________________________________________________________..
0295 std::vector<int> caloTreeGen::returnClusterTowEta(RawCluster *cluster, TowerInfoContainer *towerContainer)
0296 {
0297   RawCluster::TowerConstRange towers = cluster -> get_towers();
0298   RawCluster::TowerConstIterator toweriter;
0299
0300   std::vector<int> towerIDsEta;
0301   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerIDsEta.push_back(RawTowerDefs::decode_index1(toweriter -> first));
0302
0303   return towerIDsEta;
0304 }
0305 //____________________________________________________________________________..
0306 std::vector<int> caloTreeGen::returnClusterTowPhi(RawCluster *cluster, TowerInfoContainer *towerContainer)
0307 {
0308   RawCluster::TowerConstRange towers = cluster -> get_towers();
0309   RawCluster::TowerConstIterator toweriter;
0310
0311   std::vector<int> towerIDsPhi;
0312   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerIDsPhi.push_back(RawTowerDefs::decode_index2(toweriter -> first));
0313   return towerIDsPhi;
0314 }
0315 //____________________________________________________________________________..
0316 std::vector<float> caloTreeGen::returnClusterTowE(RawCluster *cluster, TowerInfoContainer *towerContainer)
0317 {
0318   RawCluster::TowerConstRange towers = cluster -> get_towers();
0319   RawCluster::TowerConstIterator toweriter;
0320
0321   std::vector<float> towerE;
0322   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerE.push_back(toweriter -> second);
0323
0324   return towerE;
0325 }