Type1MET/interface/PFJetMETcorrInputProducerT.h

0001 #ifndef JetMETCorrections_Type1MET_PFJetMETcorrInputProducerT_h
0002 #define JetMETCorrections_Type1MET_PFJetMETcorrInputProducerT_h
0003
0004 /** \class PFJetMETcorrInputProducerT
0005  *
0006  * Produce Type 1 + 2 MET corrections corresponding to differences
0007  * between raw PFJets and PFJets with jet energy corrections (JECs) applied
0008  *
0009  * NOTE: class is templated to that it works with reco::PFJets as well as with pat::Jets of PF-type as input
0010  *
0011  * \authors Michael Schmitt, Richard Cavanaugh, The University of Florida
0012  *          Florent Lacroix, University of Illinois at Chicago
0013  *          Christian Veelken, LLR
0014  *
0015  *
0016  *
0017  */
0018
0019 #include "FWCore/Framework/interface/stream/EDProducer.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "FWCore/Utilities/interface/Exception.h"
0025 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0027
0028 #include "DataFormats/Common/interface/Ref.h"
0029 #include "DataFormats/Common/interface/RefToBase.h"
0030 #include "DataFormats/JetReco/interface/Jet.h"
0031 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0032 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0033 #include "DataFormats/METReco/interface/CorrMETData.h"
0034 #include "DataFormats/PatCandidates/interface/Jet.h"
0035
0036 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0037 #include "DataFormats/MuonReco/interface/Muon.h"
0038
0039 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0040
0041 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0042
0043 #include <string>
0044 #include <type_traits>
0045
0046 namespace PFJetMETcorrInputProducer_namespace {
0047   template <typename T, typename Textractor>
0048   class InputTypeCheckerT {
0049   public:
0050     void operator()(const T&) const {}  // no type-checking needed for reco::PFJet input
0051     bool isPatJet(const T&) const { return std::is_base_of<class pat::Jet, T>::value; }
0052   };
0053
0054   template <typename T>
0055   class RawJetExtractorT  // this template is neccessary to support pat::Jets
0056                           // (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
0057   // But it does not handle the muon removal!!!!! MM
0058   {
0059   public:
0060     RawJetExtractorT() {}
0061     reco::Candidate::LorentzVector operator()(const T& jet) const { return jet.p4(); }
0062   };
0063
0064   // functor to retrieve additional scales of jets
0065   // general template just returns 1 because only pat::Jet keeps track of the scales
0066   // specialized template for pat::Jet returns combined factor of additional scales
0067   template <typename T>
0068   class AdditionalScalesT {
0069   public:
0070     AdditionalScalesT() {}
0071     float operator()(const T& jet) const { return 1.0; }
0072   };
0073
0074 }  // namespace PFJetMETcorrInputProducer_namespace
0075
0076 template <typename T, typename Textractor>
0077 class PFJetMETcorrInputProducerT : public edm::stream::EDProducer<> {
0078 public:
0079   explicit PFJetMETcorrInputProducerT(const edm::ParameterSet& cfg) : skipMuonSelection_(nullptr) {
0080     token_ = consumes<std::vector<T> >(cfg.getParameter<edm::InputTag>("src"));
0081     offsetCorrLabel_ = cfg.getParameter<edm::InputTag>("offsetCorrLabel");
0082     offsetCorrToken_ = consumes<reco::JetCorrector>(offsetCorrLabel_);
0083     jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel");        //for MC
0084     jetCorrLabelRes_ = cfg.getParameter<edm::InputTag>("jetCorrLabelRes");  //for data
0085     jetCorrToken_ = mayConsume<reco::JetCorrector>(jetCorrLabel_);
0086     jetCorrResToken_ = mayConsume<reco::JetCorrector>(jetCorrLabelRes_);
0087
0088     jetCorrEtaMax_ = cfg.getParameter<double>("jetCorrEtaMax");
0089
0090     type1JetPtThreshold_ = cfg.getParameter<double>("type1JetPtThreshold");
0091
0092     skipEM_ = cfg.getParameter<bool>("skipEM");
0093     if (skipEM_) {
0094       skipEMfractionThreshold_ = cfg.getParameter<double>("skipEMfractionThreshold");
0095     }
0096
0097     skipMuons_ = cfg.getParameter<bool>("skipMuons");
0098     if (skipMuons_) {
0099       std::string skipMuonSelection_string = cfg.getParameter<std::string>("skipMuonSelection");
0100       skipMuonSelection_ = new StringCutObjectSelector<reco::Candidate>(skipMuonSelection_string, true);
0101     }
0102
0103     if (cfg.exists("type2Binning")) {
0104       typedef std::vector<edm::ParameterSet> vParameterSet;
0105       vParameterSet cfgType2Binning = cfg.getParameter<vParameterSet>("type2Binning");
0106       for (vParameterSet::const_iterator cfgType2BinningEntry = cfgType2Binning.begin();
0107            cfgType2BinningEntry != cfgType2Binning.end();
0108            ++cfgType2BinningEntry) {
0109         type2Binning_.push_back(new type2BinningEntryType(*cfgType2BinningEntry));
0110       }
0111     } else {
0112       type2Binning_.push_back(new type2BinningEntryType());
0113     }
0114
0115     produces<CorrMETData>("type1");
0116     for (typename std::vector<type2BinningEntryType*>::const_iterator type2BinningEntry = type2Binning_.begin();
0117          type2BinningEntry != type2Binning_.end();
0118          ++type2BinningEntry) {
0119       produces<CorrMETData>((*type2BinningEntry)->getInstanceLabel_full("type2"));
0120       produces<CorrMETData>((*type2BinningEntry)->getInstanceLabel_full("offset"));
0121     }
0122   }
0123   ~PFJetMETcorrInputProducerT() override {
0124     delete skipMuonSelection_;
0125
0126     for (typename std::vector<type2BinningEntryType*>::const_iterator it = type2Binning_.begin();
0127          it != type2Binning_.end();
0128          ++it) {
0129       delete (*it);
0130     }
0131   }
0132
0133   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0134     edm::ParameterSetDescription desc;
0135     desc.add<edm::InputTag>("src", edm::InputTag("ak4PFJetsCHS"));
0136     desc.add<edm::InputTag>("offsetCorrLabel", edm::InputTag("ak4PFCHSL1FastjetCorrector"));
0137     desc.add<edm::InputTag>("jetCorrLabel", edm::InputTag("ak4PFCHSL1FastL2L3Corrector"));
0138     desc.add<edm::InputTag>("jetCorrLabelRes", edm::InputTag("ak4PFCHSL1FastL2L3ResidualCorrector"));
0139     desc.add<double>("jetCorrEtaMax", 9.9);
0140     desc.add<double>("type1JetPtThreshold", 15);
0141     desc.add<bool>("skipEM", true);
0142     desc.add<double>("skipEMfractionThreshold", 0.90);
0143     desc.add<bool>("skipMuons", true);
0144     desc.add<std::string>("skipMuonSelection", "isGlobalMuon | isStandAloneMuon");
0145     descriptions.add(defaultModuleLabel<PFJetMETcorrInputProducerT<T, Textractor> >(), desc);
0146   }
0147
0148 private:
0149   void produce(edm::Event& evt, const edm::EventSetup& es) override {
0150     std::unique_ptr<CorrMETData> type1Correction(new CorrMETData());
0151     for (typename std::vector<type2BinningEntryType*>::iterator type2BinningEntry = type2Binning_.begin();
0152          type2BinningEntry != type2Binning_.end();
0153          ++type2BinningEntry) {
0154       (*type2BinningEntry)->binUnclEnergySum_ = CorrMETData();
0155       (*type2BinningEntry)->binOffsetEnergySum_ = CorrMETData();
0156     }
0157
0158     edm::Handle<reco::JetCorrector> jetCorr;
0159     //automatic switch for residual corrections
0160     if (evt.isRealData()) {
0161       jetCorrLabel_ = jetCorrLabelRes_;
0162       evt.getByToken(jetCorrResToken_, jetCorr);
0163     } else {
0164       evt.getByToken(jetCorrToken_, jetCorr);
0165     }
0166
0167     typedef std::vector<T> JetCollection;
0168     edm::Handle<JetCollection> jets;
0169     evt.getByToken(token_, jets);
0170
0171     int numJets = jets->size();
0172     for (int jetIndex = 0; jetIndex < numJets; ++jetIndex) {
0173       const T& jet = jets->at(jetIndex);
0174
0175       const static PFJetMETcorrInputProducer_namespace::InputTypeCheckerT<T, Textractor> checkInputType{};
0176       checkInputType(jet);
0177
0178       double emEnergyFraction = jet.chargedEmEnergyFraction() + jet.neutralEmEnergyFraction();
0179       if (skipEM_ && emEnergyFraction > skipEMfractionThreshold_)
0180         continue;
0181
0182       const static PFJetMETcorrInputProducer_namespace::RawJetExtractorT<T> rawJetExtractor{};
0183       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
0184
0185       if (skipMuons_) {
0186         const std::vector<reco::CandidatePtr>& cands = jet.daughterPtrVector();
0187         for (std::vector<reco::CandidatePtr>::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0188           const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(cand->get());
0189           const reco::Candidate* mu =
0190               (pfcand != nullptr ? (pfcand->muonRef().isNonnull() ? pfcand->muonRef().get() : nullptr) : cand->get());
0191           if (mu != nullptr && (*skipMuonSelection_)(*mu)) {
0192             reco::Candidate::LorentzVector muonP4 = (*cand)->p4();
0193             rawJetP4 -= muonP4;
0194           }
0195         }
0196       }
0197
0198       reco::Candidate::LorentzVector corrJetP4;
0199       if (checkInputType.isPatJet(jet))
0200         corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4);
0201       else
0202         corrJetP4 = jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
0203       // retrieve additional jet energy scales in case of pat::Jets (done via specialized template defined for pat::Jets) and apply them
0204       const static PFJetMETcorrInputProducer_namespace::AdditionalScalesT<T> additionalScales{};
0205       corrJetP4 *= additionalScales(jet);
0206
0207       if (corrJetP4.pt() > type1JetPtThreshold_) {
0208         reco::Candidate::LorentzVector rawJetP4offsetCorr = rawJetP4;
0209         if (!offsetCorrLabel_.label().empty()) {
0210           edm::Handle<reco::JetCorrector> offsetCorr;
0211           evt.getByToken(offsetCorrToken_, offsetCorr);
0212           if (checkInputType.isPatJet(jet))
0213             rawJetP4offsetCorr = jetCorrExtractor_(jet, offsetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4);
0214           else
0215             rawJetP4offsetCorr = jetCorrExtractor_(jet, offsetCorr.product(), jetCorrEtaMax_, &rawJetP4);
0216
0217           for (typename std::vector<type2BinningEntryType*>::iterator type2BinningEntry = type2Binning_.begin();
0218                type2BinningEntry != type2Binning_.end();
0219                ++type2BinningEntry) {
0220             if (!(*type2BinningEntry)->binSelection_ || (*(*type2BinningEntry)->binSelection_)(corrJetP4)) {
0221               (*type2BinningEntry)->binOffsetEnergySum_.mex += (rawJetP4.px() - rawJetP4offsetCorr.px());
0222               (*type2BinningEntry)->binOffsetEnergySum_.mey += (rawJetP4.py() - rawJetP4offsetCorr.py());
0223               (*type2BinningEntry)->binOffsetEnergySum_.sumet += (rawJetP4.Et() - rawJetP4offsetCorr.Et());
0224             }
0225           }
0226         }
0227
0228         //--- MET balances momentum of reconstructed particles,
0229         //    hence correction to jets and corresponding Type 1 MET correction are of opposite sign
0230         type1Correction->mex -= (corrJetP4.px() - rawJetP4offsetCorr.px());
0231         type1Correction->mey -= (corrJetP4.py() - rawJetP4offsetCorr.py());
0232         type1Correction->sumet += (corrJetP4.Et() - rawJetP4offsetCorr.Et());
0233       } else {
0234         for (typename std::vector<type2BinningEntryType*>::iterator type2BinningEntry = type2Binning_.begin();
0235              type2BinningEntry != type2Binning_.end();
0236              ++type2BinningEntry) {
0237           if (!(*type2BinningEntry)->binSelection_ || (*(*type2BinningEntry)->binSelection_)(corrJetP4)) {
0238             (*type2BinningEntry)->binUnclEnergySum_.mex += rawJetP4.px();
0239             (*type2BinningEntry)->binUnclEnergySum_.mey += rawJetP4.py();
0240             (*type2BinningEntry)->binUnclEnergySum_.sumet += rawJetP4.Et();
0241           }
0242         }
0243       }
0244     }
0245
0246     //--- add
0247     //     o Type 1 MET correction                (difference corrected-uncorrected jet energy for jets of (corrected) Pt > 10 GeV)
0248     //     o momentum sum of "unclustered energy" (jets of (corrected) Pt < 10 GeV)
0249     //     o momentum sum of "offset energy"      (sum of energy attributed to pile-up/underlying event)
0250     //    to the event
0251     evt.put(std::move(type1Correction), "type1");
0252     for (typename std::vector<type2BinningEntryType*>::const_iterator type2BinningEntry = type2Binning_.begin();
0253          type2BinningEntry != type2Binning_.end();
0254          ++type2BinningEntry) {
0255       evt.put(std::unique_ptr<CorrMETData>(new CorrMETData((*type2BinningEntry)->binUnclEnergySum_)),
0256               (*type2BinningEntry)->getInstanceLabel_full("type2"));
0257       evt.put(std::unique_ptr<CorrMETData>(new CorrMETData((*type2BinningEntry)->binOffsetEnergySum_)),
0258               (*type2BinningEntry)->getInstanceLabel_full("offset"));
0259     }
0260   }
0261
0262   edm::EDGetTokenT<std::vector<T> > token_;
0263
0264   edm::InputTag offsetCorrLabel_;
0265   edm::EDGetTokenT<reco::JetCorrector> offsetCorrToken_;  // e.g. 'ak5CaloJetL1Fastjet'
0266   edm::InputTag jetCorrLabel_;
0267   edm::InputTag jetCorrLabelRes_;
0268   edm::EDGetTokenT<reco::JetCorrector>
0269       jetCorrToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0270   edm::EDGetTokenT<reco::JetCorrector>
0271       jetCorrResToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0272   Textractor jetCorrExtractor_;
0273
0274   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold
0275
0276   double type1JetPtThreshold_;  // threshold to distinguish between jets entering Type 1 MET correction
0277                                 // and jets entering "unclustered energy" sum
0278                                 // NOTE: threshold is applied on **corrected** jet energy (recommended default = 10 GeV)
0279
0280   bool skipEM_;  // flag to exclude jets with large fraction of electromagnetic energy (electrons/photons)
0281                  // from Type 1 + 2 MET corrections
0282   double skipEMfractionThreshold_;
0283
0284   bool skipMuons_;  // flag to subtract momentum of muons (provided muons pass selection cuts) which are within jets
0285                     // from jet energy before compute JECs/propagating JECs to Type 1 + 2 MET corrections
0286   StringCutObjectSelector<reco::Candidate>* skipMuonSelection_;
0287
0288   struct type2BinningEntryType {
0289     type2BinningEntryType() : binLabel_(""), binSelection_(nullptr) {}
0290     type2BinningEntryType(const edm::ParameterSet& cfg)
0291         : binLabel_(cfg.getParameter<std::string>("binLabel")),
0292           binSelection_(new StringCutObjectSelector<reco::Candidate::LorentzVector>(
0293               cfg.getParameter<std::string>("binSelection"))) {}
0294     ~type2BinningEntryType() { delete binSelection_; }
0295     std::string getInstanceLabel_full(const std::string& instanceLabel) {
0296       std::string retVal = instanceLabel;
0297       if (!instanceLabel.empty() && !binLabel_.empty())
0298         retVal.append("#");
0299       retVal.append(binLabel_);
0300       return retVal;
0301     }
0302     std::string binLabel_;
0303     StringCutObjectSelector<reco::Candidate::LorentzVector>* binSelection_;
0304     CorrMETData binUnclEnergySum_;
0305     CorrMETData binOffsetEnergySum_;
0306   };
0307   std::vector<type2BinningEntryType*> type2Binning_;
0308 };
0309
0310 #endif