Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-27 03:18:00 

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/ConfigurationDescriptions.h"
0022 #include "FWCore/ParameterSet/interface/FileInPath.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 
0025 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0026 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0027 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0028 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0029 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 
0032 #include "PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h"
0033 #include "PhysicsTools/PatUtils/interface/RawJetExtractorT.h"
0034 
0035 #include <string>
0036 #include <optional>
0037 
0038 template <typename T, typename Textractor>
0039 class ShiftedJetProducerT : public edm::stream::EDProducer<> {
0040   using JetCollection = std::vector<T>;
0041 
0042 public:
0043   explicit ShiftedJetProducerT(const edm::ParameterSet& cfg)
0044       : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0045         src_(cfg.getParameter<edm::InputTag>("src")),
0046         srcToken_(consumes<JetCollection>(src_)),
0047         jetCorrPayloadName_(""),
0048         jecUncertainty_(),
0049         jecUncertaintyValue_() {
0050     if (cfg.exists("jecUncertaintyValue")) {
0051       jecUncertaintyValue_.emplace(cfg.getParameter<double>("jecUncertaintyValue"));
0052     } else {
0053       jetCorrUncertaintyTag_ = cfg.getParameter<std::string>("jetCorrUncertaintyTag");
0054       if (cfg.exists("jetCorrInputFileName")) {
0055         jetCorrInputFileName_ = cfg.getParameter<edm::FileInPath>("jetCorrInputFileName");
0056         if (jetCorrInputFileName_.location() == edm::FileInPath::Unknown)
0057           throw cms::Exception("ShiftedJetProducerT")
0058               << " Failed to find JEC parameter file = " << jetCorrInputFileName_ << " !!\n";
0059         JetCorrectorParameters jetCorrParameters(jetCorrInputFileName_.fullPath(), jetCorrUncertaintyTag_);
0060         jecUncertainty_.emplace(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 (addResidualJES_) {
0069       if (cfg.exists("jetCorrLabelUpToL3")) {
0070         jetCorrLabelUpToL3_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3");
0071         jetCorrTokenUpToL3_ = consumes<reco::JetCorrector>(jetCorrLabelUpToL3_);
0072       }
0073       if (cfg.exists("jetCorrLabelUpToL3Res")) {
0074         jetCorrLabelUpToL3Res_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3Res");
0075         jetCorrTokenUpToL3Res_ = consumes<reco::JetCorrector>(jetCorrLabelUpToL3Res_);
0076       }
0077     }
0078     jetCorrEtaMax_ = cfg.getParameter<double>("jetCorrEtaMax");
0079 
0080     shiftBy_ = cfg.getParameter<double>("shiftBy");
0081 
0082     verbosity_ = cfg.getUntrackedParameter<int>("verbosity");
0083 
0084     putToken_ = produces<JetCollection>();
0085 
0086     //PATJetCorrExtractor wants a string not a product when called
0087     static_assert(not std::is_base_of<class PATJetCorrExtractor, Textractor>::value);
0088   }
0089 
0090   static void fillDescriptions(edm::ConfigurationDescriptions& iDesc) {
0091     //NOTE: the full complexity of the pset handling could be expressed using
0092     // ifValue and addNode but it might result in present configs failing
0093     edm::ParameterSetDescription ps;
0094     ps.add<edm::InputTag>("src");
0095     ps.addOptional<double>("jecUncertaintyValue");
0096     ps.addOptional<std::string>("jetCorrUncertaintyTag")->setComment("only used if 'jecUncertaintyValue' not declared");
0097     ps.addOptional<edm::FileInPath>("jetCorrInputFileName")
0098         ->setComment("only used if 'jecUncertaintyValue' not declared");
0099     ps.addOptional<std::string>("jetCorrPayloadName")
0100         ->setComment("only used if neither 'jecUncertaintyValue' nor 'jetCorrInputFileName' are declared");
0101 
0102     ps.add<bool>("addResidualJES");
0103     ps.addOptional<edm::InputTag>("jetCorrLabelUpToL3")->setComment("only used of 'addResidualJES' is set true");
0104     ps.addOptional<edm::InputTag>("jetCorrLabelUpToL3Res")->setComment("only used of 'addResidualJES' is set true");
0105 
0106     ps.add<double>("jetCorrEtaMax", 9.9);
0107     ps.add<double>("shiftBy");
0108     ps.addUntracked<int>("verbosity", 0);
0109     iDesc.addDefault(ps);
0110   }
0111 
0112 private:
0113   void produce(edm::Event& evt, const edm::EventSetup& es) override {
0114     if (verbosity_) {
0115       std::cout << "<ShiftedJetProducerT::produce>:" << std::endl;
0116       std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
0117       std::cout << " src = " << src_.label() << std::endl;
0118     }
0119 
0120     JetCollection const& originalJets = evt.get(srcToken_);
0121 
0122     edm::Handle<reco::JetCorrector> jetCorrUpToL3;
0123     edm::Handle<reco::JetCorrector> jetCorrUpToL3Res;
0124     if (evt.isRealData() && addResidualJES_) {
0125       jetCorrUpToL3 = evt.getHandle(jetCorrTokenUpToL3_);
0126       jetCorrUpToL3Res = evt.getHandle(jetCorrTokenUpToL3Res_);
0127     }
0128 
0129     if (not jecUncertaintyValue_) {
0130       if (jetCorrPayloadToken_.isInitialized()) {
0131         auto cacheID = es.get<JetCorrectionsRecord>().cacheIdentifier();
0132         if (cacheID != recordIdentifier_) {
0133           const JetCorrectorParametersCollection& jetCorrParameterSet = es.getData(jetCorrPayloadToken_);
0134           const JetCorrectorParameters& jetCorrParameters = (jetCorrParameterSet)[jetCorrUncertaintyTag_];
0135           jecUncertainty_.emplace(jetCorrParameters);
0136           recordIdentifier_ = cacheID;
0137         }
0138       }
0139     }
0140 
0141     JetCollection shiftedJets;
0142     shiftedJets.reserve(originalJets.size());
0143     for (auto const& originalJet : originalJets) {
0144       reco::Candidate::LorentzVector originalJetP4 = originalJet.p4();
0145       if (verbosity_) {
0146         std::cout << "originalJet: Pt = " << originalJetP4.pt() << ", eta = " << originalJetP4.eta()
0147                   << ", phi = " << originalJetP4.phi() << std::endl;
0148       }
0149 
0150       double shift = 0.;
0151       if (jecUncertaintyValue_) {
0152         shift = *jecUncertaintyValue_;
0153       } else {
0154         jecUncertainty_->setJetEta(originalJetP4.eta());
0155         jecUncertainty_->setJetPt(originalJetP4.pt());
0156 
0157         shift = jecUncertainty_->getUncertainty(true);
0158       }
0159       if (verbosity_) {
0160         std::cout << "shift = " << shift << std::endl;
0161       }
0162 
0163       if (evt.isRealData() && addResidualJES_) {
0164         reco::Candidate::LorentzVector rawJetP4 = pat::RawJetExtractorT<T>{}(originalJet);
0165         if (rawJetP4.E() > 1.e-1) {
0166           reco::Candidate::LorentzVector corrJetP4upToL3 =
0167               jetCorrExtractor_(originalJet, jetCorrUpToL3.product(), jetCorrEtaMax_, &rawJetP4);
0168           reco::Candidate::LorentzVector corrJetP4upToL3Res =
0169               jetCorrExtractor_(originalJet, jetCorrUpToL3Res.product(), jetCorrEtaMax_, &rawJetP4);
0170           if (corrJetP4upToL3.E() > 1.e-1 && corrJetP4upToL3Res.E() > 1.e-1) {
0171             double residualJES = (corrJetP4upToL3Res.E() / corrJetP4upToL3.E()) - 1.;
0172             shift = sqrt(shift * shift + residualJES * residualJES);
0173           }
0174         }
0175       }
0176 
0177       shift *= shiftBy_;
0178       if (verbosity_) {
0179         std::cout << "shift*shiftBy = " << shift << std::endl;
0180       }
0181 
0182       shiftedJets.emplace_back(originalJet);
0183       shiftedJets.back().setP4((1. + shift) * originalJetP4);
0184       if (verbosity_) {
0185         auto const& shiftedJet = shiftedJets.back();
0186         std::cout << "shiftedJet: Pt = " << shiftedJet.pt() << ", eta = " << shiftedJet.eta()
0187                   << ", phi = " << shiftedJet.phi() << std::endl;
0188       }
0189     }
0190 
0191     evt.emplace(putToken_, std::move(shiftedJets));
0192   }
0193 
0194   std::string moduleLabel_;
0195 
0196   edm::InputTag src_;
0197   edm::EDGetTokenT<JetCollection> srcToken_;
0198   edm::EDPutTokenT<JetCollection> putToken_;
0199 
0200   edm::FileInPath jetCorrInputFileName_;
0201   std::string jetCorrPayloadName_;
0202   edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> jetCorrPayloadToken_;
0203   std::string jetCorrUncertaintyTag_;
0204   std::optional<JetCorrectionUncertainty> jecUncertainty_;
0205   unsigned long long recordIdentifier_ = 0;
0206 
0207   bool addResidualJES_;
0208   edm::InputTag jetCorrLabelUpToL3_;                            // L1+L2+L3 correction
0209   edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3_;     // L1+L2+L3 correction
0210   edm::InputTag jetCorrLabelUpToL3Res_;                         // L1+L2+L3+Residual correction
0211   edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3Res_;  // L1+L2+L3+Residual correction
0212   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold (recommended default = 4.7),
0213                           // in order to work around problem with CMSSW_4_2_x JEC factors at high eta,
0214                           // reported in
0215                           //  https://hypernews.cern.ch/HyperNews/CMS/get/jes/270.html
0216                           //  https://hypernews.cern.ch/HyperNews/CMS/get/JetMET/1259/1.html
0217   Textractor jetCorrExtractor_;
0218 
0219   std::optional<double> jecUncertaintyValue_;
0220 
0221   double shiftBy_;  // set to +1.0/-1.0 for up/down variation of energy scale
0222 
0223   int verbosity_;  // flag to enabled/disable debug output
0224 };
0225 
0226 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!