PartonShowerVeto/test/BasicGenJetTester.cc

0001 #include <iostream>
0002
0003 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0004
0005 // essentials !!!
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "TH1.h"
0015
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017
0018 class BasicGenJetTester : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0019 public:
0020   //
0021   explicit BasicGenJetTester(const edm::ParameterSet&);
0022   ~BasicGenJetTester() override = default;  // no need to delete ROOT stuff
0023                                             // as it'll be deleted upon closing TFile
0024
0025   void analyze(const edm::Event&, const edm::EventSetup&) override;
0026   void beginJob() override;
0027   void endJob() override {}
0028
0029 private:
0030   const double fQCut;
0031   const edm::EDGetTokenT<reco::GenJetCollection> ak5GenJetToken_;
0032
0033   TH1F* fNJets;
0034   TH1F* fNJetsAboveQCut;
0035   TH1D* fLeadingJetPt;
0036   TH1D* fLeadingJetEta;
0037   TH1D* fNext2LeadingJetPt;
0038   TH1D* fNext2LeadingJetEta;
0039   TH1D* fLowestJetHt;
0040   TH1D* fLowestJetEta;
0041 };
0042
0043 BasicGenJetTester::BasicGenJetTester(const edm::ParameterSet& pset)
0044     : fQCut(pset.getParameter<double>("qcut")),
0045       ak5GenJetToken_(consumes<reco::GenJetCollection>("ak5GenJets")),
0046       fNJets(0),
0047       fNJetsAboveQCut(0),
0048       fLeadingJetPt(0),
0049       fLeadingJetEta(0),
0050       fNext2LeadingJetPt(0),
0051       fNext2LeadingJetEta(0),
0052       fLowestJetHt(0),
0053       fLowestJetEta(0) {
0054   usesResource(TFileService::kSharedResource);
0055 }
0056
0057 void BasicGenJetTester::beginJob() {
0058   edm::Service<TFileService> fs;
0059
0060   fNJets = fs->make<TH1F>("NJets", "Number of Jets (total)", 50, 0., 50.);
0061   fNJetsAboveQCut = fs->make<TH1F>("NJetsAboveQCut", "Number of Jets (above qcut)", 10, 0., 10.);
0062   fLeadingJetPt = fs->make<TH1D>("LeadingJetPt", "Leading Jet Pt", 100, 0., 250.);
0063   fLeadingJetEta = fs->make<TH1D>("LeadingJetEta", "Leading Jet Eta", 100, -5.0, 5.0);
0064   fNext2LeadingJetPt = fs->make<TH1D>("Next2LeadingJetPt", "Next to Leading Jet Pt", 100, 0., 250.);
0065   fNext2LeadingJetEta = fs->make<TH1D>("Next2LeadingJetEta", "Next to Leading Jet Eta", 100, -5.0, 5.0);
0066   fLowestJetHt = fs->make<TH1D>("LowestJetHt", "Ht (Lowest Jet above qcut)", 100, 0., 250.);
0067   fLowestJetEta = fs->make<TH1D>("LowestJetEta", "Lowest Jet Eta", 100, -5.0, 5.0);
0068
0069   return;
0070 }
0071
0072 void BasicGenJetTester::analyze(const edm::Event& e, const edm::EventSetup&) {
0073   // here's an example of accessing GenJetCollection
0074   // find initial (unsmeared, unfiltered,...) HepMCProduct
0075   //
0076   const edm::Handle<reco::GenJetCollection>& ak5GenJetHandle = e.getHandle(ak5GenJetToken_);
0077   //const edm::Handle<reco::GenJetCollection>& ak7GenJetHandle = e.getHandle(ak7GenJetToken_);
0078
0079   int NGenJets5 = ak5GenJetHandle->size();
0080   // int NGenJets7 = ak7GenJetHandle->size();
0081
0082   if (NGenJets5 <= 0)
0083     return;
0084
0085   fNJets->Fill((float)NGenJets5);
0086
0087   int NGenJets5AboveQCut = 0;
0088   reco::GenJet GJet;
0089
0090   for (unsigned int idx = 0; idx < ak5GenJetHandle->size(); ++idx) {
0091     GJet = (*ak5GenJetHandle)[idx];
0092     double pt = GJet.pt();  //cout << ": pt=" << pt;
0093     if (pt < fQCut)
0094       continue;
0095     NGenJets5AboveQCut++;
0096   }
0097
0098   if (NGenJets5AboveQCut <= 0)
0099     return;
0100
0101   fNJetsAboveQCut->Fill((float)NGenJets5AboveQCut);
0102
0103   // leading jet
0104   //
0105   GJet = (*ak5GenJetHandle)[0];
0106   fLeadingJetPt->Fill(GJet.pt());
0107   fLeadingJetEta->Fill(GJet.eta());
0108
0109   if (NGenJets5AboveQCut <= 1)
0110     return;
0111
0112   // next-to-leading jet
0113   //
0114   GJet = (*ak5GenJetHandle)[1];
0115   fNext2LeadingJetPt->Fill(GJet.pt());
0116   fNext2LeadingJetEta->Fill(GJet.eta());
0117
0118   if (NGenJets5AboveQCut <= 2)
0119     return;
0120
0121   // lowest jet (above qcut)
0122   //
0123   GJet = (*ak5GenJetHandle)[NGenJets5AboveQCut - 1];
0124   fLowestJetHt->Fill(GJet.pt());
0125   fLowestJetEta->Fill(GJet.eta());
0126
0127   return;
0128 }
0129
0130 DEFINE_FWK_MODULE(BasicGenJetTester);