File indexing completed on 2024-11-14 04:15:41
0001 #ifndef JetMETCorrections_Type1MET_JetCleanerForType1METT_h
0002 #define JetMETCorrections_Type1MET_JetCleanerForType1METT_h
0003
0004
0005
0006
0007
0008
0009
0010
0011
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 {}
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
0054
0055
0056 {
0057 public:
0058 RawJetExtractorT() {}
0059 reco::Candidate::LorentzVector operator()(const T& jet) const { return jet.p4(); }
0060 };
0061 }
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");
0076 jetCorrLabelRes_ = cfg.getParameter<edm::InputTag>("jetCorrLabelRes");
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
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
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_;
0223 edm::InputTag jetCorrLabel_;
0224 edm::InputTag jetCorrLabelRes_;
0225 edm::EDGetTokenT<reco::JetCorrector>
0226 jetCorrToken_;
0227 edm::EDGetTokenT<reco::JetCorrector>
0228 jetCorrResToken_;
0229 Textractor jetCorrExtractor_;
0230
0231 double jetCorrEtaMax_;
0232
0233 double type1JetPtThreshold_;
0234
0235
0236
0237 bool skipEM_;
0238
0239 double skipEMfractionThreshold_;
0240
0241 bool skipMuons_;
0242
0243 std::unique_ptr<StringCutObjectSelector<reco::Candidate>> skipMuonSelection_;
0244
0245 bool calcMuonSubtrRawPtAsValueMap_;
0246 };
0247
0248 #endif