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