FFTJetProducers/plugins/FFTJetCorrectionProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoJets/FFTJetProducers
0004 // Class:      FFTJetCorrectionProducer
0005 //
0006 /**\class FFTJetCorrectionProducer FFTJetCorrectionProducer.cc RecoJets/FFTJetProducers/plugins/FFTJetCorrectionProducer.cc
0007
0008  Description: producer for correcting jets created by FFTJetProducer
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Mon Aug  6 11:03:38 CDT 2012
0016 //
0017 //
0018
0019 // system include files
0020 #include <iostream>
0021 #include <memory>
0022 #include <cfloat>
0023 #include <cmath>
0024
0025 // user include files
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 // A generic switch statement based on jet type
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 }  // namespace
0083
0084 //
0085 // class declaration
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   // Label for the input collection
0105   const edm::InputTag inputLabel;
0106
0107   // Label for the output objects
0108   const std::string outputLabel;
0109
0110   // Jet type to process
0111   const JetType jetType;
0112
0113   // The names of the jet correction records
0114   const std::vector<std::string> records;
0115
0116   // Are we going to create the uncertainty collection?
0117   const bool writeUncertainties;
0118
0119   // What to do about pileup subtraction
0120   const bool subtractPileup;
0121   const bool subtractPileupAs4Vec;
0122
0123   // Print some info about jets
0124   const bool verbose;
0125
0126   // Space for level masks
0127   std::vector<int> sequenceMasks;
0128
0129   // Event counter
0130   unsigned long eventCount;
0131
0132   // tokens for data access
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   // Load the jet corrector sequences
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   // Figure out which correction levels we are applying
0173   // and create masks which will indicate this
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       // Not tracking "level 0" corrections in the status word.
0185       // Level 0 is basically reserved for uncertainty calculations.
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   // Is this data or MC?
0201   const bool isMC = !iEvent.isRealData();
0202
0203   // Load the jet collection
0204   edm::Handle<MyCollection> jets;
0205   iEvent.getByToken(input_jets_token_, jets);
0206
0207   // Create the output collection
0208   const unsigned nJets = jets->size();
0209   auto coll = std::make_unique<MyCollection>();
0210   coll->reserve(nJets);
0211
0212   // Cycle over jets and apply the corrector sequences
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     // Check that this jet has not been corrected yet
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     // Check if we need to subtract pileup first.
0235     // Pileup subtraction is not part of the corrector sequence
0236     // itself because 4-vector subtraction does not map well
0237     // into multiplication of 4-vectors by a scale factor.
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     // Apply all jet correction sequences
0256     double sigmaSquared = 0.0;
0257     for (unsigned irec = 0; irec < nRecords; ++irec) {
0258       const CorrectorResult& corr = handles[irec]->correct(corJ, isMC);
0259
0260       // Update the 4-vector
0261       FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(corJ.getFFTSpecific()));
0262       corJ.setP4(corr.vec());
0263       fftJet.setFourVec(corr.vec());
0264
0265       // Update the jet correction status
0266       fftJet.setStatus(fftJet.f_status() | sequenceMasks[irec]);
0267
0268       // Update the (systematic) uncertainty
0269       const double s = corr.sigma();
0270       sigmaSquared += s * s;
0271     }
0272
0273     // There is no place for uncertainty in the jet structure.
0274     // However, there is the unused pileup field (FFTJet maintains
0275     // the pileup separately as a 4-vector). Use this unused field
0276     // to store the uncertainty. This hack is needed because
0277     // subsequent sequence sorting by Pt can change the jet ordering.
0278     if (writeUncertainties)
0279       corJ.setPileup(sqrt(sigmaSquared));
0280
0281     coll->push_back(corJ);
0282
0283     // Check whether the sequence remains sorted by pt
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   // Create the uncertainty sequence
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 // constructors and destructor
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 // ------------ method called to produce the data  ------------
0341 void FFTJetCorrectionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0342   jet_type_switch(applyCorrections, iEvent, iSetup);
0343 }
0344
0345 //define this as a plug-in
0346 DEFINE_FWK_MODULE(FFTJetCorrectionProducer);