File indexing completed on 2023-03-17 11:10:46
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
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
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
0081
0082 JME::JetResolution resolution;
0083 JME::JetResolutionScaleFactor res_sf;
0084
0085
0086 if (m_use_conddb) {
0087
0088 resolution = JME::JetResolution::get(iSetup, m_resolutions_token);
0089 res_sf = JME::JetResolutionScaleFactor::get(iSetup, m_scale_factors_token);
0090 } else {
0091
0092 resolution = JME::JetResolution(m_resolutions_file);
0093 res_sf = JME::JetResolutionScaleFactor(m_scale_factors_file);
0094 }
0095
0096 if (m_debug) {
0097
0098 resolution.dump();
0099 }
0100
0101 if (!m_res_vs_eta) {
0102
0103
0104
0105 const std::vector<JME::Binning>& bins = resolution.getResolutionObject()->getDefinition().getBins();
0106
0107
0108 if ((!bins.empty()) && (bins[0] == JME::Binning::JetEta)) {
0109 const std::vector<JME::JetResolutionObject::Record> records = resolution.getResolutionObject()->getRecords();
0110
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
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149 JME::JetParameters parameters_1;
0150 parameters_1.setJetPt(jet.pt());
0151 parameters_1.setJetEta(jet.eta());
0152 parameters_1.setRho(*rho);
0153
0154
0155
0156 JME::JetParameters parameters_2;
0157 parameters_2.setJetPt(jet.pt()).setJetEta(jet.eta()).setRho(*rho);
0158
0159
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
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
0171 JME::JetParameters parameters_5 = {
0172 {JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}, {JME::Binning::Rho, *rho}};
0173
0174
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
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
0193 float sf = res_sf.getScaleFactor({{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}});
0194
0195
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
0210 DEFINE_FWK_MODULE(JetResolutionDemo);