isoCluster/src/isoCluster.cc

0001
0002 //____________________________________________________________________________..
0003
0004 //our base code
0005 #include "isoCluster.h"
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010
0011 //Fun4all stuff
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/Fun4AllHistoManager.h>
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/getClass.h>
0017 #include <phool/phool.h>
0018 #include <ffaobjects/EventHeader.h>
0019
0020 //ROOT stuff
0021 #include <TH1F.h>
0022 #include <TH2F.h>
0023 #include <TH3F.h>
0024 #include <TFile.h>
0025 #include <TLorentzVector.h>
0026 #include <TTree.h>
0027
0028 //for emc clusters
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterContainer.h>
0031 #include <calobase/RawClusterUtility.h>
0032 #include <calobase/RawTowerGeomContainer.h>
0033 #include <calobase/RawTower.h>
0034 #include <calobase/RawTowerContainer.h>
0035
0036 //caloEvalStack for cluster to truth matching
0037 #include <g4eval/CaloEvalStack.h>
0038 #include <g4eval/CaloRawClusterEval.h>
0039
0040 //for vetex information
0041 #include <g4vertex/GlobalVertex.h>
0042 #include <g4vertex/GlobalVertexMap.h>
0043
0044 //tracking
0045 #include <trackbase_historic/SvtxTrack.h>
0046 #include <trackbase_historic/SvtxTrackMap.h>
0047 #include <trackbase_historic/SvtxVertex.h>
0048 #include <trackbase_historic/SvtxVertexMap.h>
0049
0050 //truth information
0051 #include <g4main/PHG4TruthInfoContainer.h>
0052 #include <g4main/PHG4Particle.h>
0053 #include <g4main/PHG4VtxPoint.h>
0054 #include <phhepmc/PHHepMCGenEvent.h>
0055 #include <phhepmc/PHHepMCGenEventMap.h>
0056 #pragma GCC diagnostic push
0057 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0058 #include <HepMC/GenEvent.h>
0059 #include <HepMC/GenParticle.h>
0060 #include <HepMC/GenVertex.h>
0061 #include <HepMC/IteratorRange.h>
0062 #include <HepMC/SimpleVector.h>
0063 #include <HepMC/GenParticle.h>
0064 #pragma GCC diagnostic pop
0065
0066 //____________________________________________________________________________..
0067 isoCluster::isoCluster(const std::string &name):
0068 SubsysReco(name),
0069   m_clusterPhi(0),
0070 m_clusterEta(0),
0071   m_clusterE(0),
0072 m_clusterPt(0),
0073   m_clusterEtIso(0),
0074   m_clusterChisq(0),
0075   m_photonPhi(0),
0076   m_photonEta(0),
0077   m_photonE(0),
0078 m_photonPt(0),
0079 fout(NULL),
0080   outname(name),
0081   getEvent(-9999)
0082 {
0083 std::cout << "isoCluster::isoCluster(const std::string &name) Calling ctor" << std::endl;
0084 }
0085 //____________________________________________________________________________..
0086 isoCluster::~isoCluster()
0087 {
0088   std::cout << "isoCluster::~isoCluster() Calling dtor" << std::endl;
0089 }
0090
0091 //____________________________________________________________________________..
0092 int isoCluster::Init(PHCompositeNode *topNode)
0093 {
0094   std::cout << "isoCluster::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0095
0096   fout = new TFile(outname.c_str(),"RECREATE");
0097
0098   T = new TTree("T","T");
0099
0100   T -> Branch("clusterPhi",&m_clusterPhi);
0101   T -> Branch("clusterEta",&m_clusterEta);
0102   T -> Branch("clusterE",&m_clusterE);
0103   T -> Branch("clusterPt",&m_clusterPt);
0104   T -> Branch("clusterEtIso",&m_clusterEtIso);
0105   T -> Branch("clusterchisq",&m_clusterChisq);
0106
0107   T -> Branch("photonPhi",&m_photonPhi);
0108   T -> Branch("photonEta",&m_photonEta);
0109   T -> Branch("photonE",&m_photonE);
0110   T -> Branch("photonPt",&m_photonPt);
0111
0112   return Fun4AllReturnCodes::EVENT_OK;
0113 }
0114
0115 //____________________________________________________________________________..
0116 int isoCluster::InitRun(PHCompositeNode *topNode)
0117 {
0118   std::cout << "isoCluster::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0119   return Fun4AllReturnCodes::EVENT_OK;
0120 }
0121
0122 //____________________________________________________________________________..
0123 int isoCluster::process_event(PHCompositeNode *topNode)
0124 {
0125   //Information on clusters
0126   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0127   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0128   if(!clusterContainer)
0129     {
0130       std::cout << PHWHERE << "isoCluster::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
0131
0132       return 0;
0133     }
0134
0135   //Vertex information
0136   GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0137   if (!vtxContainer)
0138     {
0139       std::cout << PHWHERE << "isoCluster::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0140       assert(vtxContainer);  // force quit
0141
0142       return 0;
0143     }
0144
0145   if (vtxContainer->empty())
0146     {
0147       std::cout << PHWHERE << "isoCluster::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0148       return 0;
0149     }
0150
0151   //More vertex information
0152   GlobalVertex *vtx = vtxContainer->begin()->second;
0153   if(!vtx)
0154     {
0155
0156       std::cout << PHWHERE << "isoCluster::process_event Could not find vtx from vtxContainer"  << std::endl;
0157       return Fun4AllReturnCodes::ABORTEVENT;
0158     }
0159
0160   //truth particle information
0161   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0162   if(!truthinfo)
0163     {
0164       std::cout << PHWHERE << "isoCluster::process_event Could not find node G4TruthInfo"  << std::endl;
0165       return Fun4AllReturnCodes::ABORTEVENT;
0166     }
0167
0168   //For eventgen ancestory information
0169   PHHepMCGenEventMap *genEventMap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0170   if(!genEventMap)
0171     {
0172       std::cout << PHWHERE << "isoCluster::process_event Could not find PHHepMCGenEventMap"  << std::endl;
0173       return Fun4AllReturnCodes::ABORTEVENT;
0174     }
0175
0176   //event level information of the above
0177   PHHepMCGenEvent *genEvent = genEventMap -> get(getEvent);
0178   if(!genEvent)
0179     {
0180       std::cout << PHWHERE << "isoCluster::process_event Could not find PHHepMCGenEvent"  << std::endl;
0181       return Fun4AllReturnCodes::ABORTEVENT;
0182     }
0183
0184   HepMC::GenEvent *theEvent = genEvent -> getEvent();
0185
0186   CaloEvalStack caloevalstack(topNode, "CEMC");
0187   CaloRawClusterEval *clustereval = caloevalstack.get_rawcluster_eval();
0188   if(!clustereval)
0189     {
0190       std::cout << PHWHERE << "isoCluster::process_event cluster eval not available"  << std::endl;
0191       return Fun4AllReturnCodes::ABORTEVENT;
0192     }
0193
0194   //grab all the clusters in the event
0195   RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
0196   RawClusterContainer::ConstIterator clusterIter;
0197   for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0198     {
0199       RawCluster *recoCluster = clusterIter -> second;
0200       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0201       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0202       if(E_vec_cluster.mag() < 3) continue;
0203       m_clusterE.push_back(E_vec_cluster.mag());
0204       m_clusterEta.push_back(E_vec_cluster.pseudoRapidity());
0205       m_clusterPhi.push_back(E_vec_cluster.phi());
0206       m_clusterPt.push_back(E_vec_cluster.perp());
0207       m_clusterChisq.push_back(recoCluster -> get_chi2());
0208       m_clusterEtIso.push_back(recoCluster -> get_et_iso(3,1,1));
0209     }
0210
0211
0212   //grab all the direct photons in the event
0213   PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
0214   PHG4TruthInfoContainer::ConstIterator truthIter;
0215   //from the HepMC event log
0216   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0217     {
0218       //PHG4Particle* part = truthinfo->GetParticle(truthIter->second->get_trkid())
0219
0220       PHG4Particle *truthPar = truthIter->second;
0221       int previousparentid = -999;
0222       if(truthPar -> get_pid() != 22) continue;//end with a photon
0223
0224
0225       int parentid = truthPar->get_parent_id();
0226       previousparentid = truthPar->get_track_id();
0227       while (parentid != 0)
0228     {
0229       previousparentid = parentid;
0230       PHG4Particle* mommypart = truthinfo->GetParticle(parentid);
0231       parentid = mommypart->get_parent_id();
0232     }
0233
0234       if(previousparentid > 0)truthPar = truthinfo -> GetParticle(previousparentid);
0235       if(truthPar -> get_pid() != 22) continue;//start with a photon
0236
0237       for(HepMC::GenEvent::particle_const_iterator p = theEvent -> particles_begin(); p != theEvent -> particles_end(); ++p)
0238     {
0239       assert(*p);
0240
0241       if((*p) -> barcode() == truthPar -> get_barcode())
0242         {
0243           for(HepMC::GenVertex::particle_iterator mother = (*p)->production_vertex()->particles_begin(HepMC::parents);
0244           mother != (*p)->production_vertex()->particles_end(HepMC::parents); ++mother)
0245         {
0246           //std::cout << "Mother pid: " << (*mother) -> pdg_id() << std::endl;
0247           //if((*mother) -> pdg_id() != 21 || (*mother) -> pdg_id() != 1 || (*mother) -> pdg_id() != 2) continue;
0248           if((*mother) -> pdg_id() != 22) continue;
0249
0250           HepMC::FourVector moMomentum = (*mother) -> momentum();
0251
0252           m_photonE.push_back(moMomentum.e());
0253           m_photonEta.push_back(moMomentum.pseudoRapidity());
0254           m_photonPhi.push_back(moMomentum.phi());
0255           m_photonPt.push_back(sqrt(moMomentum.px()*moMomentum.px()+moMomentum.py()*moMomentum.py()));
0256
0257
0258         }
0259
0260         }
0261     }
0262     }
0263
0264
0265   T -> Fill();
0266
0267 return Fun4AllReturnCodes::EVENT_OK;
0268 }
0269
0270 //____________________________________________________________________________..
0271 int isoCluster::ResetEvent(PHCompositeNode *topNode)
0272 {
0273
0274   m_clusterPhi.clear();
0275   m_clusterEta.clear();
0276   m_clusterE.clear();
0277   m_clusterPt.clear();
0278   m_clusterEtIso.clear();
0279   m_clusterChisq.clear();
0280   m_photonPhi.clear();
0281   m_photonEta.clear();
0282   m_photonE.clear();
0283   m_photonPt.clear();
0284
0285
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 }
0288
0289 //____________________________________________________________________________..
0290 int isoCluster::EndRun(const int runnumber)
0291 {
0292   std::cout << "isoCluster::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0293   return Fun4AllReturnCodes::EVENT_OK;
0294 }
0295
0296 //____________________________________________________________________________..
0297 int isoCluster::End(PHCompositeNode *topNode)
0298 {
0299   std::cout << "isoCluster::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0300   fout -> cd();
0301   T -> Write();
0302   fout -> Close();
0303   return Fun4AllReturnCodes::EVENT_OK;
0304 }
0305
0306 //____________________________________________________________________________..
0307 int isoCluster::Reset(PHCompositeNode *topNode)
0308 {
0309   std::cout << "isoCluster::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0310   return Fun4AllReturnCodes::EVENT_OK;
0311 }
0312
0313 //____________________________________________________________________________..
0314 void isoCluster::Print(const std::string &what) const
0315 {
0316   std::cout << "isoCluster::Print(const std::string &what) const Printing info for " << what << std::endl;
0317 }