Modules/plugins/JetResolutionDemo.cc

0001 #include <memory>
0002 #include <iomanip>
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "JetMETCorrections/Modules/interface/JetResolution.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include <DataFormats/PatCandidates/interface/Jet.h>
0018
0019 #include <TH2.h>
0020
0021 //
0022 // class declaration
0023 //
0024
0025 class JetResolutionDemo : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0026 public:
0027   explicit JetResolutionDemo(const edm::ParameterSet&);
0028   ~JetResolutionDemo() override;
0029
0030 private:
0031   void analyze(const edm::Event&, const edm::EventSetup&) override;
0032   void endJob() override;
0033
0034   edm::EDGetTokenT<std::vector<pat::Jet>> m_jets_token;
0035   edm::EDGetTokenT<double> m_rho_token;
0036
0037   bool m_debug = false;
0038   bool m_use_conddb = false;
0039
0040   std::string m_payload;
0041   std::string m_resolutions_file;
0042   std::string m_scale_factors_file;
0043   JME::JetResolution::Token m_resolutions_token;
0044   JME::JetResolutionScaleFactor::Token m_scale_factors_token;
0045   edm::Service<TFileService> fs;
0046
0047   TH2* m_res_vs_eta = nullptr;
0048 };
0049 //
0050 //----------- Class Implementation ------------------------------------------
0051 //
0052 //---------------------------------------------------------------------------
0053 JetResolutionDemo::JetResolutionDemo(const edm::ParameterSet& iConfig) {
0054   m_jets_token = consumes<std::vector<pat::Jet>>(iConfig.getParameter<edm::InputTag>("jets"));
0055   m_rho_token = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0056   m_debug = iConfig.getUntrackedParameter<bool>("debug", false);
0057   m_use_conddb = iConfig.getUntrackedParameter<bool>("useCondDB", false);
0058
0059   if (m_use_conddb) {
0060     m_payload = iConfig.getParameter<std::string>("payload");
0061     m_resolutions_token = esConsumes(edm::ESInputTag("", m_payload));
0062     m_scale_factors_token = esConsumes(edm::ESInputTag("", m_payload));
0063   } else {
0064     m_resolutions_file = iConfig.getParameter<edm::FileInPath>("resolutionsFile").fullPath();
0065     m_scale_factors_file = iConfig.getParameter<edm::FileInPath>("scaleFactorsFile").fullPath();
0066   }
0067   usesResource(TFileService::kSharedResource);
0068 }
0069 //---------------------------------------------------------------------------
0070 JetResolutionDemo::~JetResolutionDemo() {}
0071
0072 //---------------------------------------------------------------------------
0073 void JetResolutionDemo::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0074   edm::Handle<std::vector<pat::Jet>> jets;
0075   iEvent.getByToken(m_jets_token, jets);
0076
0077   edm::Handle<double> rho;
0078   iEvent.getByToken(m_rho_token, rho);
0079
0080   // Access jet resolution and scale factor from the condition database
0081   // or from text files
0082   JME::JetResolution resolution;
0083   JME::JetResolutionScaleFactor res_sf;
0084
0085   // Two differents way to create a class instance
0086   if (m_use_conddb) {
0087     // First way, using the get() static method
0088     resolution = JME::JetResolution::get(iSetup, m_resolutions_token);
0089     res_sf = JME::JetResolutionScaleFactor::get(iSetup, m_scale_factors_token);
0090   } else {
0091     // Second way, using the constructor
0092     resolution = JME::JetResolution(m_resolutions_file);
0093     res_sf = JME::JetResolutionScaleFactor(m_scale_factors_file);
0094   }
0095
0096   if (m_debug) {
0097     // Dump resolution
0098     resolution.dump();
0099   }
0100
0101   if (!m_res_vs_eta) {
0102     // Advanced usage. Create the histogram by retriving the eta binning directly from the JetResolutionObject. This suppose you need that the parametrization only depends on eta.
0103
0104     // Get the list of bins of this object
0105     const std::vector<JME::Binning>& bins = resolution.getResolutionObject()->getDefinition().getBins();
0106
0107     // Check that the first bin is eta
0108     if ((!bins.empty()) && (bins[0] == JME::Binning::JetEta)) {
0109       const std::vector<JME::JetResolutionObject::Record> records = resolution.getResolutionObject()->getRecords();
0110       // Get all records from the object. Each record correspond to a different binning and different parameters
0111
0112       std::vector<float> etas;
0113       for (const auto& record : records) {
0114         if (etas.empty()) {
0115           etas.push_back(record.getBinsRange()[0].min);
0116           etas.push_back(record.getBinsRange()[0].max);
0117         } else {
0118           etas.push_back(record.getBinsRange()[0].max);
0119         }
0120       }
0121
0122       std::vector<float> res;
0123       for (std::size_t i = 0; i < 40; i++) {
0124         res.push_back(i * 0.005);
0125       }
0126
0127       m_res_vs_eta = fs->make<TH2F>("res_vs_eta", "res_vs_eta", etas.size() - 1, &etas[0], res.size() - 1, &res[0]);
0128     }
0129   }
0130
0131   for (const auto& jet : *jets) {
0132     if (m_debug) {
0133       std::cout << "New jet; pt=" << jet.pt() << "  eta=" << jet.eta() << "  phi=" << jet.phi()
0134                 << "  e=" << jet.energy() << "  rho=" << *rho << std::endl;
0135     }
0136
0137     // Get resolution for this jet
0138     // The method to use is getResolution(const JME::JetParameters& parameters) from a JetResolution object
0139
0140     // You *must* now in advance which variables are needed for getting the resolution. For the moment, only pt and eta are needed, but this will
0141     // probably change in the futur when PU dependency is added. Please keep an eye on the twiki page
0142     //     https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookJetEnergyResolution
0143     // to stay up-to-date. All currently supported parameters (ie, 'set' functions) are available here:
0144     //     https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookJetEnergyResolution#List_of_supported_parameters
0145
0146     // Three way to create a JetParameters object
0147
0148     // First, using 'set' functions
0149     JME::JetParameters parameters_1;
0150     parameters_1.setJetPt(jet.pt());
0151     parameters_1.setJetEta(jet.eta());
0152     parameters_1.setRho(*rho);
0153
0154     // You can also chain calls
0155
0156     JME::JetParameters parameters_2;
0157     parameters_2.setJetPt(jet.pt()).setJetEta(jet.eta()).setRho(*rho);
0158
0159     // Second, using the set() function
0160     JME::JetParameters parameters_3;
0161     parameters_3.set(JME::Binning::JetPt, jet.pt());
0162     parameters_3.set({JME::Binning::JetEta, jet.eta()});
0163     parameters_3.set({JME::Binning::Rho, *rho});
0164
0165     // Or
0166
0167     JME::JetParameters parameters_4;
0168     parameters_4.set(JME::Binning::JetPt, jet.pt()).set(JME::Binning::JetEta, jet.eta()).set(JME::Binning::Rho, *rho);
0169
0170     // Third, using a initializer_list
0171     JME::JetParameters parameters_5 = {
0172         {JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}, {JME::Binning::Rho, *rho}};
0173
0174     // Now, get the resolution
0175
0176     float r = resolution.getResolution(parameters_1);
0177     if (m_debug) {
0178       std::cout << "Resolution with parameters_1: " << r << std::endl;
0179       std::cout << "Resolution with parameters_2: " << resolution.getResolution(parameters_2) << std::endl;
0180       std::cout << "Resolution with parameters_3: " << resolution.getResolution(parameters_3) << std::endl;
0181       std::cout << "Resolution with parameters_4: " << resolution.getResolution(parameters_4) << std::endl;
0182       std::cout << "Resolution with parameters_5: " << resolution.getResolution(parameters_5) << std::endl;
0183
0184       // You can also use a shortcut to get the resolution
0185       float r2 = resolution.getResolution(
0186           {{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}, {JME::Binning::Rho, *rho}});
0187       std::cout << "Resolution using shortcut   : " << r2 << std::endl;
0188     }
0189
0190     m_res_vs_eta->Fill(jet.eta(), r);
0191
0192     // We do the same thing to access the scale factors
0193     float sf = res_sf.getScaleFactor({{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}});
0194
0195     // Access up and down variation of the scale factor
0196     float sf_up =
0197         res_sf.getScaleFactor({{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}}, Variation::UP);
0198     float sf_down =
0199         res_sf.getScaleFactor({{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}}, Variation::DOWN);
0200
0201     if (m_debug) {
0202       std::cout << "Scale factors (Nominal / Up / Down) : " << sf << " / " << sf_up << " / " << sf_down << std::endl;
0203     }
0204   }
0205 }
0206 //---------------------------------------------------------------------------
0207 void JetResolutionDemo::endJob() {}
0208 //---------------------------------------------------------------------------
0209 //define this as a plug-in
0210 DEFINE_FWK_MODULE(JetResolutionDemo);