analysis/EICAnalysis/Quarkonia2LeptonsMC.C

0001 #include "Quarkonia2LeptonsMC.h"
0002
0003 #include <phool/getClass.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008
0009 //#include <g4main/PHG4TruthInfoContainer.h>
0010 //#include <g4main/PHG4Particle.h>
0011
0012 #include <TLorentzVector.h>
0013 #include <TNtuple.h>
0014 #include <TFile.h>
0015 #include <TString.h>
0016 #include <TH2D.h>
0017 #include <TDatabasePDG.h>
0018
0019 //#include <cmath>
0020 //#include <iostream>
0021
0022 #include <phhepmc/PHHepMCGenEvent.h>
0023 #include <phhepmc/PHHepMCGenEventMap.h>
0024 #include <HepMC/GenEvent.h>
0025 #include <HepMC/GenVertex.h>
0026
0027 using namespace std;
0028
0029 Quarkonia2LeptonsMC::Quarkonia2LeptonsMC(std::string filename) :
0030   SubsysReco("Quarkonia2LeptonsMC" ),
0031   _embedding_id(1)
0032 {
0033   _foutname = filename;
0034   _tree_quarkonia = NULL;
0035 }
0036
0037 int
0038 Quarkonia2LeptonsMC::Init(PHCompositeNode *topNode)
0039 {
0040   _verbose = false;
0041   _ievent = 0;
0042
0043   _fout_root = new TFile(_foutname.c_str(), "RECREATE");
0044   _tree_quarkonia = new TNtuple("tree_quarkonia","quarkonia info","ev:px:py:pz:e:pt:eta:phi:invmass:ndaughter:d1d2angle");
0045
0046   return 0;
0047 }
0048
0049
0050 int
0051 Quarkonia2LeptonsMC::process_event(PHCompositeNode *topNode)
0052 {
0053
0054   _ievent ++;
0055
0056   //  cout << "Process event # " << _ievent << endl;
0057
0058   PHHepMCGenEventMap * geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0059   if (!geneventmap)
0060   {
0061     std::cout <<PHWHERE<<" - Fatal error - missing node PHHepMCGenEventMap"<<std::endl;
0062     return Fun4AllReturnCodes::ABORTRUN;
0063   }
0064
0065   PHHepMCGenEvent *genevt = geneventmap->get(_embedding_id);
0066   if (!genevt)
0067   {
0068     std::cout <<PHWHERE<<" - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "<<_embedding_id;
0069     std::cout <<". Print PHHepMCGenEventMap:";
0070     geneventmap->identify();
0071     return Fun4AllReturnCodes::ABORTRUN;
0072   }
0073
0074   HepMC::GenEvent* theEvent = genevt->getEvent();
0075
0076   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0077        p != theEvent->particles_end(); ++p)
0078     {
0079
0080       TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*p)->pdg_id() );
0081
0082       //      if ( (*p)->pdg_id() == 443 )
0083       if (TString(pdg_p->GetName()) == "J/psi" && (*p)->end_vertex() && (*p)->status() == 2 )
0084     {
0085
0086       /* quarkonia data */
0087       float px = (*p)->momentum().px();
0088       float py = (*p)->momentum().py();
0089       float pz = (*p)->momentum().pz();
0090       float pe = (*p)->momentum().e();
0091       float pt = sqrt(pow((*p)->momentum().px(),2) + pow((*p)->momentum().py(),2));
0092       float eta = (*p)->momentum().eta();
0093       float phi = (*p)->momentum().phi();
0094       float invmass = (*p)->momentum().m();
0095
0096       /* daughter = decay products */
0097       float ndaughter = (*p)->end_vertex()->particles_out_size();
0098
0099       /* calculate angle between first two daugher particles */
0100       float d1d2angle = 0;
0101       if ( ndaughter == 2 )
0102         {
0103           HepMC::GenVertex::particles_out_const_iterator p_d1 = (*p)->end_vertex()->particles_out_const_begin();
0104           HepMC::GenVertex::particles_out_const_iterator p_d2 = p_d1 + 1;
0105
0106           TLorentzVector v_d1((*p_d1)->momentum().px(),(*p_d1)->momentum().py(),(*p_d1)->momentum().pz(),(*p_d1)->momentum().e());
0107           TLorentzVector v_d2((*p_d2)->momentum().px(),(*p_d2)->momentum().py(),(*p_d2)->momentum().pz(),(*p_d2)->momentum().e());
0108
0109           d1d2angle = v_d1.Angle(v_d2.Vect()); // get angle between v_d1 and v_d2
0110         }
0111
0112       float quarkonia_data[70] = {(float)_ievent,px,py,pz,pe,pt,eta,phi,invmass,ndaughter,d1d2angle};
0113       _tree_quarkonia -> Fill(quarkonia_data);
0114
0115       //      cout << "Found a J/Psi in event " << _ievent << " with momentum " << pmom << " and which decays into " << ndecay << " particles:" << endl;
0116
0117       //HepMC::IteratorRange irange = HepMC::children;
0118 //    for ( HepMC::GenVertex::particles_out_const_iterator p2 = (*p)->end_vertex()->particles_out_const_begin();
0119 //      p2 != (*p)->end_vertex()->particles_out_const_end(); ++p2 ) {
0120 //      cout << "    -> " << (*p2)->pdg_id() << endl;
0121 //    }
0122
0123     }
0124
0125
0126 //      TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*p)->pdg_id() );
0127 //
0128 //      if (pdg_p)
0129 //        {
0130 //    if (TString(pdg_p->GetName()) == "c")
0131 //      {
0132 //        cout << "Found a c-quark!" << endl;
0133 //      }
0134 //    if (TString(pdg_p->GetName()) == "c_bar")
0135 //      {
0136 //        cout << "Found a c-antiquark!" << endl;
0137 //      }
0138 //    if (TString(pdg_p->GetName()) == "J/psi" || TString(pdg_p->GetName()) == "J/psi_di")
0139 //      {
0140 //        cout << "Found a J/Psi!" << endl;
0141 //      }
0142 //    if (TString(pdg_p->GetName()) == "Upsilon")
0143 //      {
0144 //        cout << "Found an Upsilon!" << endl;
0145 //      }
0146       //      if (TString(pdg_p->GetName()) == "J/psi" && !(*p)->end_vertex() && (*p)->status() == 1 )
0147           //            {
0148           //              if((*p)->momentum().perp()>5.) _h2->Fill((*p)->momentum().perp(), (*p)->momentum().pseudoRapidity());
0149           //              if((*p)->momentum().perp()>5.) _h1->Fill((*p)->momentum().perp());
0150           //              if((*p)->momentum().perp()>5.) {
0151           //                const float px = (*p)->momentum().x(),
0152           //                  py = (*p)->momentum().y(),
0153           //                  pz = (*p)->momentum().z(),
0154           //                  pt = (*p)->momentum().perp(),
0155           //                  eta = (*p)->momentum().pseudoRapidity(),
0156           //                  phi = (*p)->momentum().phi();
0157           //                float shower_data[70] = {(float)_ievent,px,py,pz,pt,eta,phi};
0158           //
0159           //                _ntp_gamma -> Fill(shower_data);
0160           //
0161           //                return 0;
0162           //              }
0163           //            }
0164       //        }
0165     }
0166
0167   return Fun4AllReturnCodes::ABORTEVENT;
0168 }
0169
0170 int
0171 Quarkonia2LeptonsMC::End(PHCompositeNode *topNode)
0172 {
0173   _fout_root->cd();
0174   _tree_quarkonia->Write();
0175   _fout_root->Close();
0176
0177   return 0;
0178 }
0179
0180