FFTJetProducers/plugins/FFTJetPileupProcessor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoJets/FFTJetProducers
0004 // Class:      FFTJetPileupProcessor
0005 //
0006 /**\class FFTJetPileupProcessor FFTJetPileupProcessor.cc RecoJets/FFTJetProducers/plugins/FFTJetPileupProcessor.cc
0007
0008  Description: Runs FFTJet multiscale pileup filtering code and saves the results
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 #include <fstream>
0021 #include <memory>
0022
0023 // FFTJet headers
0024 #include "fftjet/FrequencyKernelConvolver.hh"
0025 #include "fftjet/DiscreteGauss2d.hh"
0026 #include "fftjet/EquidistantSequence.hh"
0027 #include "fftjet/interpolate.hh"
0028
0029 // framework include files
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033
0034 #include "DataFormats/Common/interface/View.h"
0035 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0036
0037 #include "RecoJets/FFTJetAlgorithms/interface/gridConverters.h"
0038
0039 // parameter parser header
0040 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0041
0042 // useful utilities collected in the second base
0043 #include "RecoJets/FFTJetProducers/interface/FFTJetInterface.h"
0044
0045 // Loader for the lookup tables
0046 #include "JetMETCorrections/FFTJetObjects/interface/FFTJetLookupTableSequenceLoader.h"
0047
0048 using namespace fftjetcms;
0049
0050 //
0051 // class declaration
0052 //
0053 class FFTJetPileupProcessor : public FFTJetInterface {
0054 public:
0055   explicit FFTJetPileupProcessor(const edm::ParameterSet&);
0056   FFTJetPileupProcessor() = delete;
0057   FFTJetPileupProcessor(const FFTJetPileupProcessor&) = delete;
0058   FFTJetPileupProcessor& operator=(const FFTJetPileupProcessor&) = delete;
0059   ~FFTJetPileupProcessor() override;
0060
0061 protected:
0062   // methods
0063   void produce(edm::Event&, const edm::EventSetup&) override;
0064
0065 private:
0066   void buildKernelConvolver(const edm::ParameterSet&);
0067   void mixExtraGrid();
0068   void loadFlatteningFactors(const edm::EventSetup& iSetup);
0069
0070   // The FFT engine
0071   std::unique_ptr<MyFFTEngine> engine;
0072
0073   // The pattern recognition kernel(s)
0074   std::unique_ptr<fftjet::AbsFrequencyKernel> kernel2d;
0075
0076   // The kernel convolver
0077   std::unique_ptr<fftjet::AbsConvolverBase<Real>> convolver;
0078
0079   // Storage for convolved energy flow
0080   std::unique_ptr<fftjet::Grid2d<fftjetcms::Real>> convolvedFlow;
0081
0082   // Filtering scales
0083   std::unique_ptr<fftjet::EquidistantInLogSpace> filterScales;
0084
0085   // Eta-dependent factors to use for flattening the distribution
0086   // _after_ the filtering
0087   std::vector<double> etaFlatteningFactors;
0088
0089   // Number of percentile points to use
0090   unsigned nPercentiles;
0091
0092   // Bin range. Both values of 0 means use the whole range.
0093   unsigned convolverMinBin;
0094   unsigned convolverMaxBin;
0095
0096   // Density conversion factor
0097   double pileupEtaPhiArea;
0098
0099   // Variable related to mixing additional grids
0100   std::vector<std::string> externalGridFiles;
0101   std::ifstream gridStream;
0102   double externalGridMaxEnergy;
0103   unsigned currentFileNum;
0104
0105   // Some memory to hold the percentiles found
0106   std::vector<double> percentileData;
0107
0108   // Variables to load the flattening factors from
0109   // the calibration database (this has to be done
0110   // in sync with configuring the appropriate event
0111   // setup source)
0112   std::string flatteningTableRecord;
0113   std::string flatteningTableName;
0114   std::string flatteningTableCategory;
0115   bool loadFlatteningFactorsFromDB;
0116
0117   FFTJetLookupTableSequenceLoader esLoader;
0118 };
0119
0120 //
0121 // constructors and destructor
0122 //
0123 FFTJetPileupProcessor::FFTJetPileupProcessor(const edm::ParameterSet& ps)
0124     : FFTJetInterface(ps),
0125       etaFlatteningFactors(ps.getParameter<std::vector<double>>("etaFlatteningFactors")),
0126       nPercentiles(ps.getParameter<unsigned>("nPercentiles")),
0127       convolverMinBin(ps.getParameter<unsigned>("convolverMinBin")),
0128       convolverMaxBin(ps.getParameter<unsigned>("convolverMaxBin")),
0129       pileupEtaPhiArea(ps.getParameter<double>("pileupEtaPhiArea")),
0130       externalGridFiles(ps.getParameter<std::vector<std::string>>("externalGridFiles")),
0131       externalGridMaxEnergy(ps.getParameter<double>("externalGridMaxEnergy")),
0132       currentFileNum(externalGridFiles.size() + 1U),
0133       flatteningTableRecord(ps.getParameter<std::string>("flatteningTableRecord")),
0134       flatteningTableName(ps.getParameter<std::string>("flatteningTableName")),
0135       flatteningTableCategory(ps.getParameter<std::string>("flatteningTableCategory")),
0136       loadFlatteningFactorsFromDB(ps.getParameter<bool>("loadFlatteningFactorsFromDB")) {
0137   // Build the discretization grid
0138   energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0139   checkConfig(energyFlow, "invalid discretization grid");
0140
0141   // Copy of the grid which will be used for convolutions
0142   convolvedFlow = std::make_unique<fftjet::Grid2d<fftjetcms::Real>>(*energyFlow);
0143
0144   // Make sure the size of flattening factors is appropriate
0145   if (!etaFlatteningFactors.empty()) {
0146     if (etaFlatteningFactors.size() != convolvedFlow->nEta())
0147       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor constructor:"
0148                                                  " number of elements in the \"etaFlatteningFactors\""
0149                                                  " vector is inconsistent with the discretization grid binning"
0150                                               << std::endl;
0151   }
0152
0153   // Build the FFT engine(s), pattern recognition kernel(s),
0154   // and the kernel convolver
0155   buildKernelConvolver(ps);
0156
0157   // Build the set of pattern recognition scales
0158   const unsigned nScales = ps.getParameter<unsigned>("nScales");
0159   const double minScale = ps.getParameter<double>("minScale");
0160   const double maxScale = ps.getParameter<double>("maxScale");
0161   if (minScale <= 0.0 || maxScale < minScale || nScales == 0U)
0162     throw cms::Exception("FFTJetBadConfig") << "invalid filter scales" << std::endl;
0163
0164   filterScales = std::make_unique<fftjet::EquidistantInLogSpace>(minScale, maxScale, nScales);
0165
0166   percentileData.resize(nScales * nPercentiles);
0167
0168   produces<reco::DiscretizedEnergyFlow>(outputLabel);
0169   produces<std::pair<double, double>>(outputLabel);
0170
0171   esLoader.acquireToken(flatteningTableRecord, consumesCollector());
0172 }
0173
0174 void FFTJetPileupProcessor::buildKernelConvolver(const edm::ParameterSet& ps) {
0175   // Get the eta and phi scales for the kernel(s)
0176   double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
0177   const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
0178   if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
0179     throw cms::Exception("FFTJetBadConfig") << "invalid kernel scales" << std::endl;
0180
0181   if (convolverMinBin || convolverMaxBin)
0182     if (convolverMinBin >= convolverMaxBin || convolverMaxBin > energyFlow->nEta())
0183       throw cms::Exception("FFTJetBadConfig") << "invalid convolver bin range" << std::endl;
0184
0185   // FFT assumes that the grid extent in eta is 2*Pi. Adjust the
0186   // kernel scale in eta to compensate.
0187   kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
0188
0189   // Build the FFT engine
0190   engine = std::make_unique<MyFFTEngine>(energyFlow->nEta(), energyFlow->nPhi());
0191
0192   // 2d kernel
0193   kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
0194       new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
0195
0196   // 2d convolver
0197   convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real>>(new fftjet::FrequencyKernelConvolver<Real, Complex>(
0198       engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
0199 }
0200
0201 FFTJetPileupProcessor::~FFTJetPileupProcessor() {}
0202
0203 // ------------ method called to produce the data  ------------
0204 void FFTJetPileupProcessor::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0205   loadInputCollection(iEvent);
0206   discretizeEnergyFlow();
0207
0208   // Determine the average Et density for this event.
0209   // Needs to be done here, before mixing in another grid.
0210   const fftjet::Grid2d<Real>& g(*energyFlow);
0211   const double densityBeforeMixing = g.sum() / pileupEtaPhiArea;
0212
0213   // Mix an extra grid (if requested)
0214   double densityAfterMixing = -1.0;
0215   if (!externalGridFiles.empty()) {
0216     mixExtraGrid();
0217     densityAfterMixing = g.sum() / pileupEtaPhiArea;
0218   }
0219
0220   // Various useful variables
0221   const unsigned nScales = filterScales->size();
0222   const double* scales = &(*filterScales)[0];
0223   Real* convData = const_cast<Real*>(convolvedFlow->data());
0224   Real* sortData = convData + convolverMinBin * convolvedFlow->nPhi();
0225   const unsigned dataLen = convolverMaxBin ? (convolverMaxBin - convolverMinBin) * convolvedFlow->nPhi()
0226                                            : convolvedFlow->nEta() * convolvedFlow->nPhi();
0227
0228   // Load flattenning factors from DB (if requested)
0229   if (loadFlatteningFactorsFromDB)
0230     loadFlatteningFactors(iSetup);
0231
0232   // Go over all scales and perform the convolutions
0233   convolver->setEventData(g.data(), g.nEta(), g.nPhi());
0234   for (unsigned iscale = 0; iscale < nScales; ++iscale) {
0235     // Perform the convolution
0236     convolver->convolveWithKernel(scales[iscale], convData, convolvedFlow->nEta(), convolvedFlow->nPhi());
0237
0238     // Apply the flattening factors
0239     if (!etaFlatteningFactors.empty())
0240       convolvedFlow->scaleData(&etaFlatteningFactors[0], etaFlatteningFactors.size());
0241
0242     // Sort the convolved data
0243     std::sort(sortData, sortData + dataLen);
0244
0245     // Determine the percentile points
0246     for (unsigned iper = 0; iper < nPercentiles; ++iper) {
0247       // Map percentile 0 into point number 0,
0248       // 1 into point number dataLen-1
0249       const double q = (iper + 0.5) / nPercentiles;
0250       const double dindex = q * (dataLen - 1U);
0251       const unsigned ilow = static_cast<unsigned>(std::floor(dindex));
0252       const double percentile =
0253           fftjet::lin_interpolate_1d(ilow, ilow + 1U, sortData[ilow], sortData[ilow + 1U], dindex);
0254
0255       // Store the calculated percentile
0256       percentileData[iscale * nPercentiles + iper] = percentile;
0257     }
0258   }
0259
0260   // Convert percentile data into a more convenient storable object
0261   // and put it into the event record
0262   iEvent.put(std::make_unique<reco::DiscretizedEnergyFlow>(
0263                  &percentileData[0], "FFTJetPileupProcessor", -0.5, nScales - 0.5, 0.0, nScales, nPercentiles),
0264              outputLabel);
0265
0266   iEvent.put(std::make_unique<std::pair<double, double>>(densityBeforeMixing, densityAfterMixing), outputLabel);
0267 }
0268
0269 void FFTJetPileupProcessor::mixExtraGrid() {
0270   const unsigned nFiles = externalGridFiles.size();
0271   if (currentFileNum > nFiles) {
0272     // This is the first time this function is called
0273     currentFileNum = 0;
0274     gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
0275     if (!gridStream.is_open())
0276       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0277                                                  " failed to open external grid file "
0278                                               << externalGridFiles[currentFileNum] << std::endl;
0279   }
0280
0281   const fftjet::Grid2d<float>* g = nullptr;
0282   const unsigned maxFail = 100U;
0283   unsigned nEnergyRejected = 0;
0284
0285   while (!g) {
0286     g = fftjet::Grid2d<float>::read(gridStream);
0287
0288     // If we can't read the grid, we need to switch to another file
0289     for (unsigned ntries = 0; ntries < nFiles && g == nullptr; ++ntries) {
0290       gridStream.close();
0291       currentFileNum = (currentFileNum + 1U) % nFiles;
0292       gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
0293       if (!gridStream.is_open())
0294         throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0295                                                    " failed to open external grid file "
0296                                                 << externalGridFiles[currentFileNum] << std::endl;
0297       g = fftjet::Grid2d<float>::read(gridStream);
0298     }
0299
0300     if (g)
0301       if (g->sum() > externalGridMaxEnergy) {
0302         delete g;
0303         g = nullptr;
0304         if (++nEnergyRejected >= maxFail)
0305           throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0306                                                      " too many grids in a row ("
0307                                                   << nEnergyRejected << ") failed the maximum energy cut" << std::endl;
0308       }
0309   }
0310
0311   if (g) {
0312     add_Grid2d_data(energyFlow.get(), *g);
0313     delete g;
0314   } else {
0315     // Too bad, no useful file found
0316     throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0317                                                " no valid grid records found"
0318                                             << std::endl;
0319   }
0320 }
0321
0322 void FFTJetPileupProcessor::loadFlatteningFactors(const edm::EventSetup& iSetup) {
0323   edm::ESHandle<FFTJetLookupTableSequence> h = esLoader.load(flatteningTableRecord, iSetup);
0324   std::shared_ptr<npstat::StorableMultivariateFunctor> f = (*h)[flatteningTableCategory][flatteningTableName];
0325
0326   // Fill out the table of flattening factors as a function of eta
0327   const unsigned nEta = energyFlow->nEta();
0328   etaFlatteningFactors.clear();
0329   etaFlatteningFactors.reserve(nEta);
0330   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
0331     const double eta = energyFlow->etaBinCenter(ieta);
0332     etaFlatteningFactors.push_back((*f)(&eta, 1U));
0333   }
0334 }
0335
0336 //define this as a plug-in
0337 DEFINE_FWK_MODULE(FFTJetPileupProcessor);