File indexing completed on 2022-05-21 03:43:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <iostream>
0021 #include <memory>
0022 #include <cfloat>
0023 #include <cmath>
0024
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/stream/EDProducer.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033
0034 #include "RecoJets/FFTJetAlgorithms/interface/adjustForPileup.h"
0035 #include "RecoJets/FFTJetProducers/interface/JetType.h"
0036 #include "JetMETCorrections/FFTJetObjects/interface/FFTJetCorrectorSequenceTypemap.h"
0037
0038 #define PILEUP_CALCULATION_MASK 0x200
0039 #define PILEUP_SUBTRACTION_MASK_4VEC 0x400
0040 #define PILEUP_SUBTRACTION_MASK_PT 0x800
0041 #define PILEUP_SUBTRACTION_MASK_ANY (PILEUP_SUBTRACTION_MASK_4VEC | PILEUP_SUBTRACTION_MASK_PT)
0042
0043 using namespace fftjetcms;
0044
0045
0046
0047
0048 #define jet_type_switch(method, arg1, arg2) \
0049 do { \
0050 switch (jetType) { \
0051 case CALOJET: \
0052 method<reco::CaloJet>(arg1, arg2); \
0053 break; \
0054 case PFJET: \
0055 method<reco::PFJet>(arg1, arg2); \
0056 break; \
0057 case GENJET: \
0058 method<reco::GenJet>(arg1, arg2); \
0059 break; \
0060 case TRACKJET: \
0061 method<reco::TrackJet>(arg1, arg2); \
0062 break; \
0063 case BASICJET: \
0064 method<reco::BasicJet>(arg1, arg2); \
0065 break; \
0066 case JPTJET: \
0067 method<reco::JPTJet>(arg1, arg2); \
0068 break; \
0069 default: \
0070 assert(!"ERROR in FFTJetCorrectionProducer : invalid jet type."\
0071 " This is a bug. Please report."); \
0072 } \
0073 } while (0);
0074
0075 namespace {
0076 struct LocalSortByPt {
0077 template <class Jet>
0078 inline bool operator()(const Jet& l, const Jet& r) const {
0079 return l.pt() > r.pt();
0080 }
0081 };
0082 }
0083
0084
0085
0086
0087 class FFTJetCorrectionProducer : public edm::stream::EDProducer<> {
0088 public:
0089 explicit FFTJetCorrectionProducer(const edm::ParameterSet&);
0090 ~FFTJetCorrectionProducer() override;
0091
0092 private:
0093 void produce(edm::Event&, const edm::EventSetup&) override;
0094
0095 template <typename Jet>
0096 void makeProduces(const std::string& alias, const std::string& tag);
0097
0098 template <typename Jet>
0099 void applyCorrections(edm::Event& iEvent, const edm::EventSetup& iSetup);
0100
0101 template <typename Jet>
0102 void performPileupSubtraction(Jet&);
0103
0104
0105 const edm::InputTag inputLabel;
0106
0107
0108 const std::string outputLabel;
0109
0110
0111 const JetType jetType;
0112
0113
0114 const std::vector<std::string> records;
0115
0116
0117 const bool writeUncertainties;
0118
0119
0120 const bool subtractPileup;
0121 const bool subtractPileupAs4Vec;
0122
0123
0124 const bool verbose;
0125
0126
0127 std::vector<int> sequenceMasks;
0128
0129
0130 unsigned long eventCount;
0131
0132
0133 edm::EDGetTokenT<std::vector<reco::FFTAnyJet<reco::Jet>>> input_jets_token_;
0134 };
0135
0136 template <typename T>
0137 void FFTJetCorrectionProducer::makeProduces(const std::string& alias, const std::string& tag) {
0138 produces<std::vector<reco::FFTAnyJet<T>>>(tag).setBranchAlias(alias);
0139 }
0140
0141 template <typename Jet>
0142 void FFTJetCorrectionProducer::performPileupSubtraction(Jet& jet) {
0143 using reco::FFTJet;
0144
0145 FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(jet.getFFTSpecific()));
0146 const math::XYZTLorentzVector& new4vec = adjustForPileup(fftJet.f_vec(), fftJet.f_pileup(), subtractPileupAs4Vec);
0147 fftJet.setFourVec(new4vec);
0148 int status = fftJet.f_status();
0149 if (subtractPileupAs4Vec)
0150 status |= PILEUP_SUBTRACTION_MASK_4VEC;
0151 else
0152 status |= PILEUP_SUBTRACTION_MASK_PT;
0153 fftJet.setStatus(status);
0154 jet.setP4(new4vec);
0155 }
0156
0157 template <typename Jet>
0158 void FFTJetCorrectionProducer::applyCorrections(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0159 using reco::FFTJet;
0160 typedef reco::FFTAnyJet<Jet> MyJet;
0161 typedef std::vector<MyJet> MyCollection;
0162 typedef typename FFTJetCorrectorSequenceTypemap<MyJet>::loader Loader;
0163 typedef typename Loader::data_type CorrectorSequence;
0164 typedef typename CorrectorSequence::result_type CorrectorResult;
0165
0166
0167 const unsigned nRecords = records.size();
0168 std::vector<edm::ESHandle<CorrectorSequence>> handles(nRecords);
0169 for (unsigned irec = 0; irec < nRecords; ++irec)
0170 Loader::instance().load(iSetup, records[irec], handles[irec]);
0171
0172
0173
0174 sequenceMasks.clear();
0175 sequenceMasks.reserve(nRecords);
0176
0177 int totalMask = 0;
0178 for (unsigned irec = 0; irec < nRecords; ++irec) {
0179 int levelMask = 0;
0180 const unsigned nLevels = handles[irec]->nLevels();
0181 for (unsigned i = 0; i < nLevels; ++i) {
0182 const unsigned lev = (*handles[irec])[i].level();
0183
0184
0185
0186 if (lev) {
0187 const int mask = (1 << lev);
0188 if (totalMask & mask)
0189 throw cms::Exception("FFTJetBadConfig")
0190 << "Error in FFTJetCorrectionProducer::applyCorrections:"
0191 << " jet correction at level " << lev << " is applied more than once\n";
0192 totalMask |= mask;
0193 levelMask |= mask;
0194 }
0195 }
0196 sequenceMasks.push_back(levelMask << 12);
0197 }
0198 totalMask = (totalMask << 12);
0199
0200
0201 const bool isMC = !iEvent.isRealData();
0202
0203
0204 edm::Handle<MyCollection> jets;
0205 iEvent.getByToken(input_jets_token_, jets);
0206
0207
0208 const unsigned nJets = jets->size();
0209 auto coll = std::make_unique<MyCollection>();
0210 coll->reserve(nJets);
0211
0212
0213 bool sorted = true;
0214 double previousPt = DBL_MAX;
0215 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
0216 const MyJet& j((*jets)[ijet]);
0217
0218
0219 const int initialStatus = j.getFFTSpecific().f_status();
0220 if (initialStatus & totalMask)
0221 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0222 << "this jet collection is already corrected for some or all "
0223 << "of the specified levels\n";
0224
0225 MyJet corJ(j);
0226
0227 if (verbose) {
0228 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0229 std::cout << "++++ Evt " << eventCount << " jet " << ijet << ": pt = " << corJ.pt()
0230 << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale() << ", s = 0x" << std::hex
0231 << fj.f_status() << std::dec << std::endl;
0232 }
0233
0234
0235
0236
0237
0238 if (subtractPileup) {
0239 if (initialStatus & PILEUP_SUBTRACTION_MASK_ANY)
0240 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0241 << "this jet collection is already pileup-subtracted\n";
0242 if (!(initialStatus & PILEUP_CALCULATION_MASK))
0243 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0244 << "pileup was not calculated for this jet collection\n";
0245 performPileupSubtraction(corJ);
0246
0247 if (verbose) {
0248 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0249 std::cout << " Pileup subtract"
0250 << ": pt = " << corJ.pt() << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale()
0251 << ", s = 0x" << std::hex << fj.f_status() << std::dec << std::endl;
0252 }
0253 }
0254
0255
0256 double sigmaSquared = 0.0;
0257 for (unsigned irec = 0; irec < nRecords; ++irec) {
0258 const CorrectorResult& corr = handles[irec]->correct(corJ, isMC);
0259
0260
0261 FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(corJ.getFFTSpecific()));
0262 corJ.setP4(corr.vec());
0263 fftJet.setFourVec(corr.vec());
0264
0265
0266 fftJet.setStatus(fftJet.f_status() | sequenceMasks[irec]);
0267
0268
0269 const double s = corr.sigma();
0270 sigmaSquared += s * s;
0271 }
0272
0273
0274
0275
0276
0277
0278 if (writeUncertainties)
0279 corJ.setPileup(sqrt(sigmaSquared));
0280
0281 coll->push_back(corJ);
0282
0283
0284 const double pt = corJ.pt();
0285 if (pt > previousPt)
0286 sorted = false;
0287 previousPt = pt;
0288
0289 if (verbose) {
0290 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0291 std::cout << " Fully corrected"
0292 << ": pt = " << corJ.pt() << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale()
0293 << ", s = 0x" << std::hex << fj.f_status() << std::dec << std::endl;
0294 }
0295 }
0296
0297 if (!sorted)
0298 std::sort(coll->begin(), coll->end(), LocalSortByPt());
0299
0300
0301 if (writeUncertainties) {
0302 auto unc = std::make_unique<std::vector<float>>();
0303 unc->reserve(nJets);
0304 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
0305 MyJet& j((*coll)[ijet]);
0306 unc->push_back(j.pileup());
0307 j.setPileup(0.f);
0308 }
0309 iEvent.put(std::move(unc), outputLabel);
0310 }
0311
0312 iEvent.put(std::move(coll), outputLabel);
0313 ++eventCount;
0314 }
0315
0316
0317
0318
0319 FFTJetCorrectionProducer::FFTJetCorrectionProducer(const edm::ParameterSet& ps)
0320 : inputLabel(ps.getParameter<edm::InputTag>("src")),
0321 outputLabel(ps.getParameter<std::string>("outputLabel")),
0322 jetType(parseJetType(ps.getParameter<std::string>("jetType"))),
0323 records(ps.getParameter<std::vector<std::string>>("records")),
0324 writeUncertainties(ps.getParameter<bool>("writeUncertainties")),
0325 subtractPileup(ps.getParameter<bool>("subtractPileup")),
0326 subtractPileupAs4Vec(ps.getParameter<bool>("subtractPileupAs4Vec")),
0327 verbose(ps.getUntrackedParameter<bool>("verbose", false)),
0328 eventCount(0UL) {
0329 const std::string alias(ps.getUntrackedParameter<std::string>("alias", outputLabel));
0330 jet_type_switch(makeProduces, alias, outputLabel);
0331
0332 if (writeUncertainties)
0333 produces<std::vector<float>>(outputLabel).setBranchAlias(alias);
0334
0335 input_jets_token_ = consumes<std::vector<reco::FFTAnyJet<reco::Jet>>>(inputLabel);
0336 }
0337
0338 FFTJetCorrectionProducer::~FFTJetCorrectionProducer() {}
0339
0340
0341 void FFTJetCorrectionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0342 jet_type_switch(applyCorrections, iEvent, iSetup);
0343 }
0344
0345
0346 DEFINE_FWK_MODULE(FFTJetCorrectionProducer);