File indexing completed on 2023-03-17 11:18:18
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, esLoaderCalo); \
0053 break; \
0054 case PFJET: \
0055 method<reco::PFJet>(arg1, arg2, esLoaderPF); \
0056 break; \
0057 case GENJET: \
0058 method<reco::GenJet>(arg1, arg2, esLoaderGen); \
0059 break; \
0060 case TRACKJET: \
0061 method<reco::TrackJet>(arg1, arg2, esLoaderTrack); \
0062 break; \
0063 case BASICJET: \
0064 method<reco::BasicJet>(arg1, arg2, esLoaderBasic); \
0065 break; \
0066 case JPTJET: \
0067 method<reco::JPTJet>(arg1, arg2, esLoaderJPT); \
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,
0097 const std::string& tag,
0098 std::unique_ptr<typename FFTJetCorrectorSequenceTypemap<reco::FFTAnyJet<Jet>>::loader>& ptr);
0099
0100 template <typename Jet>
0101 void applyCorrections(edm::Event& iEvent,
0102 const edm::EventSetup& iSetup,
0103 std::unique_ptr<typename FFTJetCorrectorSequenceTypemap<reco::FFTAnyJet<Jet>>::loader>& ptr);
0104
0105 template <typename Jet>
0106 void performPileupSubtraction(Jet&);
0107
0108
0109 const edm::InputTag inputLabel;
0110
0111
0112 const std::string outputLabel;
0113
0114
0115 const JetType jetType;
0116
0117
0118 const std::vector<std::string> records;
0119
0120
0121 const bool writeUncertainties;
0122
0123
0124 const bool subtractPileup;
0125 const bool subtractPileupAs4Vec;
0126
0127
0128 const bool verbose;
0129
0130
0131 std::vector<int> sequenceMasks;
0132
0133
0134 unsigned long eventCount;
0135
0136
0137 edm::EDGetTokenT<std::vector<reco::FFTAnyJet<reco::Jet>>> input_jets_token_;
0138
0139
0140 std::unique_ptr<FFTBasicJetCorrectorSequenceLoader> esLoaderBasic;
0141 std::unique_ptr<FFTCaloJetCorrectorSequenceLoader> esLoaderCalo;
0142 std::unique_ptr<FFTGenJetCorrectorSequenceLoader> esLoaderGen;
0143 std::unique_ptr<FFTPFJetCorrectorSequenceLoader> esLoaderPF;
0144 std::unique_ptr<FFTTrackJetCorrectorSequenceLoader> esLoaderTrack;
0145 std::unique_ptr<FFTJPTJetCorrectorSequenceLoader> esLoaderJPT;
0146 };
0147
0148 template <typename Jet>
0149 void FFTJetCorrectionProducer::makeProduces(
0150 const std::string& alias,
0151 const std::string& tag,
0152 std::unique_ptr<typename FFTJetCorrectorSequenceTypemap<reco::FFTAnyJet<Jet>>::loader>& ptr) {
0153 typedef reco::FFTAnyJet<Jet> MyJet;
0154 typedef typename FFTJetCorrectorSequenceTypemap<MyJet>::loader Loader;
0155
0156 produces<std::vector<MyJet>>(tag).setBranchAlias(alias);
0157
0158 ptr = std::make_unique<Loader>();
0159 const unsigned nRecords = records.size();
0160 for (unsigned irec = 0; irec < nRecords; ++irec)
0161 ptr->acquireToken(records[irec], consumesCollector());
0162 }
0163
0164 template <typename Jet>
0165 void FFTJetCorrectionProducer::performPileupSubtraction(Jet& jet) {
0166 using reco::FFTJet;
0167
0168 FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(jet.getFFTSpecific()));
0169 const math::XYZTLorentzVector& new4vec = adjustForPileup(fftJet.f_vec(), fftJet.f_pileup(), subtractPileupAs4Vec);
0170 fftJet.setFourVec(new4vec);
0171 int status = fftJet.f_status();
0172 if (subtractPileupAs4Vec)
0173 status |= PILEUP_SUBTRACTION_MASK_4VEC;
0174 else
0175 status |= PILEUP_SUBTRACTION_MASK_PT;
0176 fftJet.setStatus(status);
0177 jet.setP4(new4vec);
0178 }
0179
0180 template <typename Jet>
0181 void FFTJetCorrectionProducer::applyCorrections(
0182 edm::Event& iEvent,
0183 const edm::EventSetup& iSetup,
0184 std::unique_ptr<typename FFTJetCorrectorSequenceTypemap<reco::FFTAnyJet<Jet>>::loader>& loader) {
0185 using reco::FFTJet;
0186 typedef reco::FFTAnyJet<Jet> MyJet;
0187 typedef std::vector<MyJet> MyCollection;
0188 typedef typename FFTJetCorrectorSequenceTypemap<MyJet>::loader Loader;
0189 typedef typename Loader::data_type CorrectorSequence;
0190 typedef typename CorrectorSequence::result_type CorrectorResult;
0191
0192
0193 const unsigned nRecords = records.size();
0194 std::vector<edm::ESHandle<CorrectorSequence>> handles(nRecords);
0195 for (unsigned irec = 0; irec < nRecords; ++irec)
0196 handles[irec] = loader->load(records[irec], iSetup);
0197
0198
0199
0200 sequenceMasks.clear();
0201 sequenceMasks.reserve(nRecords);
0202
0203 int totalMask = 0;
0204 for (unsigned irec = 0; irec < nRecords; ++irec) {
0205 int levelMask = 0;
0206 const unsigned nLevels = handles[irec]->nLevels();
0207 for (unsigned i = 0; i < nLevels; ++i) {
0208 const unsigned lev = (*handles[irec])[i].level();
0209
0210
0211
0212 if (lev) {
0213 const int mask = (1 << lev);
0214 if (totalMask & mask)
0215 throw cms::Exception("FFTJetBadConfig")
0216 << "Error in FFTJetCorrectionProducer::applyCorrections:"
0217 << " jet correction at level " << lev << " is applied more than once\n";
0218 totalMask |= mask;
0219 levelMask |= mask;
0220 }
0221 }
0222 sequenceMasks.push_back(levelMask << 12);
0223 }
0224 totalMask = (totalMask << 12);
0225
0226
0227 const bool isMC = !iEvent.isRealData();
0228
0229
0230 edm::Handle<MyCollection> jets;
0231 iEvent.getByToken(input_jets_token_, jets);
0232
0233
0234 const unsigned nJets = jets->size();
0235 auto coll = std::make_unique<MyCollection>();
0236 coll->reserve(nJets);
0237
0238
0239 bool sorted = true;
0240 double previousPt = DBL_MAX;
0241 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
0242 const MyJet& j((*jets)[ijet]);
0243
0244
0245 const int initialStatus = j.getFFTSpecific().f_status();
0246 if (initialStatus & totalMask)
0247 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0248 << "this jet collection is already corrected for some or all "
0249 << "of the specified levels\n";
0250
0251 MyJet corJ(j);
0252
0253 if (verbose) {
0254 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0255 std::cout << "++++ Evt " << eventCount << " jet " << ijet << ": pt = " << corJ.pt()
0256 << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale() << ", s = 0x" << std::hex
0257 << fj.f_status() << std::dec << std::endl;
0258 }
0259
0260
0261
0262
0263
0264 if (subtractPileup) {
0265 if (initialStatus & PILEUP_SUBTRACTION_MASK_ANY)
0266 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0267 << "this jet collection is already pileup-subtracted\n";
0268 if (!(initialStatus & PILEUP_CALCULATION_MASK))
0269 throw cms::Exception("FFTJetBadConfig") << "Error in FFTJetCorrectionProducer::applyCorrections: "
0270 << "pileup was not calculated for this jet collection\n";
0271 performPileupSubtraction(corJ);
0272
0273 if (verbose) {
0274 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0275 std::cout << " Pileup subtract"
0276 << ": pt = " << corJ.pt() << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale()
0277 << ", s = 0x" << std::hex << fj.f_status() << std::dec << std::endl;
0278 }
0279 }
0280
0281
0282 double sigmaSquared = 0.0;
0283 for (unsigned irec = 0; irec < nRecords; ++irec) {
0284 const CorrectorResult& corr = handles[irec]->correct(corJ, isMC);
0285
0286
0287 FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(corJ.getFFTSpecific()));
0288 corJ.setP4(corr.vec());
0289 fftJet.setFourVec(corr.vec());
0290
0291
0292 fftJet.setStatus(fftJet.f_status() | sequenceMasks[irec]);
0293
0294
0295 const double s = corr.sigma();
0296 sigmaSquared += s * s;
0297 }
0298
0299
0300
0301
0302
0303
0304 if (writeUncertainties)
0305 corJ.setPileup(sqrt(sigmaSquared));
0306
0307 coll->push_back(corJ);
0308
0309
0310 const double pt = corJ.pt();
0311 if (pt > previousPt)
0312 sorted = false;
0313 previousPt = pt;
0314
0315 if (verbose) {
0316 const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
0317 std::cout << " Fully corrected"
0318 << ": pt = " << corJ.pt() << ", eta = " << fj.f_vec().eta() << ", R = " << fj.f_recoScale()
0319 << ", s = 0x" << std::hex << fj.f_status() << std::dec << std::endl;
0320 }
0321 }
0322
0323 if (!sorted)
0324 std::sort(coll->begin(), coll->end(), LocalSortByPt());
0325
0326
0327 if (writeUncertainties) {
0328 auto unc = std::make_unique<std::vector<float>>();
0329 unc->reserve(nJets);
0330 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
0331 MyJet& j((*coll)[ijet]);
0332 unc->push_back(j.pileup());
0333 j.setPileup(0.f);
0334 }
0335 iEvent.put(std::move(unc), outputLabel);
0336 }
0337
0338 iEvent.put(std::move(coll), outputLabel);
0339 ++eventCount;
0340 }
0341
0342
0343
0344
0345 FFTJetCorrectionProducer::FFTJetCorrectionProducer(const edm::ParameterSet& ps)
0346 : inputLabel(ps.getParameter<edm::InputTag>("src")),
0347 outputLabel(ps.getParameter<std::string>("outputLabel")),
0348 jetType(parseJetType(ps.getParameter<std::string>("jetType"))),
0349 records(ps.getParameter<std::vector<std::string>>("records")),
0350 writeUncertainties(ps.getParameter<bool>("writeUncertainties")),
0351 subtractPileup(ps.getParameter<bool>("subtractPileup")),
0352 subtractPileupAs4Vec(ps.getParameter<bool>("subtractPileupAs4Vec")),
0353 verbose(ps.getUntrackedParameter<bool>("verbose", false)),
0354 eventCount(0UL) {
0355 const std::string alias(ps.getUntrackedParameter<std::string>("alias", outputLabel));
0356 jet_type_switch(makeProduces, alias, outputLabel);
0357
0358 if (writeUncertainties)
0359 produces<std::vector<float>>(outputLabel).setBranchAlias(alias);
0360
0361 input_jets_token_ = consumes<std::vector<reco::FFTAnyJet<reco::Jet>>>(inputLabel);
0362 }
0363
0364 FFTJetCorrectionProducer::~FFTJetCorrectionProducer() {}
0365
0366
0367 void FFTJetCorrectionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0368 jet_type_switch(applyCorrections, iEvent, iSetup);
0369 }
0370
0371
0372 DEFINE_FWK_MODULE(FFTJetCorrectionProducer);