Pythia8Interface/test/ZJetsAnalyzer.cc

0001 #include <iostream>
0002 #include <fstream>
0003
0004 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0005 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0007
0008 // essentials !!!
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/Run.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0018
0019 #include "TH1.h"
0020
0021 #include <fastjet/JetDefinition.hh>
0022 #include <fastjet/PseudoJet.hh>
0023 #include <fastjet/ClusterSequence.hh>
0024
0025 #include "GeneratorInterface/Pythia8Interface/test/analyserhepmc/LeptonAnalyserHepMC.h"
0026 #include "GeneratorInterface/Pythia8Interface/test/analyserhepmc/JetInputHepMC.h"
0027
0028 struct ParticlePtGreater {
0029   double operator()(const HepMC::GenParticle* v1, const HepMC::GenParticle* v2) {
0030     return v1->momentum().perp() > v2->momentum().perp();
0031   }
0032 };
0033
0034 class ZJetsAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0035 public:
0036   //
0037   explicit ZJetsAnalyzer(const edm::ParameterSet&);
0038   virtual ~ZJetsAnalyzer() = default;  // no need to delete ROOT stuff
0039                                        // as it'll be deleted upon closing TFile
0040
0041   void analyze(const edm::Event&, const edm::EventSetup&) override;
0042   void beginJob() override;
0043   void beginRun(const edm::Run&, const edm::EventSetup&) override {}
0044   void endRun(const edm::Run&, const edm::EventSetup&) override;
0045
0046 private:
0047   const edm::EDGetTokenT<GenEventInfoProduct> tokenGenEvent_;
0048   const edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0049   const edm::EDGetTokenT<GenRunInfoProduct> tokenGenRun_;
0050
0051   LeptonAnalyserHepMC LA;
0052   JetInputHepMC JetInput;
0053   fastjet::Strategy strategy;
0054   fastjet::RecombinationScheme recombScheme;
0055   fastjet::JetDefinition* jetDef;
0056
0057   int icategories[6];
0058
0059   TH1D* fHist2muMass;
0060 };
0061
0062 ZJetsAnalyzer::ZJetsAnalyzer(const edm::ParameterSet& pset)
0063     : tokenGenEvent_(consumes<GenEventInfoProduct>(
0064           edm::InputTag(pset.getUntrackedParameter("moduleLabel", std::string("generator")), ""))),
0065       tokenHepMC_(consumes<edm::HepMCProduct>(
0066           edm::InputTag(pset.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0067       tokenGenRun_(consumes<GenRunInfoProduct, edm::InRun>(
0068           edm::InputTag(pset.getUntrackedParameter("moduleLabel", std::string("generator")), ""))),
0069       fHist2muMass(0) {
0070   usesResource(TFileService::kSharedResource);
0071   // actually, pset is NOT in use - we keep it here just for illustratory putposes
0072 }
0073
0074 void ZJetsAnalyzer::beginJob() {
0075   edm::Service<TFileService> fs;
0076   fHist2muMass = fs->make<TH1D>("Hist2muMass", "2-mu inv. mass", 100, 60., 120.);
0077
0078   double Rparam = 0.5;
0079   strategy = fastjet::Best;
0080   recombScheme = fastjet::E_scheme;
0081   jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, Rparam, recombScheme, strategy);
0082
0083   for (int ind = 0; ind < 6; ind++) {
0084     icategories[ind] = 0;
0085   }
0086
0087   return;
0088 }
0089
0090 void ZJetsAnalyzer::endRun(const edm::Run& r, const edm::EventSetup&) {
0091   std::ofstream testi("testi.dat");
0092   double val, errval;
0093
0094   const edm::Handle<GenRunInfoProduct>& genRunInfoProduct = r.getHandle(tokenGenRun_);
0095
0096   val = (double)genRunInfoProduct->crossSection();
0097   std::cout << std::endl;
0098   std::cout << "cross section = " << val << std::endl;
0099   std::cout << std::endl;
0100
0101   errval = 0.;
0102   if (icategories[0] > 0)
0103     errval = val / sqrt((double)(icategories[0]));
0104   testi << "pythia8_test1  1   " << val << " " << errval << " " << std::endl;
0105
0106   std::cout << std::endl;
0107   std::cout << " Events with at least 1 isolated lepton  :                     "
0108             << ((double)icategories[1]) / ((double)icategories[0]) << std::endl;
0109   std::cout << " Events with at least 2 isolated leptons :                     "
0110             << ((double)icategories[2]) / ((double)icategories[0]) << std::endl;
0111   std::cout << " Events with at least 2 isolated leptons and at least 1 jet  : "
0112             << ((double)icategories[3]) / ((double)icategories[0]) << std::endl;
0113   std::cout << " Events with at least 2 isolated leptons and at least 2 jets : "
0114             << ((double)icategories[4]) / ((double)icategories[0]) << std::endl;
0115   std::cout << std::endl;
0116
0117   val = ((double)icategories[4]) / ((double)icategories[0]);
0118   errval = 0.;
0119   if (icategories[4] > 0)
0120     errval = val / sqrt((double)icategories[4]);
0121   testi << "pythia8_test1  2   " << val << " " << errval << " " << std::endl;
0122 }
0123
0124 void ZJetsAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
0125   icategories[0]++;
0126
0127   // here's an example of accessing GenEventInfoProduct
0128   const edm::Handle<GenEventInfoProduct>& GenInfoHandle = e.getHandle(tokenGenEvent_);
0129
0130   double qScale = GenInfoHandle->qScale();
0131   double pthat = (GenInfoHandle->hasBinningValues() ? (GenInfoHandle->binningValues())[0] : 0.0);
0132   std::cout << " qScale = " << qScale << " pthat = " << pthat << std::endl;
0133   //
0134   // this (commented out) code below just exemplifies how to access certain info
0135   //
0136   //double evt_weight1 = GenInfoHandle->weights()[0]; // this is "stanrd Py6 evt weight;
0137   // corresponds to PYINT1/VINT(97)
0138   //double evt_weight2 = GenInfoHandle->weights()[1]; // in case you run in CSA mode or otherwise
0139   // use PYEVWT routine, this will be weight
0140   // as returned by PYEVWT, i.e. PYINT1/VINT(99)
0141   //std::cout << " evt_weight1 = " << evt_weight1 << std::endl;
0142   //std::cout << " evt_weight2 = " << evt_weight2 << std::endl;
0143   //double weight = GenInfoHandle->weight();
0144   //std::cout << " as returned by the weight() method, integrated event weight = " << weight << std::endl;
0145
0146   // here's an example of accessing particles in the event record (HepMCProduct)
0147   //
0148   const edm::Handle<edm::HepMCProduct>& EvtHandle = e.getHandle(tokenHepMC_);
0149
0150   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0151
0152   int nisolep = LA.nIsolatedLeptons(Evt);
0153
0154   //std::cout << "Number of leptons = " << nisolep << std::endl;
0155   if (nisolep > 0)
0156     icategories[1]++;
0157   if (nisolep > 1)
0158     icategories[2]++;
0159
0160   JetInputHepMC::ParticleVector jetInput = JetInput(Evt);
0161   std::sort(jetInput.begin(), jetInput.end(), ParticlePtGreater());
0162
0163   // Fastjet input
0164   std::vector<fastjet::PseudoJet> jfInput;
0165   jfInput.reserve(jetInput.size());
0166   for (JetInputHepMC::ParticleVector::const_iterator iter = jetInput.begin(); iter != jetInput.end(); ++iter) {
0167     jfInput.push_back(fastjet::PseudoJet(
0168         (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e()));
0169     jfInput.back().set_user_index(iter - jetInput.begin());
0170   }
0171
0172   // Run Fastjet algorithm
0173   std::vector<fastjet::PseudoJet> inclusiveJets, sortedJets, cleanedJets;
0174   fastjet::ClusterSequence clustSeq(jfInput, *jetDef);
0175
0176   // Extract inclusive jets sorted by pT (note minimum pT in GeV)
0177   inclusiveJets = clustSeq.inclusive_jets(20.0);
0178   sortedJets = sorted_by_pt(inclusiveJets);
0179
0180   cleanedJets = LA.removeLeptonsFromJets(sortedJets, Evt);
0181
0182   if (nisolep > 1) {
0183     if (cleanedJets.size() > 0)
0184       icategories[3]++;
0185     if (cleanedJets.size() > 1)
0186       icategories[4]++;
0187   }
0188
0189   return;
0190 }
0191
0192 typedef ZJetsAnalyzer ZJetsTest;
0193 DEFINE_FWK_MODULE(ZJetsTest);