File indexing completed on 2024-04-06 12:19:26
0001 #ifndef JetMETCorrections_Type1MET_PFJetMETcorrInputProducerT_h
0002 #define JetMETCorrections_Type1MET_PFJetMETcorrInputProducerT_h
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
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 {}
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
0056
0057
0058 {
0059 public:
0060 RawJetExtractorT() {}
0061 reco::Candidate::LorentzVector operator()(const T& jet) const { return jet.p4(); }
0062 };
0063
0064
0065
0066
0067 template <typename T>
0068 class AdditionalScalesT {
0069 public:
0070 AdditionalScalesT() {}
0071 float operator()(const T& jet) const { return 1.0; }
0072 };
0073
0074 }
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");
0084 jetCorrLabelRes_ = cfg.getParameter<edm::InputTag>("jetCorrLabelRes");
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
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
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
0229
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
0247
0248
0249
0250
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_;
0266 edm::InputTag jetCorrLabel_;
0267 edm::InputTag jetCorrLabelRes_;
0268 edm::EDGetTokenT<reco::JetCorrector>
0269 jetCorrToken_;
0270 edm::EDGetTokenT<reco::JetCorrector>
0271 jetCorrResToken_;
0272 Textractor jetCorrExtractor_;
0273
0274 double jetCorrEtaMax_;
0275
0276 double type1JetPtThreshold_;
0277
0278
0279
0280 bool skipEM_;
0281
0282 double skipEMfractionThreshold_;
0283
0284 bool skipMuons_;
0285
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