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, 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 }  // 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,
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   // Label for the input collection
0109   const edm::InputTag inputLabel;
0110
0111   // Label for the output objects
0112   const std::string outputLabel;
0113
0114   // Jet type to process
0115   const JetType jetType;
0116
0117   // The names of the jet correction records
0118   const std::vector<std::string> records;
0119
0120   // Are we going to create the uncertainty collection?
0121   const bool writeUncertainties;
0122
0123   // What to do about pileup subtraction
0124   const bool subtractPileup;
0125   const bool subtractPileupAs4Vec;
0126
0127   // Print some info about jets
0128   const bool verbose;
0129
0130   // Space for level masks
0131   std::vector<int> sequenceMasks;
0132
0133   // Event counter
0134   unsigned long eventCount;
0135
0136   // tokens for data access
0137   edm::EDGetTokenT<std::vector<reco::FFTAnyJet<reco::Jet>>> input_jets_token_;
0138
0139   // Tools for accessing event setup
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   // Load the jet corrector sequences
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   // Figure out which correction levels we are applying
0199   // and create masks which will indicate this
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       // Not tracking "level 0" corrections in the status word.
0211       // Level 0 is basically reserved for uncertainty calculations.
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   // Is this data or MC?
0227   const bool isMC = !iEvent.isRealData();
0228
0229   // Load the jet collection
0230   edm::Handle<MyCollection> jets;
0231   iEvent.getByToken(input_jets_token_, jets);
0232
0233   // Create the output collection
0234   const unsigned nJets = jets->size();
0235   auto coll = std::make_unique<MyCollection>();
0236   coll->reserve(nJets);
0237
0238   // Cycle over jets and apply the corrector sequences
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     // Check that this jet has not been corrected yet
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     // Check if we need to subtract pileup first.
0261     // Pileup subtraction is not part of the corrector sequence
0262     // itself because 4-vector subtraction does not map well
0263     // into multiplication of 4-vectors by a scale factor.
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     // Apply all jet correction sequences
0282     double sigmaSquared = 0.0;
0283     for (unsigned irec = 0; irec < nRecords; ++irec) {
0284       const CorrectorResult& corr = handles[irec]->correct(corJ, isMC);
0285
0286       // Update the 4-vector
0287       FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(corJ.getFFTSpecific()));
0288       corJ.setP4(corr.vec());
0289       fftJet.setFourVec(corr.vec());
0290
0291       // Update the jet correction status
0292       fftJet.setStatus(fftJet.f_status() | sequenceMasks[irec]);
0293
0294       // Update the (systematic) uncertainty
0295       const double s = corr.sigma();
0296       sigmaSquared += s * s;
0297     }
0298
0299     // There is no place for uncertainty in the jet structure.
0300     // However, there is the unused pileup field (FFTJet maintains
0301     // the pileup separately as a 4-vector). Use this unused field
0302     // to store the uncertainty. This hack is needed because
0303     // subsequent sequence sorting by Pt can change the jet ordering.
0304     if (writeUncertainties)
0305       corJ.setPileup(sqrt(sigmaSquared));
0306
0307     coll->push_back(corJ);
0308
0309     // Check whether the sequence remains sorted by pt
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   // Create the uncertainty sequence
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 // constructors and destructor
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 // ------------ method called to produce the data  ------------
0367 void FFTJetCorrectionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0368   jet_type_switch(applyCorrections, iEvent, iSetup);
0369 }
0370
0371 //define this as a plug-in
0372 DEFINE_FWK_MODULE(FFTJetCorrectionProducer);