hbcup/src/hijbkg_upc.cc

0001 #include "hijbkg_upc.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/getClass.h>
0007
0008 #include <g4main/PHG4TruthInfoContainer.h>
0009 #include <g4main/PHG4Particle.h>
0010 #include <ffaobjects/EventHeader.h>
0011
0012 //#include <phhepmc/PHHepMCGenEvent.h>
0013 //#include <phhepmc/PHHepMCGenEventMap.h>
0014 #pragma GCC diagnostic push
0015 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0016 //#include <HepMC/GenEvent.h>
0017 //#include <HepMC/GenParticle.h>
0018 //#include <HepMC/GenVertex.h>
0019 //#include <HepMC/IteratorRange.h>
0020 //#include <HepMC/SimpleVector.h>
0021 //#include <HepMC/GenParticle.h>
0022 #pragma GCC diagnostic pop
0023
0024 #include <TTree.h>
0025 #include <TFile.h>
0026
0027 //____________________________________________________________________________..
0028 hijbkg_upc::hijbkg_upc(const std::string &name):
0029   SubsysReco(name),
0030   Outfile(name)
0031 {
0032   std::cout << "hijbkg_upc::hijbkg_upc(const std::string &name) Calling ctor" << std::endl;
0033 }
0034 //____________________________________________________________________________..
0035 hijbkg_upc::~hijbkg_upc()
0036 {
0037   std::cout << "hijbkg_upc::~hijbkg_upc() Calling dtor" << std::endl;
0038 }
0039
0040 //____________________________________________________________________________..
0041 int hijbkg_upc::Init(PHCompositeNode * /*topNode*/)
0042 {
0043   std::cout << "hijbkg_upc::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0044
0045   out = new TFile("hijbkg.root","RECREATE");
0046
0047   T = new TTree("T","T");
0048   T -> Branch("evt",&m_evt);
0049   T -> Branch("b",&m_b);
0050   T -> Branch("cent",&m_cent);
0051   T -> Branch("part0",&m_part[0]);
0052   T -> Branch("part1",&m_part[1]);
0053   /*
0054   T -> Branch("pid",&m_pid);
0055   T -> Branch("pt",&m_pt);
0056   T -> Branch("eta",&m_eta);
0057   T -> Branch("phi",&m_phi);
0058   T -> Branch("e",&m_e);
0059   T -> Branch("psi2",&m_psi2);
0060   T -> Branch("p",&m_p);
0061   */
0062
0063
0064   return Fun4AllReturnCodes::EVENT_OK;
0065 }
0066
0067 //____________________________________________________________________________..
0068 int hijbkg_upc::InitRun(PHCompositeNode * /*topNode*/)
0069 {
0070   std::cout << "hijbkg_upc::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0071   return Fun4AllReturnCodes::EVENT_OK;
0072 }
0073
0074 //____________________________________________________________________________..
0075 int hijbkg_upc::process_event(PHCompositeNode *topNode)
0076 {
0077
0078   // Get Event Header
0079   EventHeader *evtheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0080   m_evt = evtheader->get_EvtSequence();
0081   if ( m_evt%1 == 0 ) std::cout << "Event " << m_evt << std::endl;
0082
0083
0084   //truth particle information
0085   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0086   if(!truthinfo)
0087     {
0088       std::cout << PHWHERE << "hijingTruthCheck::process_event Could not find node G4TruthInfo"  << std::endl;
0089       return Fun4AllReturnCodes::ABORTEVENT;
0090     }
0091
0092   PHG4TruthInfoContainer::Range truthRange = truthinfo -> GetPrimaryParticleRange();
0093   PHG4TruthInfoContainer::ConstIterator truthIter;
0094
0095   int ntracks = 0;  // stable, charged, and in acceptance
0096
0097   // Loop over the monte carlo particles (eg hijing) in the truth container
0098   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0099   {
0100     PHG4Particle *truthPar = truthIter->second;
0101     double px = truthPar->get_px();
0102     double py = truthPar->get_py();
0103     double pz = truthPar->get_pz();
0104     double e = truthPar->get_e();
0105
0106     // return 1 if particle is stable and charged
0107     int pid = truthPar->get_pid();
0108
0109     if ( isStableCharged( pid ) )
0110     {
0111         // Stop processing the event since there are already more than two tracks
0112         if ( ntracks > 2 ) return Fun4AllReturnCodes::EVENT_OK;
0113
0114         float eta = getEta( truthPar );
0115         float pt = getpT( truthPar );
0116         if ( fabs(eta) < 1.1 && pt > 0.3 )  // in sPHENIX acceptance
0117         {
0118             m_part[ntracks].SetPdgCode( pid );
0119             m_part[ntracks].SetMomentum( px, py, pz, e );
0120             ntracks++;
0121         }
0122
0123     }
0124     //std::cout << truthPar->get_pid() << "\t" << px << "\t" << py << "\t" << pz << std::endl;
0125
0126
0127     //if(abs(truthPar -> get_pid()) >= 12 && abs(truthPar -> get_pid()) <= 16) continue;
0128
0129     /*
0130        m_pid.push_back(truthPar -> get_pid());
0131
0132        m_p.push_back(getP(truthPar));
0133
0134        m_e.push_back(truthPar -> get_e());
0135
0136        m_eta.push_back(getEta(truthPar));
0137
0138        m_phi.push_back(getPhi(truthPar));
0139
0140        m_pt.push_back(getpT(truthPar));
0141        */
0142
0143     //if(getpT(truthPar) < 0 || getpT(truthPar) > 10)std::cout << "Found particle with pt: " << getpT(truthPar) << std::endl;
0144
0145   } // end of loop over truthcontainer
0146
0147   //std::cout << "Saw " << npart << " particles" << std::endl;
0148   if ( ntracks==2 ) T->Fill();
0149
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 //____________________________________________________________________________..
0154 int hijbkg_upc::ResetEvent(PHCompositeNode * /*topNode*/)
0155 {
0156
0157     m_pid.clear();
0158     m_e.clear();
0159     m_eta.clear();
0160     m_pt.clear();
0161     m_phi.clear();
0162
0163     return Fun4AllReturnCodes::EVENT_OK;
0164 }
0165
0166 //____________________________________________________________________________..
0167 int hijbkg_upc::EndRun(const int runnumber)
0168 {
0169     std::cout << "hijbkg_upc::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0170     return Fun4AllReturnCodes::EVENT_OK;
0171 }
0172
0173 //____________________________________________________________________________..
0174 int hijbkg_upc::End(PHCompositeNode * /*topNode*/)
0175 {
0176     out -> cd();
0177     T -> Write();
0178     out -> Close();
0179
0180     std::cout << "hijbkg_upc::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0181     return Fun4AllReturnCodes::EVENT_OK;
0182 }
0183
0184 //____________________________________________________________________________..
0185 int hijbkg_upc::Reset(PHCompositeNode * /*topNode*/)
0186 {
0187     std::cout << "hijbkg_upc::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0188     return Fun4AllReturnCodes::EVENT_OK;
0189 }
0190
0191 //____________________________________________________________________________..
0192 void hijbkg_upc::Print(const std::string &what) const
0193 {
0194     std::cout << "hijbkg_upc::Print(const std::string &what) const Printing info for " << what << std::endl;
0195 }
0196
0197 // check if particle is stable and charged
0198 int hijbkg_upc::isStableCharged(int pid)
0199 {
0200     if ( abs(pid) == 211 ) return 1;      // pi+/pi-
0201     if ( abs(pid) == 211 ) return 1;      // pi+/pi-
0202     if ( abs(pid) == 321 ) return 1;      // K+/K-
0203     if ( abs(pid) == 2212 ) return 1;     // p+/p-
0204     if ( abs(pid) == 11 ) return 1;       // e+/e-
0205     if ( abs(pid) == 13 ) return 1;       // mu+/mu-
0206
0207
0208     return 0;
0209 }
0210
0211 //____________________________________________________________________________..
0212 float hijbkg_upc::getEta(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 0.5*log((p+pz)/(p-pz));
0220 }
0221
0222 //____________________________________________________________________________..
0223 float hijbkg_upc::getPhi(PHG4Particle *particle)
0224 {
0225     float phi = atan2(particle -> get_py(),particle -> get_px());
0226     return phi;
0227 }
0228
0229 //____________________________________________________________________________..
0230 float hijbkg_upc::getpT(PHG4Particle *particle)
0231 {
0232     float px = particle -> get_px();
0233     float py = particle -> get_py();
0234
0235     float pt = sqrt(pow(px,2) + pow(py,2));
0236     return pt;
0237 }
0238
0239 //____________________________________________________________________________..
0240 float hijbkg_upc::getP(PHG4Particle *particle)
0241 {
0242     float px = particle -> get_px();
0243     float py = particle -> get_py();
0244     float pz = particle -> get_pz();
0245     float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0246
0247     return p;
0248 }