Type1MET/interface/JetCleanerForType1METT.h

0001 #ifndef JetMETCorrections_Type1MET_JetCleanerForType1METT_h
0002 #define JetMETCorrections_Type1MET_JetCleanerForType1METT_h
0003
0004 /** \class JetCleanerForType1METT
0005  *
0006  * Clean jets for MET corrections and uncertainties (pt/eta/EM fraction and muons)
0007  * apply also JECs
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 Matthieu Marionneau ETH
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/Utilities/interface/InputTag.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0025
0026 #include "DataFormats/Common/interface/Ref.h"
0027 #include "DataFormats/Common/interface/RefToBase.h"
0028 #include "DataFormats/JetReco/interface/Jet.h"
0029 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0030 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0031 #include "DataFormats/METReco/interface/CorrMETData.h"
0032 #include "DataFormats/PatCandidates/interface/Jet.h"
0033
0034 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0035 #include "DataFormats/MuonReco/interface/Muon.h"
0036
0037 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0038
0039 #include <memory>
0040
0041 #include <string>
0042 #include <type_traits>
0043
0044 namespace JetCleanerForType1MET_namespace {
0045   template <typename T, typename Textractor>
0046   class InputTypeCheckerT {
0047   public:
0048     void operator()(const T&) const {}  // no type-checking needed for reco::PFJet input
0049     bool isPatJet(const T&) const { return std::is_base_of<class pat::Jet, T>::value; }
0050   };
0051
0052   template <typename T>
0053   class RawJetExtractorT  // this template is neccessary to support pat::Jets
0054   // (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
0055   // But it does not handle the muon removal!!!!! MM
0056   {
0057   public:
0058     RawJetExtractorT() {}
0059     reco::Candidate::LorentzVector operator()(const T& jet) const { return jet.p4(); }
0060   };
0061 }  // namespace JetCleanerForType1MET_namespace
0062
0063 template <typename T, typename Textractor>
0064 class JetCleanerForType1METT : public edm::stream::EDProducer<> {
0065 public:
0066   explicit JetCleanerForType1METT(const edm::ParameterSet& cfg)
0067       : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0068         offsetCorrLabel_(""),
0069         skipMuonSelection_(nullptr) {
0070     token_ = consumes<std::vector<T>>(cfg.getParameter<edm::InputTag>("src"));
0071
0072     offsetCorrLabel_ = cfg.getParameter<edm::InputTag>("offsetCorrLabel");
0073     offsetCorrToken_ = consumes<reco::JetCorrector>(offsetCorrLabel_);
0074
0075     jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel");        //for MC
0076     jetCorrLabelRes_ = cfg.getParameter<edm::InputTag>("jetCorrLabelRes");  //for data
0077     jetCorrToken_ = mayConsume<reco::JetCorrector>(jetCorrLabel_);
0078     jetCorrResToken_ = mayConsume<reco::JetCorrector>(jetCorrLabelRes_);
0079
0080     jetCorrEtaMax_ = cfg.getParameter<double>("jetCorrEtaMax");
0081
0082     type1JetPtThreshold_ = cfg.getParameter<double>("type1JetPtThreshold");
0083
0084     skipEM_ = cfg.getParameter<bool>("skipEM");
0085     if (skipEM_) {
0086       skipEMfractionThreshold_ = cfg.getParameter<double>("skipEMfractionThreshold");
0087     }
0088
0089     skipMuons_ = cfg.getParameter<bool>("skipMuons");
0090     if (skipMuons_) {
0091       std::string skipMuonSelection_string = cfg.getParameter<std::string>("skipMuonSelection");
0092       skipMuonSelection_ = std::make_unique<StringCutObjectSelector<reco::Candidate>>(skipMuonSelection_string, true);
0093     }
0094
0095     calcMuonSubtrRawPtAsValueMap_ = cfg.getParameter<bool>("calcMuonSubtrRawPtAsValueMap");
0096
0097     produces<std::vector<T>>();
0098     if (calcMuonSubtrRawPtAsValueMap_)
0099       produces<edm::ValueMap<float>>("MuonSubtrRawPt");
0100   }
0101
0102   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0103     edm::ParameterSetDescription desc;
0104     desc.add<std::string>("@module_label");
0105     desc.add<edm::InputTag>("src");
0106     desc.add<edm::InputTag>("offsetCorrLabel");
0107     desc.add<edm::InputTag>("jetCorrLabel");
0108     desc.add<edm::InputTag>("jetCorrLabelRes");
0109     desc.add<double>("jetCorrEtaMax", 9.9);
0110     desc.add<double>("type1JetPtThreshold");
0111     desc.add<bool>("skipEM");
0112     desc.add<double>("skipEMfractionThreshold");
0113     desc.add<bool>("skipMuons");
0114     desc.add<std::string>("skipMuonSelection");
0115     desc.add<bool>("calcMuonSubtrRawPtAsValueMap", false)
0116         ->setComment(
0117             "calculate muon-subtracted raw pt as a ValueMap for the input collection (only for selected jets, zero for "
0118             "others)");
0119     descriptions.addDefault(desc);
0120   }
0121
0122 private:
0123   void produce(edm::Event& evt, const edm::EventSetup& es) override {
0124     edm::Handle<reco::JetCorrector> jetCorr;
0125     //automatic switch for residual corrections
0126     if (evt.isRealData()) {
0127       jetCorrLabel_ = jetCorrLabelRes_;
0128       evt.getByToken(jetCorrResToken_, jetCorr);
0129     } else {
0130       evt.getByToken(jetCorrToken_, jetCorr);
0131     }
0132
0133     typedef std::vector<T> JetCollection;
0134     edm::Handle<JetCollection> jets;
0135     evt.getByToken(token_, jets);
0136
0137     //new collection
0138     std::unique_ptr<std::vector<T>> cleanedJets(new std::vector<T>());
0139
0140     int numJets = jets->size();
0141     std::vector<float> muonSubtrRawPt(numJets, 0);
0142
0143     for (int jetIndex = 0; jetIndex < numJets; ++jetIndex) {
0144       const T& jet = jets->at(jetIndex);
0145
0146       const static JetCleanerForType1MET_namespace::InputTypeCheckerT<T, Textractor> checkInputType{};
0147       checkInputType(jet);
0148
0149       double emEnergyFraction = jet.chargedEmEnergyFraction() + jet.neutralEmEnergyFraction();
0150       if (skipEM_ && emEnergyFraction > skipEMfractionThreshold_)
0151         continue;
0152
0153       const static JetCleanerForType1MET_namespace::RawJetExtractorT<T> rawJetExtractor{};
0154       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
0155
0156       if (skipMuons_) {
0157         const std::vector<reco::CandidatePtr>& cands = jet.daughterPtrVector();
0158         for (std::vector<reco::CandidatePtr>::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0159           const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(cand->get());
0160           const reco::Candidate* mu =
0161               (pfcand != nullptr ? (pfcand->muonRef().isNonnull() ? pfcand->muonRef().get() : nullptr) : cand->get());
0162           if (mu != nullptr && (*skipMuonSelection_)(*mu)) {
0163             reco::Candidate::LorentzVector muonP4 = (*cand)->p4();
0164             rawJetP4 -= muonP4;
0165           }
0166         }
0167       }
0168
0169       reco::Candidate::LorentzVector corrJetP4;
0170       if (checkInputType.isPatJet(jet)) {
0171         corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4);
0172       } else {
0173         corrJetP4 = jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
0174         if (corrJetP4.pt() < type1JetPtThreshold_)
0175           continue;
0176       }
0177
0178       if (corrJetP4.pt() < type1JetPtThreshold_)
0179         continue;
0180
0181       cleanedJets->push_back(jet);
0182       if (calcMuonSubtrRawPtAsValueMap_)
0183         muonSubtrRawPt[jetIndex] = rawJetP4.Pt();
0184     }
0185
0186     evt.put(std::move(cleanedJets));
0187
0188     if (calcMuonSubtrRawPtAsValueMap_) {
0189       std::unique_ptr<edm::ValueMap<float>> muonSubtrRawPtV(new edm::ValueMap<float>());
0190       edm::ValueMap<float>::Filler fillerMuonSubtrRawPt(*muonSubtrRawPtV);
0191       fillerMuonSubtrRawPt.insert(jets, muonSubtrRawPt.begin(), muonSubtrRawPt.end());
0192       fillerMuonSubtrRawPt.fill();
0193       evt.put(std::move(muonSubtrRawPtV), "MuonSubtrRawPt");
0194     }
0195   }
0196
0197   std::string moduleLabel_;
0198
0199   edm::EDGetTokenT<std::vector<T>> token_;
0200
0201   edm::InputTag offsetCorrLabel_;
0202   edm::EDGetTokenT<reco::JetCorrector> offsetCorrToken_;  // e.g. 'ak5CaloJetL1Fastjet'
0203   edm::InputTag jetCorrLabel_;
0204   edm::InputTag jetCorrLabelRes_;
0205   edm::EDGetTokenT<reco::JetCorrector>
0206       jetCorrToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0207   edm::EDGetTokenT<reco::JetCorrector>
0208       jetCorrResToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0209   Textractor jetCorrExtractor_;
0210
0211   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold
0212
0213   double type1JetPtThreshold_;  // threshold to distinguish between jets entering Type 1 MET correction
0214   // and jets entering "unclustered energy" sum
0215   // NOTE: threshold is applied on **corrected** jet energy (recommended default = 10 GeV)
0216
0217   bool skipEM_;  // flag to exclude jets with large fraction of electromagnetic energy (electrons/photons)
0218   // from Type 1 + 2 MET corrections
0219   double skipEMfractionThreshold_;
0220
0221   bool skipMuons_;  // flag to subtract momentum of muons (provided muons pass selection cuts) which are within jets
0222   // from jet energy before compute JECs/propagating JECs to Type 1 + 2 MET corrections
0223   std::unique_ptr<StringCutObjectSelector<reco::Candidate>> skipMuonSelection_;
0224
0225   bool calcMuonSubtrRawPtAsValueMap_;  // calculate muon-subtracted raw pt as a ValueMap for the input collection (only for selected jets, zero for others)
0226 };
0227
0228 #endif