EMCal-analysis/DirectPhoton/DirectPhotonPythia.C

0001 #include "DirectPhotonPythia.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 #include <g4main/PHG4Hit.h>
0012
0013 #include <TLorentzVector.h>
0014 #include <TTree.h>
0015 #include <TFile.h>
0016 #include <TString.h>
0017 #include <TH2D.h>
0018 #include <TDatabasePDG.h>
0019
0020 #include <cmath>
0021 #include <iostream>
0022
0023 #include <g4jets/JetMap.h>
0024 #include <g4jets/Jet.h>
0025
0026 #include <phhepmc/PHHepMCGenEvent.h>
0027 #include <phhepmc/PHHepMCGenEventMap.h>
0028 #include <HepMC/GenEvent.h>
0029 #include <HepMC/GenVertex.h>
0030 using namespace std;
0031
0032 DirectPhotonPythia::DirectPhotonPythia(std::string filename) :
0033     SubsysReco("DirectPhoton" ),
0034     _embedding_id(1)
0035 {
0036   _foutname = filename;
0037
0038   _pt_min = 25;
0039   _pt_max = 100;
0040
0041   _eta_min = -.6;
0042   _eta_max = +.6;
0043
0044   _rejection_action = Fun4AllReturnCodes::DISCARDEVENT;
0045 }
0046
0047 int
0048 DirectPhotonPythia::Init(PHCompositeNode *topNode)
0049 {
0050
0051   _verbose = true;
0052
0053   _ievent = 0;
0054
0055   _total_pass = 0;
0056
0057   _f = new TFile(_foutname.c_str(), "RECREATE");
0058
0059   _h1 = new TH1D("h1", "", 100, 0, 100.0);
0060   _h2 = new TH2D("h2", "", 100, 0, 100.0, 40, -2, +2);
0061   _h2_b = new TH2D("h2_b", "", 100, 0, 100.0, 40, -2, +2);
0062   _h2_c = new TH2D("h2_c", "", 100, 0, 100.0, 40, -2, +2);
0063   _h2all = new TH2D("h2all", "", 100, 0, 100.0, 40, -2, +2);
0064
0065  _ntp_gamma = new TNtuple("ntp_gamma","truth shower => best cluster","ev:px:py:pz:pt:eta:phi");
0066
0067   return 0;
0068
0069 }
0070
0071 int
0072 DirectPhotonPythia::process_event(PHCompositeNode *topNode)
0073 {
0074   _ievent ++;
0075
0076   PHHepMCGenEventMap * geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0077   if (!geneventmap)
0078   {
0079     std::cout <<PHWHERE<<" - Fatal error - missing node PHHepMCGenEventMap"<<std::endl;
0080     return Fun4AllReturnCodes::ABORTRUN;
0081   }
0082
0083   PHHepMCGenEvent *genevt = geneventmap->get(_embedding_id);
0084   if (!genevt)
0085   {
0086     std::cout <<PHWHERE<<" - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "<<_embedding_id;
0087     std::cout <<". Print PHHepMCGenEventMap:";
0088     geneventmap->identify();
0089     return Fun4AllReturnCodes::ABORTRUN;
0090   }
0091
0092   HepMC::GenEvent* theEvent = genevt->getEvent();
0093
0094
0095 for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0096       p != theEvent->particles_end(); ++p)
0097     {
0098       TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle(
0099           (*p)->pdg_id());
0100       if (pdg_p)
0101         {
0102           if (TString(pdg_p->GetName()) == "gamma" && !(*p)->end_vertex() && (*p)->status() == 1 )
0103             {
0104           if((*p)->momentum().perp()>5.) _h2->Fill((*p)->momentum().perp(), (*p)->momentum().pseudoRapidity());
0105               if((*p)->momentum().perp()>5.) _h1->Fill((*p)->momentum().perp());
0106
0107               if((*p)->momentum().perp()>5.) {
0108         const float px = (*p)->momentum().x(),
0109                             py = (*p)->momentum().y(),
0110                             pz = (*p)->momentum().z(),
0111                     pt = (*p)->momentum().perp(),
0112                            eta = (*p)->momentum().pseudoRapidity(),
0113                            phi = (*p)->momentum().phi();
0114                 float shower_data[70] = {(float)_ievent,px,py,pz,pt,eta,phi};
0115
0116                 _ntp_gamma -> Fill(shower_data);
0117
0118                  return 0;
0119              }
0120             }
0121         }
0122       }
0123
0124       return Fun4AllReturnCodes::ABORTEVENT;
0125 }
0126
0127 int
0128 DirectPhotonPythia::End(PHCompositeNode *topNode)
0129 {
0130
0131   if (verbosity >= DirectPhotonPythia::VERBOSITY_SOME)
0132     std::cout << __PRETTY_FUNCTION__ << " DVP PASSED " << _total_pass
0133         << " events" << std::endl;
0134
0135   _f->cd();
0136   _f->Write();
0137   //_f->Close();
0138
0139   return 0;
0140 }
0141
0142