File indexing completed on 2024-11-27 03:18:00
0001 #ifndef PhysicsTools_PatUtils_ShiftedJetProducerT_h
0002 #define PhysicsTools_PatUtils_ShiftedJetProducerT_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/ParameterSet/interface/ConfigurationDescriptions.h"
0022 #include "FWCore/ParameterSet/interface/FileInPath.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024
0025 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0026 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0027 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
0028 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0029 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031
0032 #include "PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h"
0033 #include "PhysicsTools/PatUtils/interface/RawJetExtractorT.h"
0034
0035 #include <string>
0036 #include <optional>
0037
0038 template <typename T, typename Textractor>
0039 class ShiftedJetProducerT : public edm::stream::EDProducer<> {
0040 using JetCollection = std::vector<T>;
0041
0042 public:
0043 explicit ShiftedJetProducerT(const edm::ParameterSet& cfg)
0044 : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0045 src_(cfg.getParameter<edm::InputTag>("src")),
0046 srcToken_(consumes<JetCollection>(src_)),
0047 jetCorrPayloadName_(""),
0048 jecUncertainty_(),
0049 jecUncertaintyValue_() {
0050 if (cfg.exists("jecUncertaintyValue")) {
0051 jecUncertaintyValue_.emplace(cfg.getParameter<double>("jecUncertaintyValue"));
0052 } else {
0053 jetCorrUncertaintyTag_ = cfg.getParameter<std::string>("jetCorrUncertaintyTag");
0054 if (cfg.exists("jetCorrInputFileName")) {
0055 jetCorrInputFileName_ = cfg.getParameter<edm::FileInPath>("jetCorrInputFileName");
0056 if (jetCorrInputFileName_.location() == edm::FileInPath::Unknown)
0057 throw cms::Exception("ShiftedJetProducerT")
0058 << " Failed to find JEC parameter file = " << jetCorrInputFileName_ << " !!\n";
0059 JetCorrectorParameters jetCorrParameters(jetCorrInputFileName_.fullPath(), jetCorrUncertaintyTag_);
0060 jecUncertainty_.emplace(jetCorrParameters);
0061 } else {
0062 jetCorrPayloadName_ = cfg.getParameter<std::string>("jetCorrPayloadName");
0063 jetCorrPayloadToken_ = esConsumes(edm::ESInputTag("", jetCorrPayloadName_));
0064 }
0065 }
0066
0067 addResidualJES_ = cfg.getParameter<bool>("addResidualJES");
0068 if (addResidualJES_) {
0069 if (cfg.exists("jetCorrLabelUpToL3")) {
0070 jetCorrLabelUpToL3_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3");
0071 jetCorrTokenUpToL3_ = consumes<reco::JetCorrector>(jetCorrLabelUpToL3_);
0072 }
0073 if (cfg.exists("jetCorrLabelUpToL3Res")) {
0074 jetCorrLabelUpToL3Res_ = cfg.getParameter<edm::InputTag>("jetCorrLabelUpToL3Res");
0075 jetCorrTokenUpToL3Res_ = consumes<reco::JetCorrector>(jetCorrLabelUpToL3Res_);
0076 }
0077 }
0078 jetCorrEtaMax_ = cfg.getParameter<double>("jetCorrEtaMax");
0079
0080 shiftBy_ = cfg.getParameter<double>("shiftBy");
0081
0082 verbosity_ = cfg.getUntrackedParameter<int>("verbosity");
0083
0084 putToken_ = produces<JetCollection>();
0085
0086
0087 static_assert(not std::is_base_of<class PATJetCorrExtractor, Textractor>::value);
0088 }
0089
0090 static void fillDescriptions(edm::ConfigurationDescriptions& iDesc) {
0091
0092
0093 edm::ParameterSetDescription ps;
0094 ps.add<edm::InputTag>("src");
0095 ps.addOptional<double>("jecUncertaintyValue");
0096 ps.addOptional<std::string>("jetCorrUncertaintyTag")->setComment("only used if 'jecUncertaintyValue' not declared");
0097 ps.addOptional<edm::FileInPath>("jetCorrInputFileName")
0098 ->setComment("only used if 'jecUncertaintyValue' not declared");
0099 ps.addOptional<std::string>("jetCorrPayloadName")
0100 ->setComment("only used if neither 'jecUncertaintyValue' nor 'jetCorrInputFileName' are declared");
0101
0102 ps.add<bool>("addResidualJES");
0103 ps.addOptional<edm::InputTag>("jetCorrLabelUpToL3")->setComment("only used of 'addResidualJES' is set true");
0104 ps.addOptional<edm::InputTag>("jetCorrLabelUpToL3Res")->setComment("only used of 'addResidualJES' is set true");
0105
0106 ps.add<double>("jetCorrEtaMax", 9.9);
0107 ps.add<double>("shiftBy");
0108 ps.addUntracked<int>("verbosity", 0);
0109 iDesc.addDefault(ps);
0110 }
0111
0112 private:
0113 void produce(edm::Event& evt, const edm::EventSetup& es) override {
0114 if (verbosity_) {
0115 std::cout << "<ShiftedJetProducerT::produce>:" << std::endl;
0116 std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
0117 std::cout << " src = " << src_.label() << std::endl;
0118 }
0119
0120 JetCollection const& originalJets = evt.get(srcToken_);
0121
0122 edm::Handle<reco::JetCorrector> jetCorrUpToL3;
0123 edm::Handle<reco::JetCorrector> jetCorrUpToL3Res;
0124 if (evt.isRealData() && addResidualJES_) {
0125 jetCorrUpToL3 = evt.getHandle(jetCorrTokenUpToL3_);
0126 jetCorrUpToL3Res = evt.getHandle(jetCorrTokenUpToL3Res_);
0127 }
0128
0129 if (not jecUncertaintyValue_) {
0130 if (jetCorrPayloadToken_.isInitialized()) {
0131 auto cacheID = es.get<JetCorrectionsRecord>().cacheIdentifier();
0132 if (cacheID != recordIdentifier_) {
0133 const JetCorrectorParametersCollection& jetCorrParameterSet = es.getData(jetCorrPayloadToken_);
0134 const JetCorrectorParameters& jetCorrParameters = (jetCorrParameterSet)[jetCorrUncertaintyTag_];
0135 jecUncertainty_.emplace(jetCorrParameters);
0136 recordIdentifier_ = cacheID;
0137 }
0138 }
0139 }
0140
0141 JetCollection shiftedJets;
0142 shiftedJets.reserve(originalJets.size());
0143 for (auto const& originalJet : originalJets) {
0144 reco::Candidate::LorentzVector originalJetP4 = originalJet.p4();
0145 if (verbosity_) {
0146 std::cout << "originalJet: Pt = " << originalJetP4.pt() << ", eta = " << originalJetP4.eta()
0147 << ", phi = " << originalJetP4.phi() << std::endl;
0148 }
0149
0150 double shift = 0.;
0151 if (jecUncertaintyValue_) {
0152 shift = *jecUncertaintyValue_;
0153 } else {
0154 jecUncertainty_->setJetEta(originalJetP4.eta());
0155 jecUncertainty_->setJetPt(originalJetP4.pt());
0156
0157 shift = jecUncertainty_->getUncertainty(true);
0158 }
0159 if (verbosity_) {
0160 std::cout << "shift = " << shift << std::endl;
0161 }
0162
0163 if (evt.isRealData() && addResidualJES_) {
0164 reco::Candidate::LorentzVector rawJetP4 = pat::RawJetExtractorT<T>{}(originalJet);
0165 if (rawJetP4.E() > 1.e-1) {
0166 reco::Candidate::LorentzVector corrJetP4upToL3 =
0167 jetCorrExtractor_(originalJet, jetCorrUpToL3.product(), jetCorrEtaMax_, &rawJetP4);
0168 reco::Candidate::LorentzVector corrJetP4upToL3Res =
0169 jetCorrExtractor_(originalJet, jetCorrUpToL3Res.product(), jetCorrEtaMax_, &rawJetP4);
0170 if (corrJetP4upToL3.E() > 1.e-1 && corrJetP4upToL3Res.E() > 1.e-1) {
0171 double residualJES = (corrJetP4upToL3Res.E() / corrJetP4upToL3.E()) - 1.;
0172 shift = sqrt(shift * shift + residualJES * residualJES);
0173 }
0174 }
0175 }
0176
0177 shift *= shiftBy_;
0178 if (verbosity_) {
0179 std::cout << "shift*shiftBy = " << shift << std::endl;
0180 }
0181
0182 shiftedJets.emplace_back(originalJet);
0183 shiftedJets.back().setP4((1. + shift) * originalJetP4);
0184 if (verbosity_) {
0185 auto const& shiftedJet = shiftedJets.back();
0186 std::cout << "shiftedJet: Pt = " << shiftedJet.pt() << ", eta = " << shiftedJet.eta()
0187 << ", phi = " << shiftedJet.phi() << std::endl;
0188 }
0189 }
0190
0191 evt.emplace(putToken_, std::move(shiftedJets));
0192 }
0193
0194 std::string moduleLabel_;
0195
0196 edm::InputTag src_;
0197 edm::EDGetTokenT<JetCollection> srcToken_;
0198 edm::EDPutTokenT<JetCollection> putToken_;
0199
0200 edm::FileInPath jetCorrInputFileName_;
0201 std::string jetCorrPayloadName_;
0202 edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> jetCorrPayloadToken_;
0203 std::string jetCorrUncertaintyTag_;
0204 std::optional<JetCorrectionUncertainty> jecUncertainty_;
0205 unsigned long long recordIdentifier_ = 0;
0206
0207 bool addResidualJES_;
0208 edm::InputTag jetCorrLabelUpToL3_;
0209 edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3_;
0210 edm::InputTag jetCorrLabelUpToL3Res_;
0211 edm::EDGetTokenT<reco::JetCorrector> jetCorrTokenUpToL3Res_;
0212 double jetCorrEtaMax_;
0213
0214
0215
0216
0217 Textractor jetCorrExtractor_;
0218
0219 std::optional<double> jecUncertaintyValue_;
0220
0221 double shiftBy_;
0222
0223 int verbosity_;
0224 };
0225
0226 #endif