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       produces<edm::ValueMap<float>>("MuonSubtrRawEta");
0101       produces<edm::ValueMap<float>>("MuonSubtrRawPhi");
0102     }
0103   }
0104
0105   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0106     edm::ParameterSetDescription desc;
0107     desc.add<std::string>("@module_label");
0108     desc.add<edm::InputTag>("src");
0109     desc.add<edm::InputTag>("offsetCorrLabel");
0110     desc.add<edm::InputTag>("jetCorrLabel");
0111     desc.add<edm::InputTag>("jetCorrLabelRes");
0112     desc.add<double>("jetCorrEtaMax", 9.9);
0113     desc.add<double>("type1JetPtThreshold");
0114     desc.add<bool>("skipEM");
0115     desc.add<double>("skipEMfractionThreshold");
0116     desc.add<bool>("skipMuons");
0117     desc.add<std::string>("skipMuonSelection");
0118     desc.add<bool>("calcMuonSubtrRawPtAsValueMap", false)
0119         ->setComment(
0120             "calculate muon-subtracted raw pt as a ValueMap for the input collection (only for selected jets, zero for "
0121             "others)");
0122     descriptions.addDefault(desc);
0123   }
0124
0125 private:
0126   void produce(edm::Event& evt, const edm::EventSetup& es) override {
0127     edm::Handle<reco::JetCorrector> jetCorr;
0128     //automatic switch for residual corrections
0129     if (evt.isRealData()) {
0130       jetCorrLabel_ = jetCorrLabelRes_;
0131       evt.getByToken(jetCorrResToken_, jetCorr);
0132     } else {
0133       evt.getByToken(jetCorrToken_, jetCorr);
0134     }
0135
0136     typedef std::vector<T> JetCollection;
0137     edm::Handle<JetCollection> jets;
0138     evt.getByToken(token_, jets);
0139
0140     //new collection
0141     std::unique_ptr<std::vector<T>> cleanedJets(new std::vector<T>());
0142
0143     int numJets = jets->size();
0144     std::vector<float> muonSubtrRawPt(numJets, 0);
0145     std::vector<float> muonSubtrRawEta(numJets, 0);
0146     std::vector<float> muonSubtrRawPhi(numJets, 0);
0147
0148     for (int jetIndex = 0; jetIndex < numJets; ++jetIndex) {
0149       const T& jet = jets->at(jetIndex);
0150
0151       const static JetCleanerForType1MET_namespace::InputTypeCheckerT<T, Textractor> checkInputType{};
0152       checkInputType(jet);
0153
0154       double emEnergyFraction = jet.chargedEmEnergyFraction() + jet.neutralEmEnergyFraction();
0155       if (skipEM_ && emEnergyFraction > skipEMfractionThreshold_)
0156         continue;
0157
0158       const static JetCleanerForType1MET_namespace::RawJetExtractorT<T> rawJetExtractor{};
0159       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
0160
0161       if (skipMuons_) {
0162         const std::vector<reco::CandidatePtr>& cands = jet.daughterPtrVector();
0163         for (std::vector<reco::CandidatePtr>::const_iterator cand = cands.begin(); cand != cands.end(); ++cand) {
0164           const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(cand->get());
0165           const reco::Candidate* mu =
0166               (pfcand != nullptr ? (pfcand->muonRef().isNonnull() ? pfcand->muonRef().get() : nullptr) : cand->get());
0167           if (mu != nullptr && (*skipMuonSelection_)(*mu)) {
0168             reco::Candidate::LorentzVector muonP4 = (*cand)->p4();
0169             rawJetP4 -= muonP4;
0170           }
0171         }
0172       }
0173
0174       reco::Candidate::LorentzVector corrJetP4;
0175       if (checkInputType.isPatJet(jet)) {
0176         corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4);
0177       } else {
0178         corrJetP4 = jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
0179         if (corrJetP4.pt() < type1JetPtThreshold_)
0180           continue;
0181       }
0182
0183       if (corrJetP4.pt() < type1JetPtThreshold_)
0184         continue;
0185
0186       cleanedJets->push_back(jet);
0187       if (calcMuonSubtrRawPtAsValueMap_) {
0188         muonSubtrRawPt[jetIndex] = rawJetP4.Pt();
0189         muonSubtrRawEta[jetIndex] = rawJetP4.Eta();
0190         muonSubtrRawPhi[jetIndex] = rawJetP4.Phi();
0191       }
0192     }
0193
0194     evt.put(std::move(cleanedJets));
0195
0196     if (calcMuonSubtrRawPtAsValueMap_) {
0197       std::unique_ptr<edm::ValueMap<float>> muonSubtrRawPtV(new edm::ValueMap<float>());
0198       edm::ValueMap<float>::Filler fillerMuonSubtrRawPt(*muonSubtrRawPtV);
0199       fillerMuonSubtrRawPt.insert(jets, muonSubtrRawPt.begin(), muonSubtrRawPt.end());
0200       fillerMuonSubtrRawPt.fill();
0201       evt.put(std::move(muonSubtrRawPtV), "MuonSubtrRawPt");
0202
0203       std::unique_ptr<edm::ValueMap<float>> muonSubtrRawEtaV(new edm::ValueMap<float>());
0204       edm::ValueMap<float>::Filler fillerMuonSubtrRawEta(*muonSubtrRawEtaV);
0205       fillerMuonSubtrRawEta.insert(jets, muonSubtrRawEta.begin(), muonSubtrRawEta.end());
0206       fillerMuonSubtrRawEta.fill();
0207       evt.put(std::move(muonSubtrRawEtaV), "MuonSubtrRawEta");
0208
0209       std::unique_ptr<edm::ValueMap<float>> muonSubtrRawPhiV(new edm::ValueMap<float>());
0210       edm::ValueMap<float>::Filler fillerMuonSubtrRawPhi(*muonSubtrRawPhiV);
0211       fillerMuonSubtrRawPhi.insert(jets, muonSubtrRawPhi.begin(), muonSubtrRawPhi.end());
0212       fillerMuonSubtrRawPhi.fill();
0213       evt.put(std::move(muonSubtrRawPhiV), "MuonSubtrRawPhi");
0214     }
0215   }
0216
0217   std::string moduleLabel_;
0218
0219   edm::EDGetTokenT<std::vector<T>> token_;
0220
0221   edm::InputTag offsetCorrLabel_;
0222   edm::EDGetTokenT<reco::JetCorrector> offsetCorrToken_;  // e.g. 'ak5CaloJetL1Fastjet'
0223   edm::InputTag jetCorrLabel_;
0224   edm::InputTag jetCorrLabelRes_;
0225   edm::EDGetTokenT<reco::JetCorrector>
0226       jetCorrToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0227   edm::EDGetTokenT<reco::JetCorrector>
0228       jetCorrResToken_;  // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
0229   Textractor jetCorrExtractor_;
0230
0231   double jetCorrEtaMax_;  // do not use JEC factors for |eta| above this threshold
0232
0233   double type1JetPtThreshold_;  // threshold to distinguish between jets entering Type 1 MET correction
0234   // and jets entering "unclustered energy" sum
0235   // NOTE: threshold is applied on **corrected** jet energy (recommended default = 10 GeV)
0236
0237   bool skipEM_;  // flag to exclude jets with large fraction of electromagnetic energy (electrons/photons)
0238   // from Type 1 + 2 MET corrections
0239   double skipEMfractionThreshold_;
0240
0241   bool skipMuons_;  // flag to subtract momentum of muons (provided muons pass selection cuts) which are within jets
0242   // from jet energy before compute JECs/propagating JECs to Type 1 + 2 MET corrections
0243   std::unique_ptr<StringCutObjectSelector<reco::Candidate>> skipMuonSelection_;
0244
0245   bool calcMuonSubtrRawPtAsValueMap_;  // calculate muon-subtracted raw pt as a ValueMap for the input collection (only for selected jets, zero for others)
0246 };
0247
0248 #endif