FFTJetProducers/plugins/FFTJetPileupEstimator.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoJets/FFTJetProducers
0004 // Class:      FFTJetPileupEstimator
0005 //
0006 /**\class FFTJetPileupEstimator FFTJetPileupEstimator.cc RecoJets/FFTJetProducers/plugins/FFTJetPileupEstimator.cc
0007
0008  Description: applies calibration curve and estimates the actual pileup
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Wed Apr 20 13:52:23 CDT 2011
0016 //
0017 //
0018
0019 #include <cmath>
0020
0021 // Framework include files
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026
0027 // Data formats
0028 #include "DataFormats/Common/interface/View.h"
0029 #include "DataFormats/Common/interface/Handle.h"
0030 // #include "DataFormats/Histograms/interface/MEtoEDMFormat.h"
0031 #include "DataFormats/JetReco/interface/FFTJetPileupSummary.h"
0032 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0033
0034 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0035
0036 // Loader for the lookup tables
0037 #include "JetMETCorrections/FFTJetObjects/interface/FFTJetLookupTableSequenceLoader.h"
0038
0039 #define init_param(type, varname) varname(ps.getParameter<type>(#varname))
0040
0041 using namespace fftjetcms;
0042
0043 //
0044 // class declaration
0045 //
0046 class FFTJetPileupEstimator : public edm::stream::EDProducer<> {
0047 public:
0048   explicit FFTJetPileupEstimator(const edm::ParameterSet&);
0049   FFTJetPileupEstimator() = delete;
0050   FFTJetPileupEstimator(const FFTJetPileupEstimator&) = delete;
0051   FFTJetPileupEstimator& operator=(const FFTJetPileupEstimator&) = delete;
0052   ~FFTJetPileupEstimator() override;
0053
0054 protected:
0055   // methods
0056   void produce(edm::Event&, const edm::EventSetup&) override;
0057
0058 private:
0059   std::unique_ptr<reco::FFTJetPileupSummary> calibrateFromConfig(double uncalibrated) const;
0060
0061   std::unique_ptr<reco::FFTJetPileupSummary> calibrateFromDB(double uncalibrated, const edm::EventSetup& iSetup) const;
0062
0063   template <class Ptr>
0064   inline void checkConfig(const Ptr& ptr, const char* message) {
0065     if (ptr.get() == nullptr)
0066       throw cms::Exception("FFTJetBadConfig") << message << std::endl;
0067   }
0068
0069   edm::InputTag inputLabel;
0070   edm::EDGetTokenT<reco::DiscretizedEnergyFlow> inputToken;
0071
0072   std::string outputLabel;
0073   double cdfvalue;
0074   double ptToDensityFactor;
0075   unsigned filterNumber;
0076   std::vector<double> uncertaintyZones;
0077   std::unique_ptr<fftjet::Functor1<double, double> > calibrationCurve;
0078   std::unique_ptr<fftjet::Functor1<double, double> > uncertaintyCurve;
0079
0080   // Alternative method to calibrate the pileup.
0081   // We will fetch three lookup tables from the database:
0082   //
0083   // 1. calibration curve table
0084   //
0085   // 2. uncertainty curve table
0086   //
0087   // 3. the table that will lookup the uncertainty code
0088   //    given the uncalibrated pileup value
0089   //
0090   // It will be assumed that all these tables will have
0091   // the same record and category but different names.
0092   //
0093   std::string calibTableRecord;
0094   std::string calibTableCategory;
0095   std::string uncertaintyZonesName;
0096   std::string calibrationCurveName;
0097   std::string uncertaintyCurveName;
0098   bool loadCalibFromDB;
0099
0100   FFTJetLookupTableSequenceLoader esLoader;
0101 };
0102
0103 //
0104 // constructors and destructor
0105 //
0106 FFTJetPileupEstimator::FFTJetPileupEstimator(const edm::ParameterSet& ps)
0107     : init_param(edm::InputTag, inputLabel),
0108       init_param(std::string, outputLabel),
0109       init_param(double, cdfvalue),
0110       init_param(double, ptToDensityFactor),
0111       init_param(unsigned, filterNumber),
0112       init_param(std::vector<double>, uncertaintyZones),
0113       init_param(std::string, calibTableRecord),
0114       init_param(std::string, calibTableCategory),
0115       init_param(std::string, uncertaintyZonesName),
0116       init_param(std::string, calibrationCurveName),
0117       init_param(std::string, uncertaintyCurveName),
0118       init_param(bool, loadCalibFromDB) {
0119   calibrationCurve = fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("calibrationCurve"));
0120   checkConfig(calibrationCurve, "bad calibration curve definition");
0121
0122   uncertaintyCurve = fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("uncertaintyCurve"));
0123   checkConfig(uncertaintyCurve, "bad uncertainty curve definition");
0124
0125   inputToken = consumes<reco::DiscretizedEnergyFlow>(inputLabel);
0126
0127   produces<reco::FFTJetPileupSummary>(outputLabel);
0128
0129   esLoader.acquireToken(calibTableRecord, consumesCollector());
0130 }
0131
0132 FFTJetPileupEstimator::~FFTJetPileupEstimator() {}
0133
0134 //
0135 // member functions
0136 //
0137
0138 // ------------ method called to for each event  ------------
0139 void FFTJetPileupEstimator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0140   edm::Handle<reco::DiscretizedEnergyFlow> input;
0141   iEvent.getByToken(inputToken, input);
0142
0143   const reco::DiscretizedEnergyFlow& h(*input);
0144   const unsigned nScales = h.nEtaBins();
0145   const unsigned nCdfvalues = h.nPhiBins();
0146
0147   const unsigned fixedCdfvalueBin = static_cast<unsigned>(std::floor(cdfvalue * nCdfvalues));
0148   if (fixedCdfvalueBin >= nCdfvalues) {
0149     throw cms::Exception("FFTJetBadConfig") << "Bad cdf value" << std::endl;
0150   }
0151   if (filterNumber >= nScales) {
0152     throw cms::Exception("FFTJetBadConfig") << "Bad filter number" << std::endl;
0153   }
0154
0155   // Simple fixed-point pile-up estimate
0156   const double curve = h.data()[filterNumber * nCdfvalues + fixedCdfvalueBin];
0157
0158   std::unique_ptr<reco::FFTJetPileupSummary> summary;
0159   if (loadCalibFromDB)
0160     summary = calibrateFromDB(curve, iSetup);
0161   else
0162     summary = calibrateFromConfig(curve);
0163   iEvent.put(std::move(summary), outputLabel);
0164 }
0165
0166 std::unique_ptr<reco::FFTJetPileupSummary> FFTJetPileupEstimator::calibrateFromConfig(const double curve) const {
0167   const double pileupRho = ptToDensityFactor * (*calibrationCurve)(curve);
0168   const double rhoUncert = ptToDensityFactor * (*uncertaintyCurve)(curve);
0169
0170   // Determine the uncertainty zone of the estimate. The "curve"
0171   // has to be above or equal to uncertaintyZones[i]  but below
0172   // uncertaintyZones[i + 1] (the second condition is also satisfied
0173   // by i == uncertaintyZones.size() - 1). Of course, it is assumed
0174   // that the vector of zones is configured appropriately -- the zone
0175   // boundaries must be presented in the increasing order.
0176   int uncertaintyCode = -1;
0177   if (!uncertaintyZones.empty()) {
0178     const unsigned nZones = uncertaintyZones.size();
0179     for (unsigned i = 0; i < nZones; ++i)
0180       if (curve >= uncertaintyZones[i]) {
0181         if (i == nZones - 1U) {
0182           uncertaintyCode = i;
0183           break;
0184         } else if (curve < uncertaintyZones[i + 1]) {
0185           uncertaintyCode = i;
0186           break;
0187         }
0188       }
0189   }
0190
0191   return std::make_unique<reco::FFTJetPileupSummary>(curve, pileupRho, rhoUncert, uncertaintyCode);
0192 }
0193
0194 std::unique_ptr<reco::FFTJetPileupSummary> FFTJetPileupEstimator::calibrateFromDB(const double curve,
0195                                                                                   const edm::EventSetup& iSetup) const {
0196   edm::ESHandle<FFTJetLookupTableSequence> h = esLoader.load(calibTableRecord, iSetup);
0197   std::shared_ptr<npstat::StorableMultivariateFunctor> uz = (*h)[calibTableCategory][uncertaintyZonesName];
0198   std::shared_ptr<npstat::StorableMultivariateFunctor> cc = (*h)[calibTableCategory][calibrationCurveName];
0199   std::shared_ptr<npstat::StorableMultivariateFunctor> uc = (*h)[calibTableCategory][uncertaintyCurveName];
0200
0201   const double pileupRho = ptToDensityFactor * (*cc)(&curve, 1U);
0202   const double rhoUncert = ptToDensityFactor * (*uc)(&curve, 1U);
0203   const int uncertaintyCode = round((*uz)(&curve, 1U));
0204
0205   return std::make_unique<reco::FFTJetPileupSummary>(curve, pileupRho, rhoUncert, uncertaintyCode);
0206 }
0207
0208 //define this as a plug-in
0209 DEFINE_FWK_MODULE(FFTJetPileupEstimator);