quickHIJING/src/quickHIJING.cc

0001 #include "quickHIJING.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/getClass.h>
0007
0008 #include <centrality/CentralityInfo.h>
0009
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011 #include <g4main/PHG4Particle.h>
0012
0013 #include <phhepmc/PHHepMCGenEvent.h>
0014 #include <phhepmc/PHHepMCGenEventMap.h>
0015 #pragma GCC diagnostic push
0016 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0017 #include <HepMC/GenEvent.h>
0018 #include <HepMC/GenParticle.h>
0019 #include <HepMC/GenVertex.h>
0020 #include <HepMC/IteratorRange.h>
0021 #include <HepMC/SimpleVector.h>
0022 #include <HepMC/GenParticle.h>
0023 #pragma GCC diagnostic pop
0024
0025 //____________________________________________________________________________..
0026 quickHIJING::quickHIJING(const std::string &name):
0027   SubsysReco(name)
0028   ,Outfile(name)
0029 {
0030   std::cout << "quickHIJING::quickHIJING(const std::string &name) Calling ctor" << std::endl;
0031 }
0032 //____________________________________________________________________________..
0033 quickHIJING::~quickHIJING()
0034 {
0035   std::cout << "quickHIJING::~quickHIJING() Calling dtor" << std::endl;
0036 }
0037
0038 //____________________________________________________________________________..
0039 int quickHIJING::Init(PHCompositeNode *topNode)
0040 {
0041   std::cout << "quickHIJING::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0042
0043   out = new TFile("qhTest.root","RECREATE");
0044
0045   T = new TTree("T","T");
0046   T -> Branch("pid",&m_pid);
0047   T -> Branch("pt",&m_pt);
0048   T -> Branch("eta",&m_eta);
0049   T -> Branch("phi",&m_phi);
0050   T -> Branch("e",&m_e);
0051   T -> Branch("psi2",&m_psi2);
0052   T -> Branch("b",&m_b);
0053   T -> Branch("cent",&m_cent);
0054   T -> Branch("p",&m_p);
0055
0056
0057   return Fun4AllReturnCodes::EVENT_OK;
0058 }
0059
0060 //____________________________________________________________________________..
0061 int quickHIJING::InitRun(PHCompositeNode *topNode)
0062 {
0063   std::cout << "quickHIJING::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0064   return Fun4AllReturnCodes::EVENT_OK;
0065 }
0066
0067 //____________________________________________________________________________..
0068 int quickHIJING::process_event(PHCompositeNode *topNode)
0069 {
0070
0071   CentralityInfo *cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0072   if(!cent_node)
0073     {
0074       std::cout << "No cent node" << std::endl;
0075       return Fun4AllReturnCodes::ABORTEVENT;
0076     }
0077
0078   m_cent = cent_node -> get_centile(CentralityInfo::PROP::bimp);
0079   m_b = cent_node -> get_quantity(CentralityInfo::PROP::bimp);
0080
0081   //truth particle information
0082   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0083   if(!truthinfo)
0084     {
0085       std::cout << PHWHERE << "hijingTruthCheck::process_event Could not find node G4TruthInfo"  << std::endl;
0086       return Fun4AllReturnCodes::ABORTEVENT;
0087     }
0088
0089
0090
0091   //For pythia ancestory information
0092   PHHepMCGenEventMap *genEventMap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0093   if(!genEventMap)
0094     {
0095       std::cout << PHWHERE << "cemc::process_event Could not find PHHepMCGenEventMap"  << std::endl;
0096       return Fun4AllReturnCodes::ABORTEVENT;
0097     }
0098   //event level information of the above
0099   PHHepMCGenEvent *genEvent = genEventMap -> get(0);
0100   if(!genEvent)
0101     {
0102       std::cout << PHWHERE << "cemc::process_event Could not find PHHepMCGenEvent"  << std::endl;
0103       return Fun4AllReturnCodes::ABORTEVENT;
0104     }
0105
0106   HepMC::GenEvent *evt = genEvent -> getEvent();
0107
0108   HepMC::HeavyIon *hi = evt -> heavy_ion();
0109
0110   m_psi2 = hi -> event_plane_angle();
0111
0112
0113   PHG4TruthInfoContainer::Range truthRange = truthinfo -> GetPrimaryParticleRange();
0114   PHG4TruthInfoContainer::ConstIterator truthIter;
0115
0116   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0117     {
0118       PHG4Particle *truthPar = truthIter->second;
0119
0120       if(abs(truthPar -> get_pid()) >= 12 && abs(truthPar -> get_pid()) <= 16) continue;
0121
0122       m_pid.push_back(truthPar -> get_pid());
0123
0124       m_p.push_back(getP(truthPar));
0125
0126       m_e.push_back(truthPar -> get_e());
0127
0128       m_eta.push_back(getEta(truthPar));
0129
0130       m_phi.push_back(getPhi(truthPar));
0131
0132       m_pt.push_back(getpT(truthPar));
0133
0134       //if(getpT(truthPar) < 0 || getpT(truthPar) > 10)std::cout << "Found particle with pt: " << getpT(truthPar) << std::endl;
0135
0136     }
0137   //std::cout << "Saw " << npart << " particles" << std::endl;
0138   T -> Fill();
0139
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 //____________________________________________________________________________..
0144 int quickHIJING::ResetEvent(PHCompositeNode *topNode)
0145 {
0146
0147   m_pid.clear();
0148   m_e.clear();
0149   m_eta.clear();
0150   m_pt.clear();
0151   m_phi.clear();
0152
0153   return Fun4AllReturnCodes::EVENT_OK;
0154 }
0155
0156 //____________________________________________________________________________..
0157 int quickHIJING::EndRun(const int runnumber)
0158 {
0159   std::cout << "quickHIJING::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0160   return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 //____________________________________________________________________________..
0164 int quickHIJING::End(PHCompositeNode *topNode)
0165 {
0166   out -> cd();
0167   T -> Write();
0168   out -> Close();
0169
0170   std::cout << "quickHIJING::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 //____________________________________________________________________________..
0175 int quickHIJING::Reset(PHCompositeNode *topNode)
0176 {
0177  std::cout << "quickHIJING::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0178   return Fun4AllReturnCodes::EVENT_OK;
0179 }
0180
0181 //____________________________________________________________________________..
0182 void quickHIJING::Print(const std::string &what) const
0183 {
0184   std::cout << "quickHIJING::Print(const std::string &what) const Printing info for " << what << std::endl;
0185 }
0186 //____________________________________________________________________________..
0187 float quickHIJING::getEta(PHG4Particle *particle)
0188 {
0189   float px = particle -> get_px();
0190   float py = particle -> get_py();
0191   float pz = particle -> get_pz();
0192   float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0193
0194   return 0.5*log((p+pz)/(p-pz));
0195 }
0196 //____________________________________________________________________________..
0197 float quickHIJING::getPhi(PHG4Particle *particle)
0198 {
0199   float phi = atan2(particle -> get_py(),particle -> get_px());
0200   return phi;
0201 }
0202 //____________________________________________________________________________..
0203 float quickHIJING::getpT(PHG4Particle *particle)
0204 {
0205   float px = particle -> get_px();
0206   float py = particle -> get_py();
0207
0208   float pt = sqrt(pow(px,2) + pow(py,2));
0209   return pt;
0210 }
0211 //____________________________________________________________________________..
0212 float quickHIJING::getP(PHG4Particle *particle)
0213 {
0214   float px = particle -> get_px();
0215   float py = particle -> get_py();
0216   float pz = particle -> get_pz();
0217   float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0218
0219   return p;
0220 }