pythiaInfo/src/jetHistogrammer.cc

0001 //____________________________________________________________________________..
0002 // This module is intended to look at truth level and pseudo jet level jets
0003 // and measure the spectrum from min bias pythia. The output will be used
0004 // to determine the point of 100% efficiency between two independent data sets
0005 // e.g. the 10 and 30GeV triggered data sets
0006 //____________________________________________________________________________..
0007
0008 #include "jetHistogrammer.h"
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/PHTFileServer.h>
0012
0013 //fast jet stuff
0014 #include <fastjet/ClusterSequence.hh>
0015 #include <fastjet/JetDefinition.hh>
0016 #include <fastjet/PseudoJet.hh>
0017
0018 //jet stuff
0019 #include <g4jets/JetMap.h>
0020
0021 #include <algorithm>
0022 #include <cmath>
0023 #include <cstdlib>
0024 #include <iostream>
0025 #include <memory>
0026 #include <utility>
0027 #include <vector>
0028
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/getClass.h>
0031
0032
0033 #include <phhepmc/PHHepMCGenEvent.h>
0034 #include <phhepmc/PHHepMCGenEventMap.h>
0035 #pragma GCC diagnostic push
0036 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0037 #include <HepMC/GenEvent.h>
0038 #include <HepMC/GenParticle.h>
0039 #include <HepMC/GenVertex.h>
0040 #include <HepMC/IteratorRange.h>
0041 #include <HepMC/SimpleVector.h>
0042 #include <HepMC/GenParticle.h>
0043 #pragma GCC diagnostic pop
0044
0045 //____________________________________________________________________________..
0046 jetHistogrammer::jetHistogrammer(const std::string &name, const std::string &fileout):
0047   SubsysReco(name),
0048   m_outputFilename(fileout)
0049 {
0050   std::cout << "jetHistogrammer::jetHistogrammer(const std::string &name) Calling ctor" << std::endl;
0051 }
0052
0053 //____________________________________________________________________________..
0054 jetHistogrammer::~jetHistogrammer()
0055 {
0056
0057   for(int i = 0; i < nEtaBins; i++)ptGJet[i] = 0;
0058   ptPJet = 0;
0059   std::cout << "jetHistogrammer::~jetHistogrammer() Calling dtor" << std::endl;
0060 }
0061
0062 //____________________________________________________________________________..
0063 int jetHistogrammer::Init(PHCompositeNode *topNode)
0064 {
0065
0066   std::cout << "jetHistogrammer::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0067
0068
0069   int nBins = 100;
0070   Float_t bins[nBins+1];
0071   for(int i = 0; i < nBins+1; i++) bins[i] = 80./nBins*i;
0072
0073   for(int i = 0; i < nEtaBins; i++) ptGJet[i] = new TH1F(Form("geantJets_%dEta",i),"jets reco'd in GEANT world",nBins,bins);
0074   //0.7,1.5,inclusive
0075
0076   ptPJet = new TH1F("pythiaJets","jets reco'd in pythia world",nBins,bins);
0077
0078
0079
0080   return Fun4AllReturnCodes::EVENT_OK;
0081 }
0082
0083 //____________________________________________________________________________..
0084 int jetHistogrammer::InitRun(PHCompositeNode *topNode)
0085 {
0086   std::cout << "jetHistogrammer::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0087   return Fun4AllReturnCodes::EVENT_OK;
0088 }
0089
0090 //____________________________________________________________________________..
0091 int jetHistogrammer::process_event(PHCompositeNode *topNode)
0092 {
0093
0094   JetMap* jetsMC = findNode::getClass<JetMap>(topNode, "AntiKt_Truth_r04");
0095   if(!jetsMC)
0096     {
0097       std::cout << "jetHistogrammer::process_event Cannot find truth node!" << std::endl;
0098       exit(-1);
0099     }
0100
0101   PHHepMCGenEventMap *genEventMap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0102   if(!genEventMap)
0103     {
0104       std::cout << "jetHistogrammer::process_event Cannot find genEventMap!" << std::endl;
0105       exit(-1);
0106     }
0107
0108   PHHepMCGenEvent *genEvent = genEventMap -> get(1);
0109   if(!genEvent)
0110     {
0111       std::cout << "jetHistogrammer::process_event Cannot find genEvent!" << std::endl;
0112       exit(-1);
0113     }
0114
0115   HepMC::GenEvent *theEvent = genEvent -> getEvent();
0116
0117
0118   float maxJetPt = 0;
0119   float maxEta = -1;
0120   for(JetMap::Iter iter = jetsMC -> begin(); iter != jetsMC -> end(); ++iter)
0121     {
0122       Jet *truthjet = iter -> second;
0123
0124       if(truthjet -> get_pt() > maxJetPt)
0125     {
0126       maxJetPt = truthjet -> get_pt();
0127       maxEta = truthjet -> get_eta();
0128     }
0129     }
0130
0131
0132   if(getEtaBin(maxEta) != -1)
0133     {
0134       ptGJet[3] -> Fill(maxJetPt);
0135       for(int i = (int)getEtaBin(maxEta);  i < nEtaBins - 1; i++) ptGJet[i] -> Fill(maxJetPt);
0136     }
0137   else
0138     {
0139       ptGJet[3] -> Fill(maxJetPt);
0140     }
0141
0142   std::vector<fastjet::PseudoJet> pseudojets;
0143   unsigned int i = 0;
0144   for(HepMC::GenEvent::particle_const_iterator p = theEvent -> particles_begin(); p != theEvent -> particles_end(); ++p)
0145     {
0146       if(!(*p) -> status()) continue; //only stable particles
0147
0148       // remove some particles (muons, taus, neutrinos)...
0149       // 12 == nu_e
0150       // 13 == muons
0151       // 14 == nu_mu
0152       // 15 == taus
0153       // 16 == nu_tau
0154       if (abs((*p) -> pdg_id() >= 12) && abs((*p) -> pdg_id() <= 16))  continue;
0155
0156       // remove acceptance... _etamin,_etamax
0157       if (((*p) -> momentum().px() == 0.0) && ((*p) -> momentum().py() == 0.0)) continue;  // avoid pt=0
0158
0159       if (abs((*p) -> momentum().pseudoRapidity()) > 1.5) continue;
0160
0161       fastjet::PseudoJet pseudojet((*p)->momentum().px(),
0162                                    (*p)->momentum().py(),
0163                                    (*p)->momentum().pz(),
0164                                    (*p)->momentum().e());
0165       pseudojet.set_user_index(i);
0166       pseudojets.push_back(pseudojet);
0167       i++;
0168     }
0169   //call fastjet
0170   fastjet::JetDefinition *jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4, fastjet::E_scheme, fastjet::Best);
0171   fastjet::ClusterSequence jetFinder(pseudojets, *jetDef);
0172   std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
0173   delete jetDef;
0174
0175   maxJetPt = 0;
0176   for(unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
0177     {
0178       const double pt = sqrt(fastjets[ijet].px() * fastjets[ijet].px() + fastjets[ijet].py() * fastjets[ijet].py());
0179       if(pt > maxJetPt) maxJetPt = pt;
0180
0181     }
0182   ptPJet -> Fill(maxJetPt);
0183
0184
0185   return Fun4AllReturnCodes::EVENT_OK;
0186 }
0187
0188 //____________________________________________________________________________..
0189 int jetHistogrammer::ResetEvent(PHCompositeNode *topNode)
0190 {
0191   return Fun4AllReturnCodes::EVENT_OK;
0192 }
0193
0194 //____________________________________________________________________________..
0195 int jetHistogrammer::EndRun(const int runnumber)
0196 {
0197   std::cout << "jetHistogrammer::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0198   return Fun4AllReturnCodes::EVENT_OK;
0199 }
0200
0201 //____________________________________________________________________________..
0202 int jetHistogrammer::End(PHCompositeNode *topNode)
0203 {
0204
0205   std::cout << "jetHistogrammer opening output file: " << m_outputFilename << std::endl;
0206   PHTFileServer::get().open(m_outputFilename,"RECREATE");
0207   for(int i = 0; i < nEtaBins; i++)ptGJet[i] -> Write();
0208   ptPJet -> Write();
0209
0210   std::cout << "jetHistogrammer::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0211   return Fun4AllReturnCodes::EVENT_OK;
0212 }
0213
0214 //____________________________________________________________________________..
0215 int jetHistogrammer::Reset(PHCompositeNode *topNode)
0216 {
0217   std::cout << "jetHistogrammer::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0218   return Fun4AllReturnCodes::EVENT_OK;
0219 }
0220
0221 //____________________________________________________________________________..
0222 void jetHistogrammer::Print(const std::string &what) const
0223 {
0224   std::cout << "jetHistogrammer::Print(const std::string &what) const Printing info for " << what << std::endl;
0225 }
0226 //____________________________________________________________________________..
0227 int jetHistogrammer::getEtaBin(float eta)
0228 {
0229   int etaBin = -1;
0230   for(int i = 0; i < nEtaBins-1; i++)
0231     {
0232       if(sqrt(pow(eta,2)) >= etaBins[i] && sqrt(pow(eta,2)) < etaBins[i+1]) etaBin = i;
0233     }
0234   return etaBin;
0235
0236 }