JetAnalyzers/src/JetAnaPythia.cc

0001 // Name: JetAnaPythia
0002 // Description:  Example of analysis of Pythia produced partons & jets
0003 //               Based on Kostas Kousouris' templated JetPlotsExample.
0004 //               Plots are tailored to needs of dijet mass and ratio analysis.
0005 // Author: R. Harris
0006 // Date:  28 - Oct - 2008
0007 #include "RecoJets/JetAnalyzers/interface/JetAnaPythia.h"
0008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0009 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0010 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0011 #include "DataFormats/JetReco/interface/CaloJet.h"
0012 #include "DataFormats/JetReco/interface/PFJet.h"
0013 #include "DataFormats/JetReco/interface/GenJet.h"
0014 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0015 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0016 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/Run.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include <TFile.h>
0021 #include <cmath>
0022 using namespace edm;
0023 using namespace reco;
0024 using namespace std;
0025 ////////////////////////////////////////////////////////////////////////////////////////
0026 template <class Jet>
0027 JetAnaPythia<Jet>::JetAnaPythia(edm::ParameterSet const& cfg) {
0028   JetAlgorithm = cfg.getParameter<std::string>("JetAlgorithm");
0029   HistoFileName = cfg.getParameter<std::string>("HistoFileName");
0030   NJets = cfg.getParameter<int>("NJets");
0031   debug = cfg.getParameter<bool>("debug");
0032   eventsGen = cfg.getParameter<int>("eventsGen");
0033   anaLevel = cfg.getParameter<std::string>("anaLevel");
0034   xsecGen = cfg.getParameter<vector<double> >("xsecGen");
0035   ptHatEdges = cfg.getParameter<vector<double> >("ptHatEdges");
0036 }
0037 ////////////////////////////////////////////////////////////////////////////////////////
0038 template <class Jet>
0039 void JetAnaPythia<Jet>::beginJob() {
0040   TString hname;
0041   m_file = new TFile(HistoFileName.c_str(), "RECREATE");
0042   /////////// Booking histograms //////////////////////////
0043   const int nMassBins = 103;
0044   double massBoundaries[nMassBins + 1] = {
0045       1,    3,    6,    10,    16,    23,    31,    40,    50,    61,    74,    88,    103,   119,  137,
0046       156,  176,  197,  220,   244,   270,   296,   325,   354,   386,   419,   453,   489,   526,  565,
0047       606,  649,  693,  740,   788,   838,   890,   944,   1000,  1058,  1118,  1181,  1246,  1313, 1383,
0048       1455, 1530, 1607, 1687,  1770,  1856,  1945,  2037,  2132,  2231,  2332,  2438,  2546,  2659, 2775,
0049       2895, 3019, 3147, 3279,  3416,  3558,  3704,  3854,  4010,  4171,  4337,  4509,  4686,  4869, 5058,
0050       5253, 5455, 5663, 5877,  6099,  6328,  6564,  6808,  7060,  7320,  7589,  7866,  8152,  8447, 8752,
0051       9067, 9391, 9726, 10072, 10430, 10798, 11179, 11571, 11977, 12395, 12827, 13272, 13732, 14000};
0052
0053   hname = "JetPt";
0054   m_HistNames1D[hname] = new TH1F(hname, hname, 500, 0, 5000);
0055
0056   hname = "JetEta";
0057   m_HistNames1D[hname] = new TH1F(hname, hname, 120, -6, 6);
0058
0059   hname = "JetPhi";
0060   m_HistNames1D[hname] = new TH1F(hname, hname, 100, -M_PI, M_PI);
0061
0062   hname = "NumberOfJets";
0063   m_HistNames1D[hname] = new TH1F(hname, hname, 100, 0, 100);
0064
0065   hname = "DijetMass";
0066   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0067
0068   hname = "DijetMassWt";
0069   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0070   m_HistNames1D.find(hname)->second->Sumw2();
0071
0072   hname = "DijetMassIn";
0073   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0074
0075   hname = "DijetMassInWt";
0076   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0077   m_HistNames1D.find(hname)->second->Sumw2();
0078
0079   hname = "DijetMassOut";
0080   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0081
0082   hname = "DijetMassOutWt";
0083   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0084   m_HistNames1D.find(hname)->second->Sumw2();
0085
0086   hname = "ResonanceMass";
0087   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0088
0089   hname = "DipartonMass";
0090   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0091
0092   hname = "DipartonMassWt";
0093   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0094   m_HistNames1D.find(hname)->second->Sumw2();
0095
0096   hname = "DipartonMassIn";
0097   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0098
0099   hname = "DipartonMassInWt";
0100   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0101   m_HistNames1D.find(hname)->second->Sumw2();
0102
0103   hname = "DipartonMassOut";
0104   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0105
0106   hname = "DipartonMassOutWt";
0107   m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
0108   m_HistNames1D.find(hname)->second->Sumw2();
0109
0110   hname = "PtHat";
0111   m_HistNames1D[hname] = new TH1F(hname, hname, 1000, 0, 5000);
0112
0113   hname = "PtHatWt";
0114   m_HistNames1D[hname] = new TH1F(hname, hname, 1000, 0, 5000);
0115   m_HistNames1D.find(hname)->second->Sumw2();
0116
0117   hname = "PtHatFine";
0118   m_HistNames1D[hname] = new TH1F(hname, hname, 5000, 0, 5000);
0119
0120   hname = "PtHatFineWt";
0121   m_HistNames1D[hname] = new TH1F(hname, hname, 5000, 0, 5000);
0122   m_HistNames1D.find(hname)->second->Sumw2();
0123
0124   mcTruthTree_ = new TTree("mcTruthTree", "mcTruthTree");
0125   mcTruthTree_->Branch("xsec", &xsec, "xsec/F");
0126   mcTruthTree_->Branch("weight", &weight, "weight/F");
0127   mcTruthTree_->Branch("pt_hat", &pt_hat, "pt_hat/F");
0128   mcTruthTree_->Branch("nJets", &nJets, "nJets/I");
0129   mcTruthTree_->Branch("etaJet1", &etaJet1, "etaJet1/F");
0130   mcTruthTree_->Branch("etaJet2", &etaJet2, "etaJet2/F");
0131   mcTruthTree_->Branch("ptJet1", &ptJet1, "ptJet1/F");
0132   mcTruthTree_->Branch("ptJet2", &ptJet2, "ptJet2/F");
0133   mcTruthTree_->Branch("diJetMass", &diJetMass, "diJetMass/F");
0134   mcTruthTree_->Branch("etaPart1", &etaPart1, "etaPart1/F");
0135   mcTruthTree_->Branch("etaPart2", &etaPart2, "etaPart2/F");
0136   mcTruthTree_->Branch("ptPart1", &ptPart1, "ptPart1/F");
0137   mcTruthTree_->Branch("ptPart2", &ptPart2, "ptPart2/F");
0138   mcTruthTree_->Branch("diPartMass", &diPartMass, "diPartMass/F");
0139 }
0140 ////////////////////////////////////////////////////////////////////////////////////////
0141 template <class Jet>
0142 void JetAnaPythia<Jet>::analyze(edm::Event const& evt, edm::EventSetup const& iSetup) {
0143   int notDone = 1;
0144   while (notDone) {  //while loop to allow us to tailor the analysis level for faster running.
0145     TString hname;
0146
0147     // Process Info
0148
0149     //edm::Handle< double > genEventScale;
0150     //evt.getByLabel("genEventScale", genEventScale );
0151     //pt_hat = *genEventScale;
0152
0153     edm::Handle<edm::HepMCProduct> MCevt;
0154     evt.getByLabel("generatorSmeared", MCevt);
0155     HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(MCevt->GetEvent()));
0156
0157     double pthat = myGenEvent->event_scale();
0158     pt_hat = float(pthat);
0159
0160     delete myGenEvent;
0161
0162     if (anaLevel != "generating") {  //We are not generating events, so xsec is there
0163                                      //edm::Handle< GenRunInfoProduct > genInfoProduct;
0164                                      ///evt.getRun().getByLabel("generator", genInfoProduct );
0165                                      //xsec = (double)genInfoProduct->externalXSecLO();
0166       xsec = 0.0;
0167       if (ptHatEdges.size() > xsecGen.size()) {
0168         for (unsigned int i_pthat = 0; i_pthat < xsecGen.size(); ++i_pthat) {
0169           if (pthat >= ptHatEdges[i_pthat] && pthat < ptHatEdges[i_pthat + 1])
0170             xsec = float(xsecGen[i_pthat]);
0171         }
0172       } else {
0173         std::cout << "Number of PtHat bin edges too small. Xsec set to zero" << std::endl;
0174       }
0175     } else {
0176       xsec = xsecGen[0];  //Generating events, no xsec in event, get xsec from user input
0177     }
0178     if (debug)
0179       std::cout << "cross section=" << xsec << " pb" << std::endl;
0180     weight = xsec / eventsGen;
0181
0182     if (debug)
0183       std::cout << "pt_hat=" << pt_hat << std::endl;
0184     hname = "PtHat";
0185     FillHist1D(hname, pt_hat, 1.0);
0186     hname = "PtHatFine";
0187     FillHist1D(hname, pt_hat, 1.0);
0188     hname = "PtHatWt";
0189     FillHist1D(hname, pt_hat, weight);
0190     hname = "PtHatFineWt";
0191     FillHist1D(hname, pt_hat, weight);
0192     if (anaLevel == "PtHatOnly")
0193       break;  //ptHatOnly should be very fast
0194
0195     // Jet Info
0196     math::XYZTLorentzVector p4jet[2];
0197     float etajet[2];
0198     /////////// Get the jet collection //////////////////////
0199     Handle<JetCollection> jets;
0200     evt.getByLabel(JetAlgorithm, jets);
0201     typename JetCollection::const_iterator i_jet;
0202     int index = 0;
0203
0204     /////////// Count the jets in the event /////////////////
0205     hname = "NumberOfJets";
0206     nJets = jets->size();
0207     FillHist1D(hname, nJets, 1.0);
0208
0209     // Two Leading Jet Info
0210     for (i_jet = jets->begin(); i_jet != jets->end() && index < 2; ++i_jet) {
0211       hname = "JetPt";
0212       FillHist1D(hname, i_jet->pt(), 1.0);
0213       hname = "JetEta";
0214       FillHist1D(hname, i_jet->eta(), 1.0);
0215       hname = "JetPhi";
0216       FillHist1D(hname, i_jet->phi(), 1.0);
0217       p4jet[index] = i_jet->p4();
0218       etajet[index] = i_jet->eta();
0219       if (debug)
0220         std::cout << "jet " << index + 1 << ": pt=" << i_jet->pt() << ", eta=" << etajet[index] << std::endl;
0221       index++;
0222     }
0223
0224     // TTree variables //
0225     etaJet1 = etajet[0];
0226     etaJet2 = etajet[1];
0227     ptJet1 = p4jet[0].pt();
0228     ptJet2 = p4jet[1].pt();
0229     diJetMass = (p4jet[0] + p4jet[1]).mass();
0230
0231     ///  Histograms for Dijet Mass Analysis  ////
0232     if (index == 2 && abs(etaJet1) < 1.3 && abs(etaJet2) < 1.3) {
0233       hname = "DijetMass";
0234       FillHist1D(hname, diJetMass, 1.0);
0235       hname = "DijetMassWt";
0236       FillHist1D(hname, diJetMass, weight);
0237     }
0238
0239     /// Histograms for Dijet Ratio Analysis: Inner region ///
0240     if (index == 2 && abs(etaJet1) < 0.7 && abs(etaJet2) < 0.7) {
0241       hname = "DijetMassIn";
0242       FillHist1D(hname, diJetMass, 1.0);
0243       hname = "DijetMassInWt";
0244       FillHist1D(hname, diJetMass, weight);
0245     }
0246     /// Histograms for Dijet Ratio Analysis: Outer region ////
0247     if (index == 2 && (abs(etaJet1) > 0.7 && abs(etaJet1) < 1.3) && (abs(etaJet2) > 0.7 && abs(etaJet2) < 1.3)) {
0248       hname = "DijetMassOut";
0249       FillHist1D(hname, diJetMass, 1.0);
0250       hname = "DijetMassOutWt";
0251       FillHist1D(hname, diJetMass, weight);
0252     }
0253     if (anaLevel == "Jets")
0254       break;  //Jets level for samples without genParticles
0255
0256     // Parton Info
0257     edm::Handle<std::vector<reco::GenParticle> > genParticlesHandle_;
0258     evt.getByLabel("genParticles", genParticlesHandle_);
0259     if (debug)
0260       for (size_t i = 0; i < genParticlesHandle_->size(); ++i) {
0261         const reco::GenParticle& p = (*genParticlesHandle_)[i];
0262         int id = p.pdgId();
0263         int st = p.status();
0264         const math::XYZTLorentzVector& genP4 = p.p4();
0265         if (i >= 2 && i <= 8)
0266           std::cout << "particle " << i << ": id=" << id << ", status=" << st << ", mass=" << genP4.mass()
0267                     << ", pt=" << genP4.pt() << ", eta=" << genP4.eta() << std::endl;
0268       }
0269     // Examine the 7th particle in pythia.
0270     // It should be either a resonance (abs(id)>=32) or the first outgoing parton
0271     // for the processes we will consider: dijet resonances, QCD, or QCD +contact interactions.
0272     const reco::GenParticle& p = (*genParticlesHandle_)[6];
0273     int id = p.pdgId();
0274     math::XYZTLorentzVector resonance_p, parton1_p, parton2_p;
0275     if (abs(id) >= 32) {
0276       /// We are looking at dijet resonances. ////
0277       resonance_p = p.p4();
0278       hname = "ResonanceMass";
0279       FillHist1D(hname, resonance_p.mass(), 1.0);
0280       const reco::GenParticle& q = (*genParticlesHandle_)[7];
0281       parton1_p = q.p4();
0282       const reco::GenParticle& r = (*genParticlesHandle_)[8];
0283       parton2_p = r.p4();
0284       if (debug)
0285         std::cout << "Resonance mass=" << resonance_p.mass() << ", parton 1 pt=" << parton1_p.pt()
0286                   << ", parton 2 pt=" << parton2_p.pt() << ", diparton mass=" << (parton1_p + parton2_p).mass()
0287                   << std::endl;
0288     } else {
0289       ///  We are looking at QCD   ////
0290       parton1_p = p.p4();
0291       const reco::GenParticle& q = (*genParticlesHandle_)[7];
0292       parton2_p = q.p4();
0293       if (debug)
0294         std::cout << "parton 1 pt=" << parton1_p.pt() << ", parton 2 pt=" << parton2_p.pt()
0295                   << ", diparton mass=" << (parton1_p + parton2_p).mass() << std::endl;
0296     }
0297
0298     etaPart1 = parton1_p.eta();
0299     etaPart2 = parton2_p.eta();
0300     ptPart1 = parton1_p.pt();
0301     ptPart2 = parton2_p.pt();
0302     diPartMass = (parton1_p + parton2_p).mass();
0303     /// Diparton mass for dijet mass analysis  ////
0304     if (abs(etaPart1) < 1.3 && abs(etaPart2) < 1.3) {
0305       hname = "DipartonMass";
0306       FillHist1D(hname, diPartMass, 1.0);
0307       hname = "DipartonMassWt";
0308       FillHist1D(hname, diPartMass, weight);
0309     }
0310     /// Diparton mass for dijet ratio analysis: inner region ///
0311     if (abs(etaPart1) < 0.7 && abs(etaPart2) < 0.7) {
0312       hname = "DipartonMassIn";
0313       FillHist1D(hname, diPartMass, 1.0);
0314       hname = "DipartonMassInWt";
0315       FillHist1D(hname, diPartMass, weight);
0316     }
0317     /// Diparton mass for dijet ratio analysis: outer region ///
0318     if ((abs(etaPart1) > 0.7 && abs(etaPart1) < 1.3) && (abs(etaPart2) > 0.7 && abs(etaPart2) < 1.3)) {
0319       hname = "DipartonMassOut";
0320       FillHist1D(hname, diPartMass, 1.0);
0321       hname = "DipartonMassOutWt";
0322       FillHist1D(hname, diPartMass, weight);
0323     }
0324
0325     // Fill the TTree //
0326     mcTruthTree_->Fill();
0327
0328     notDone = 0;  //We are done, exit the while loop
0329   }               //end of while
0330 }
0331 ////////////////////////////////////////////////////////////////////////////////////////
0332 template <class Jet>
0333 void JetAnaPythia<Jet>::endJob() {
0334   /////////// Write Histograms in output ROOT file ////////
0335   if (m_file != nullptr) {
0336     m_file->cd();
0337     mcTruthTree_->Write();
0338     for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
0339       hid->second->Write();
0340     delete m_file;
0341     m_file = nullptr;
0342   }
0343 }
0344 ////////////////////////////////////////////////////////////////////////////////////////
0345 template <class Jet>
0346 void JetAnaPythia<Jet>::FillHist1D(const TString& histName, const Double_t& value, const Double_t& wt) {
0347   std::map<TString, TH1*>::iterator hid = m_HistNames1D.find(histName);
0348   if (hid == m_HistNames1D.end())
0349     std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
0350   else
0351     hid->second->Fill(value, wt);
0352 }
0353 /////////// Register Modules ////////
0354 #include "FWCore/Framework/interface/MakerMacros.h"
0355 /////////// Calo Jet Instance ////////
0356 typedef JetAnaPythia<CaloJet> CaloJetAnaPythia;
0357 DEFINE_FWK_MODULE(CaloJetAnaPythia);
0358 /////////// Cen Jet Instance ////////
0359 typedef JetAnaPythia<GenJet> GenJetAnaPythia;
0360 DEFINE_FWK_MODULE(GenJetAnaPythia);
0361 /////////// PF Jet Instance ////////
0362 typedef JetAnaPythia<PFJet> PFJetAnaPythia;
0363 DEFINE_FWK_MODULE(PFJetAnaPythia);