PatUtils/interface/ShiftedJetProducerT.h

0001 #ifndef PhysicsTools_PatUtils_ShiftedJetProducerT_h
0002 #define PhysicsTools_PatUtils_ShiftedJetProducerT_h
0003
0004 /** \class ShiftedJetProducerT
0005  *
0006  * Vary energy of jets by +/- 1 standard deviation,
0007  * in order to estimate resulting uncertainty on MET
0008  *
0009  * NOTE: energy scale uncertainties are taken from the Database
0010  *
0011  * \author Christian Veelken, LLR
0012  *
0013  *
0014  *
0015  */
0016
0017 #include "FWCore/Framework/interface/stream/EDProducer.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ParameterSet/interface/FileInPath.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023
0024 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0025 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0026 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0027 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0028 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030
0031 #include "PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h"
0032 #include "PhysicsTools/PatUtils/interface/RawJetExtractorT.h"
0033
0034 #include <string>
0035
0036 template <typename T, typename Textractor>
0037 class ShiftedJetProducerT : public edm::stream::EDProducer<> {
0038   typedef std::vector<T> JetCollection;
0039
0040 public:
0041   explicit ShiftedJetProducerT(const edm::ParameterSet& cfg)
0042       : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0043         src_(cfg.getParameter<edm::InputTag>("src")),
0044         srcToken_(consumes<JetCollection>(src_)),
0045         jetCorrPayloadName_(""),
0046         jetCorrParameters_(nullptr),
0047         jecUncertainty_(nullptr),
0048         jecUncertaintyValue_(-1.) {
0049     if (cfg.exists("jecUncertaintyValue")) {
0050       jecUncertaintyValue_ = cfg.getParameter<double>("jecUncertaintyValue");
0051     } else {
0052       jetCorrUncertaintyTag_ = cfg.getParameter<std::string>("jetCorrUncertaintyTag");
0053       if (cfg.exists("jetCorrInputFileName")) {
0054         jetCorrInputFileName_ = cfg.getParameter<edm::FileInPath>("jetCorrInputFileName");
0055         if (jetCorrInputFileName_.location() == edm::FileInPath::Unknown)
0056           throw cms::Exception("ShiftedJetProducerT")
0057               << " Failed to find JEC parameter file = " << jetCorrInputFileName_ << " !!\n";
0058         jetCorrParameters_ =
0059             std::make_unique<JetCorrectorParameters>(jetCorrInputFileName_.fullPath(), jetCorrUncertaintyTag_);
0060         jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(*jetCorrParameters_);
0061       } else {
0062         jetCorrPayloadName_ = cfg.getParameter<std::string>("jetCorrPayloadName");
0063         jetCorrPayloadToken_ = esConsumes(edm::ESInputTag("", jetCorrPayloadName_));
0064       }
0065     }
0066
0067     addResidualJES_ = cfg.getParameter<bool>("addResidualJES");
0068     if (cfg.exists("jetCorrLabelUpToL3")) {
0069       jetCorrLabelUpToL3_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3");
0070       jetCorrTokenUpToL3_ = mayConsume<reco::JetCorrector>(jetCorrLabelUpToL3_);
0071     }
0072     if (cfg.exists("jetCorrLabelUpToL3Res") && addResidualJES_) {
0073       jetCorrLabelUpToL3Res_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3Res");
0074       jetCorrTokenUpToL3Res_ = mayConsume<reco::JetCorrector>(jetCorrLabelUpToL3Res_);
0075     }
0076     jetCorrEtaMax_ = (cfg.exists("jetCorrEtaMax")) ? cfg.getParameter<double>("jetCorrEtaMax") : 9.9;
0077
0078     shiftBy_ = cfg.getParameter<double>("shiftBy");
0079
0080     verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0081
0082     produces<JetCollection>();
0083   }
0084
0085 private:
0086   void produce(edm::Event& evt, const edm::EventSetup& es) override {
0087     if (verbosity_) {
0088       std::cout << "<ShiftedJetProducerT::produce>:" << std::endl;
0089       std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
0090       std::cout << " src = " << src_.label() << std::endl;
0091     }
0092
0093     edm::Handle<JetCollection> originalJets;
0094     evt.getByToken(srcToken_, originalJets);
0095     edm::Handle<reco::JetCorrector> jetCorrUpToL3;
0096     evt.getByToken(jetCorrTokenUpToL3_, jetCorrUpToL3);
0097     edm::Handle<reco::JetCorrector> jetCorrUpToL3Res;
0098     if (evt.isRealData() && addResidualJES_) {
0099       evt.getByToken(jetCorrTokenUpToL3Res_, jetCorrUpToL3Res);
0100     }
0101     auto shiftedJets = std::make_unique<JetCollection>();
0102
0103     if (!jetCorrPayloadName_.empty()) {
0104       const JetCorrectorParametersCollection& jetCorrParameterSet = es.getData(jetCorrPayloadToken_);
0105       const JetCorrectorParameters& jetCorrParameters = (jetCorrParameterSet)[jetCorrUncertaintyTag_];
0106       jecUncertainty_ = std::make_unique<JetCorrectionUncertainty>(jetCorrParameters);
0107     }
0108
0109     for (typename JetCollection::const_iterator originalJet = originalJets->begin(); originalJet != originalJets->end();
0110          ++originalJet) {
0111       reco::Candidate::LorentzVector originalJetP4 = originalJet->p4();
0112       if (verbosity_) {
0113         std::cout << "originalJet: Pt = " << originalJetP4.pt() << ", eta = " << originalJetP4.eta()
0114                   << ", phi = " << originalJetP4.phi() << std::endl;
0115       }
0116
0117       double shift = 0.;
0118       if (jecUncertaintyValue_ != -1.) {
0119         shift = jecUncertaintyValue_;
0120       } else {
0121         jecUncertainty_->setJetEta(originalJetP4.eta());
0122         jecUncertainty_->setJetPt(originalJetP4.pt());
0123
0124         shift = jecUncertainty_->getUncertainty(true);
0125       }
0126       if (verbosity_) {
0127         std::cout << "shift = " << shift << std::endl;
0128       }
0129
0130       if (evt.isRealData() && addResidualJES_) {
0131         const static pat::RawJetExtractorT<T> rawJetExtractor{};
0132         reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(*originalJet);
0133         if (rawJetP4.E() > 1.e-1) {
0134           reco::Candidate::LorentzVector corrJetP4upToL3 =
0135               std::is_base_of<class PATJetCorrExtractor, Textractor>::value
0136                   ? jetCorrExtractor_(*originalJet, jetCorrLabelUpToL3_.label(), jetCorrEtaMax_, &rawJetP4)
0137                   : jetCorrExtractor_(*originalJet, jetCorrUpToL3.product(), jetCorrEtaMax_, &rawJetP4);
0138           reco::Candidate::LorentzVector corrJetP4upToL3Res =
0139               std::is_base_of<class PATJetCorrExtractor, Textractor>::value
0140                   ? jetCorrExtractor_(*originalJet, jetCorrLabelUpToL3Res_.label(), jetCorrEtaMax_, &rawJetP4)
0141                   : jetCorrExtractor_(*originalJet, jetCorrUpToL3Res.product(), jetCorrEtaMax_, &rawJetP4);
0142           if (corrJetP4upToL3.E() > 1.e-1 && corrJetP4upToL3Res.E() > 1.e-1) {
0143             double residualJES = (corrJetP4upToL3Res.E() / corrJetP4upToL3.E()) - 1.;
0144             shift = sqrt(shift * shift + residualJES * residualJES);
0145           }
0146         }
0147       }
0148
0149       shift *= shiftBy_;
0150       if (verbosity_) {
0151         std::cout << "shift*shiftBy = " << shift << std::endl;
0152       }
0153
0154       T shiftedJet(*originalJet);
0155       shiftedJet.setP4((1. + shift) * originalJetP4);
0156       if (verbosity_) {
0157         std::cout << "shiftedJet: Pt = " << shiftedJet.pt() << ", eta = " << shiftedJet.eta()
0158                   << ", phi = " << shiftedJet.phi() << std::endl;
0159       }
0160
0161       shiftedJets->push_back(shiftedJet);
0162     }
0163
0164     evt.put(std::move(shiftedJets));
0165   }
0166
0167   std::string moduleLabel_;
0168
0169   edm::InputTag src_;
0170   edm::EDGetTokenT<JetCollection> srcToken_;
0171
0172   edm::FileInPath jetCorrInputFileName_;
0173   std::string jetCorrPayloadName_;
0174   edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> jetCorrPayloadToken_;
0175   std::string jetCorrUncertaintyTag_;
0176   std::unique_ptr<JetCorrectorParameters> jetCorrParameters_;
0177   std::unique_ptr<JetCorrectionUncertainty> jecUncertainty_;
0178
0179   bool addResidualJES_;
0180   edm::InputTag jetCorrLabelUpToL3_;                            // L1+L2+L3 correction
0181   edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3_;     // L1+L2+L3 correction
0182   edm::InputTag jetCorrLabelUpToL3Res_;                         // L1+L2+L3+Residual correction
0183   edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3Res_;  // L1+L2+L3+Residual correction
0184   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold (recommended default = 4.7),
0185                           // in order to work around problem with CMSSW_4_2_x JEC factors at high eta,
0186                           // reported in
0187                           //  https://hypernews.cern.ch/HyperNews/CMS/get/jes/270.html
0188                           //  https://hypernews.cern.ch/HyperNews/CMS/get/JetMET/1259/1.html
0189   Textractor jetCorrExtractor_;
0190
0191   double jecUncertaintyValue_;
0192
0193   double shiftBy_;  // set to +1.0/-1.0 for up/down variation of energy scale
0194
0195   int verbosity_;  // flag to enabled/disable debug output
0196 };
0197
0198 #endif