Type1MET/interface/CaloJetMETcorrInputProducerT.h

0001 #ifndef JetMETCorrections_Type1MET_CaloJetMETcorrInputProducer_h
0002 #define JetMETCorrections_Type1MET_CaloJetMETcorrInputProducer_h
0003
0004 /** \class CaloJetMETcorrInputProducer
0005  *
0006  * Produce Type 1 + 2 MET corrections corresponding to differences
0007  * between raw CaloJets and CaloJets with jet energy corrections (JECs) applied
0008  *
0009  * \authors Michael Schmitt, Richard Cavanaugh, The University of Florida
0010  *          Florent Lacroix, University of Illinois at Chicago
0011  *          Christian Veelken, LLR
0012  *
0013  *
0014  *
0015  */
0016
0017 #include "FWCore/Framework/interface/global/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/Utilities/interface/InputTag.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 #include "DataFormats/Common/interface/Ref.h"
0025 #include "DataFormats/Common/interface/RefToBase.h"
0026 #include "DataFormats/JetReco/interface/Jet.h"
0027 #include "DataFormats/Common/interface/View.h"
0028 #include "DataFormats/METReco/interface/MET.h"
0029 #include "DataFormats/PatCandidates/interface/Jet.h"
0030
0031 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0032
0033 #include <string>
0034 #include <type_traits>
0035
0036 namespace CaloJetMETcorrInputProducer_namespace {
0037   template <typename T>
0038   class InputTypeCheckerT {
0039   public:
0040     void operator()(const T&) const {}  // no type-checking needed for reco::CaloJet input
0041     bool isPatJet(const T&) const { return std::is_base_of<class pat::Jet, T>::value; }
0042   };
0043
0044   template <typename T>
0045   class RawJetExtractorT  // this template is neccessary to support pat::Jets
0046                           // (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
0047   {
0048   public:
0049     RawJetExtractorT() {}
0050
0051     reco::Candidate::LorentzVector operator()(const T& jet) const { return jet.p4(); }
0052   };
0053 }  // namespace CaloJetMETcorrInputProducer_namespace
0054
0055 template <typename T, typename Textractor>
0056 class CaloJetMETcorrInputProducerT final : public edm::global::EDProducer<> {
0057 public:
0058   explicit CaloJetMETcorrInputProducerT(const edm::ParameterSet& cfg)
0059       : moduleLabel_(cfg.getParameter<std::string>("@module_label")), offsetCorrLabel_("") {
0060     token_ = consumes<std::vector<T> >(cfg.getParameter<edm::InputTag>("src"));
0061
0062     if (cfg.exists("offsetCorrLabel")) {
0063       offsetCorrLabel_ = cfg.getParameter<edm::InputTag>("offsetCorrLabel");
0064       offsetCorrToken_ = consumes<reco::JetCorrector>(offsetCorrLabel_);
0065     }
0066     jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel");
0067     jetCorrToken_ = consumes<reco::JetCorrector>(jetCorrLabel_);
0068
0069     jetCorrEtaMax_ = (cfg.exists("jetCorrEtaMax")) ? cfg.getParameter<double>("jetCorrEtaMax") : 9.9;
0070
0071     type1JetPtThreshold_ = cfg.getParameter<double>("type1JetPtThreshold");
0072
0073     skipEM_ = cfg.getParameter<bool>("skipEM");
0074     if (skipEM_) {
0075       skipEMfractionThreshold_ = cfg.getParameter<double>("skipEMfractionThreshold");
0076     }
0077
0078     if (cfg.exists("srcMET")) {
0079       srcMET_ = cfg.getParameter<edm::InputTag>("srcMET");
0080       metToken_ = consumes<edm::View<reco::MET> >(cfg.getParameter<edm::InputTag>("srcMET"));
0081     }
0082
0083     produces<CorrMETData>("type1");
0084     produces<CorrMETData>("type2");
0085     produces<CorrMETData>("offset");
0086   }
0087   ~CaloJetMETcorrInputProducerT() override {}
0088
0089 private:
0090   void produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const override {
0091     std::unique_ptr<CorrMETData> type1Correction(new CorrMETData());
0092     std::unique_ptr<CorrMETData> unclEnergySum(new CorrMETData());
0093     std::unique_ptr<CorrMETData> offsetEnergySum(new CorrMETData());
0094
0095     edm::Handle<reco::JetCorrector> jetCorr;
0096     evt.getByToken(jetCorrToken_, jetCorr);
0097
0098     typedef std::vector<T> JetCollection;
0099     edm::Handle<JetCollection> jets;
0100     evt.getByToken(token_, jets);
0101
0102     typedef edm::View<reco::MET> METView;
0103     edm::Handle<METView> met;
0104     if (!srcMET_.label().empty()) {
0105       evt.getByToken(metToken_, met);
0106       if (met->size() != 1)
0107         throw cms::Exception("CaloJetMETcorrInputProducer::produce")
0108             << "Failed to find unique MET in the event, src = " << srcMET_.label() << " !!\n";
0109
0110       //--- compute "unclustered energy" by sutracting from the reconstructed MET
0111       //   (i.e. from the negative vectorial sum of all particles reconstructed in the event)
0112       //    the momenta of (high Pt) jets which enter Type 1 MET corrections
0113       //
0114       //    NOTE: MET = -(jets + muons + "unclustered energy"),
0115       //          so "unclustered energy" = -(MET + jets + muons),
0116       //          i.e. (high Pt) jets enter the sum of "unclustered energy" with negative sign.
0117       //
0118       unclEnergySum->mex = -met->front().px();
0119       unclEnergySum->mey = -met->front().py();
0120       unclEnergySum->sumet = met->front().sumEt();
0121     }
0122
0123     int numJets = jets->size();
0124     for (int jetIndex = 0; jetIndex < numJets; ++jetIndex) {
0125       const T& rawJet = jets->at(jetIndex);
0126
0127       const CaloJetMETcorrInputProducer_namespace::InputTypeCheckerT<T> checkInputType{};
0128       checkInputType(rawJet);
0129
0130       const CaloJetMETcorrInputProducer_namespace::RawJetExtractorT<T> rawJetExtractor;
0131       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(rawJet);
0132
0133       reco::Candidate::LorentzVector corrJetP4;
0134       if (checkInputType.isPatJet(rawJet))
0135         corrJetP4 = jetCorrExtractor_(rawJet, jetCorrLabel_.label(), jetCorrEtaMax_);
0136       else
0137         corrJetP4 = jetCorrExtractor_(rawJet, jetCorr.product(), jetCorrEtaMax_);
0138
0139       if (corrJetP4.pt() > type1JetPtThreshold_) {
0140         unclEnergySum->mex -= rawJetP4.px();
0141         unclEnergySum->mey -= rawJetP4.py();
0142         unclEnergySum->sumet -= rawJetP4.Et();
0143
0144         if (skipEM_ && rawJet.emEnergyFraction() > skipEMfractionThreshold_)
0145           continue;
0146
0147         reco::Candidate::LorentzVector rawJetP4offsetCorr = rawJetP4;
0148         if (!offsetCorrLabel_.label().empty()) {
0149           edm::Handle<reco::JetCorrector> offsetCorr;
0150           evt.getByToken(offsetCorrToken_, offsetCorr);
0151           if (checkInputType.isPatJet(rawJet))
0152             rawJetP4offsetCorr = jetCorrExtractor_(rawJet, offsetCorrLabel_.label(), jetCorrEtaMax_);
0153           else
0154             rawJetP4offsetCorr = jetCorrExtractor_(rawJet, offsetCorr.product(), jetCorrEtaMax_);
0155
0156           offsetEnergySum->mex += (rawJetP4.px() - rawJetP4offsetCorr.px());
0157           offsetEnergySum->mey += (rawJetP4.py() - rawJetP4offsetCorr.py());
0158           offsetEnergySum->sumet += (rawJetP4.Et() - rawJetP4offsetCorr.Et());
0159         }
0160
0161         //--- MET balances momentum of reconstructed particles,
0162         //    hence correction to jets and corresponding Type 1 MET correction are of opposite sign
0163         type1Correction->mex -= (corrJetP4.px() - rawJetP4offsetCorr.px());
0164         type1Correction->mey -= (corrJetP4.py() - rawJetP4offsetCorr.py());
0165         type1Correction->sumet += (corrJetP4.Et() - rawJetP4offsetCorr.Et());
0166       }
0167     }
0168
0169     //--- add
0170     //     o Type 1 MET correction                (difference corrected-uncorrected jet energy for jets of (corrected) Pt > 20 GeV)
0171     //     o momentum sum of "unclustered energy" (jets of (corrected) Pt < 20 GeV)
0172     //     o momentum sum of "offset energy"      (sum of energy attributed to pile-up/underlying event)
0173     //    to the event
0174     evt.put(std::move(type1Correction), "type1");
0175     evt.put(std::move(unclEnergySum), "type2");
0176     evt.put(std::move(offsetEnergySum), "offset");
0177   }
0178
0179   std::string moduleLabel_;
0180
0181   edm::EDGetTokenT<std::vector<T> > token_;
0182
0183   edm::InputTag offsetCorrLabel_;
0184   edm::EDGetTokenT<reco::JetCorrector> offsetCorrToken_;  // e.g. 'ak5CaloJetL1Fastjet'
0185   edm::InputTag jetCorrLabel_;
0186   edm::EDGetTokenT<reco::JetCorrector>
0187       jetCorrToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0188   Textractor jetCorrExtractor_;
0189
0190   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold (recommended default = 4.7),
0191                           // in order to work around problem with CMSSW_4_2_x JEC factors at high eta,
0192                           // reported in
0193                           //  https://hypernews.cern.ch/HyperNews/CMS/get/jes/270.html
0194                           //  https://hypernews.cern.ch/HyperNews/CMS/get/JetMET/1259/1.html
0195
0196   double type1JetPtThreshold_;  // threshold to distinguish between jets entering Type 1 MET correction
0197                                 // and jets entering "unclustered energy" sum
0198                                 // NOTE: threshold is applied on **corrected** jet energy (recommended default = 10 GeV)
0199
0200   bool skipEM_;  // flag to exclude jets with large fraction of electromagnetic energy (electrons/photons)
0201                  // from Type 1 + 2 MET corrections
0202   double skipEMfractionThreshold_;
0203
0204   edm::InputTag srcMET_;  // MET input, needed to compute "unclustered energy" sum for Type 2 MET correction
0205   edm::EDGetTokenT<edm::View<reco::MET> > metToken_;
0206 };
0207
0208 #endif