PatExamples/plugins/PatJetAnalyzer.cc

0001 #include <map>
0002 #include <string>
0003 #include <iomanip>
0004 #include <sstream>
0005 #include <iostream>
0006
0007 #include "TH1.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015
0016 #include "DataFormats/PatCandidates/interface/Jet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 /// maximal number of bins used for the jet
0020 /// response plots
0021 static const unsigned int MAXBIN = 8;
0022 /// binning used for the jet response plots
0023 /// (NOTE BINS must have a length of MAXBIN
0024 /// +1)
0025 static const float BINS[] = {30., 40., 50., 60., 70., 80., 100., 125., 150.};
0026
0027 /**
0028    \class   PatJetAnalyzer PatJetAnalyzer.h "PhysicsTools/PatAlgos/plugins/PatJetAnalyzer.h"
0029
0030    \brief   Module to analyze pat::Jets in the context of a more complex exercise.
0031
0032    Basic quantities of jets like the transverse momentum, eta and phi as well as the
0033    invariant dijet mass are plotted. Basic histograms for a jet energy response plot
0034    as a function of the pt of the reference object are filled. As reference matched
0035    partons are chosen. Input parameters are:
0036
0037     - src       --> input for the pat jet collection (edm::InputTag).
0038
0039     - corrLevel --> string for the pat jet correction level.
0040 */
0041
0042 class PatJetAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0043 public:
0044   /// default contructor
0045   explicit PatJetAnalyzer(const edm::ParameterSet& cfg);
0046   /// default destructor
0047   ~PatJetAnalyzer() override {}
0048
0049 private:
0050   /// everything that needs to be done during the even loop
0051   void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
0052
0053   /// check if histogram was booked
0054   bool booked(const std::string histName) const { return hists_.find(histName) != hists_.end(); };
0055   /// fill histogram if it had been booked before
0056   void fill(const std::string histName, double value) const {
0057     if (booked(histName))
0058       hists_.find(histName)->second->Fill(value);
0059   };
0060   // print jet pt for each level of energy corrections
0061   void print(edm::View<pat::Jet>::const_iterator& jet, unsigned int idx);
0062
0063 private:
0064   /// correction level for pat jet
0065   std::string corrLevel_;
0066   /// pat jets
0067   edm::EDGetTokenT<edm::View<pat::Jet> > jetsToken_;
0068   /// management of 1d histograms
0069   std::map<std::string, TH1F*> hists_;
0070 };
0071
0072 PatJetAnalyzer::PatJetAnalyzer(const edm::ParameterSet& cfg)
0073     : corrLevel_(cfg.getParameter<std::string>("corrLevel")),
0074       jetsToken_(consumes<edm::View<pat::Jet> >(cfg.getParameter<edm::InputTag>("src"))) {
0075   usesResource(TFileService::kSharedResource);
0076   // register TFileService
0077   edm::Service<TFileService> fs;
0078
0079   // jet multiplicity
0080   hists_["mult"] = fs->make<TH1F>("mult", "N_{Jet}", 15, 0., 15.);
0081   // jet pt (for all jets)
0082   hists_["pt"] = fs->make<TH1F>("pt", "p_{T}(Jet) [GeV]", 60, 0., 300.);
0083   // jet eta (for all jets)
0084   hists_["eta"] = fs->make<TH1F>("eta", "#eta (Jet)", 60, -3., 3.);
0085   // jet phi (for all jets)
0086   hists_["phi"] = fs->make<TH1F>("phi", "#phi (Jet)", 60, 3.2, 3.2);
0087   // dijet mass (if available)
0088   hists_["mass"] = fs->make<TH1F>("mass", "M_{jj} [GeV]", 50, 0., 500.);
0089   // basic histograms for jet energy response
0090   for (unsigned int idx = 0; idx < MAXBIN; ++idx) {
0091     char buffer[10];
0092     sprintf(buffer, "jes_%i", idx);
0093     char title[50];
0094     sprintf(title, "p_{T}^{rec}/p_{T}^{gen} [%i GeV - %i GeV]", (int)BINS[idx], (int)BINS[idx + 1]);
0095     hists_[buffer] = fs->make<TH1F>(buffer, title, 100, 0., 2.);
0096   }
0097 }
0098
0099 void PatJetAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0100   // recieve jet collection label
0101   edm::Handle<edm::View<pat::Jet> > jets;
0102   event.getByToken(jetsToken_, jets);
0103
0104   // loop jets
0105   for (edm::View<pat::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0106     // print jec factors
0107     // print(jet, jet-jets->begin());
0108
0109     // fill basic kinematics
0110     fill("pt", jet->correctedJet(corrLevel_).pt());
0111     fill("eta", jet->eta());
0112     fill("phi", jet->phi());
0113     // basic plots for jet responds plot as a function of pt
0114     if (jet->genJet()) {
0115       double resp = jet->correctedJet(corrLevel_).pt() / jet->genJet()->pt();
0116       for (unsigned int idx = 0; idx < MAXBIN; ++idx) {
0117         if (BINS[idx] <= jet->genJet()->pt() && jet->genJet()->pt() < BINS[idx + 1]) {
0118           char buffer[10];
0119           sprintf(buffer, "jes_%i", idx);
0120           fill(buffer, resp);
0121         }
0122       }
0123     }
0124   }
0125   // jet multiplicity
0126   fill("mult", jets->size());
0127   // invariant dijet mass for first two leading jets
0128   if (jets->size() > 1) {
0129     fill("mass", ((*jets)[0].p4() + (*jets)[1].p4()).mass());
0130   }
0131 }
0132
0133 void PatJetAnalyzer::print(edm::View<pat::Jet>::const_iterator& jet, unsigned int idx) {
0134   //edm::LogInfo log("JEC");
0135   std::cout << "[" << idx << "] :: eta=" << std::setw(10) << jet->eta() << " phi=" << std::setw(10) << jet->phi()
0136             << " size: " << jet->availableJECLevels().size() << std::endl;
0137   for (unsigned int idx = 0; idx < jet->availableJECLevels().size(); ++idx) {
0138     pat::Jet correctedJet;
0139     if (jet->availableJECLevels()[idx].find("L5Flavor") != std::string::npos ||
0140         jet->availableJECLevels()[idx].find("L7Parton") != std::string::npos) {
0141       correctedJet = jet->correctedJet(idx, pat::JetCorrFactors::UDS);
0142     } else {
0143       correctedJet = jet->correctedJet(idx, pat::JetCorrFactors::NONE);
0144     }
0145     std::cout << std::setw(10) << correctedJet.currentJECLevel() << " pt=" << std::setw(10) << correctedJet.pt()
0146               << std::endl;
0147   }
0148 }
0149
0150 #include "FWCore/Framework/interface/MakerMacros.h"
0151 DEFINE_FWK_MODULE(PatJetAnalyzer);