Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-21 03:43:05

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 
0118 //
0119 // constructors and destructor
0120 //
0121 FFTJetPileupProcessor::FFTJetPileupProcessor(const edm::ParameterSet& ps)
0122     : FFTJetInterface(ps),
0123       etaFlatteningFactors(ps.getParameter<std::vector<double>>("etaFlatteningFactors")),
0124       nPercentiles(ps.getParameter<unsigned>("nPercentiles")),
0125       convolverMinBin(ps.getParameter<unsigned>("convolverMinBin")),
0126       convolverMaxBin(ps.getParameter<unsigned>("convolverMaxBin")),
0127       pileupEtaPhiArea(ps.getParameter<double>("pileupEtaPhiArea")),
0128       externalGridFiles(ps.getParameter<std::vector<std::string>>("externalGridFiles")),
0129       externalGridMaxEnergy(ps.getParameter<double>("externalGridMaxEnergy")),
0130       currentFileNum(externalGridFiles.size() + 1U),
0131       flatteningTableRecord(ps.getParameter<std::string>("flatteningTableRecord")),
0132       flatteningTableName(ps.getParameter<std::string>("flatteningTableName")),
0133       flatteningTableCategory(ps.getParameter<std::string>("flatteningTableCategory")),
0134       loadFlatteningFactorsFromDB(ps.getParameter<bool>("loadFlatteningFactorsFromDB")) {
0135   // Build the discretization grid
0136   energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0137   checkConfig(energyFlow, "invalid discretization grid");
0138 
0139   // Copy of the grid which will be used for convolutions
0140   convolvedFlow = std::make_unique<fftjet::Grid2d<fftjetcms::Real>>(*energyFlow);
0141 
0142   // Make sure the size of flattening factors is appropriate
0143   if (!etaFlatteningFactors.empty()) {
0144     if (etaFlatteningFactors.size() != convolvedFlow->nEta())
0145       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor constructor:"
0146                                                  " number of elements in the \"etaFlatteningFactors\""
0147                                                  " vector is inconsistent with the discretization grid binning"
0148                                               << std::endl;
0149   }
0150 
0151   // Build the FFT engine(s), pattern recognition kernel(s),
0152   // and the kernel convolver
0153   buildKernelConvolver(ps);
0154 
0155   // Build the set of pattern recognition scales
0156   const unsigned nScales = ps.getParameter<unsigned>("nScales");
0157   const double minScale = ps.getParameter<double>("minScale");
0158   const double maxScale = ps.getParameter<double>("maxScale");
0159   if (minScale <= 0.0 || maxScale < minScale || nScales == 0U)
0160     throw cms::Exception("FFTJetBadConfig") << "invalid filter scales" << std::endl;
0161 
0162   filterScales = std::make_unique<fftjet::EquidistantInLogSpace>(minScale, maxScale, nScales);
0163 
0164   percentileData.resize(nScales * nPercentiles);
0165 
0166   produces<reco::DiscretizedEnergyFlow>(outputLabel);
0167   produces<std::pair<double, double>>(outputLabel);
0168 }
0169 
0170 void FFTJetPileupProcessor::buildKernelConvolver(const edm::ParameterSet& ps) {
0171   // Get the eta and phi scales for the kernel(s)
0172   double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
0173   const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
0174   if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
0175     throw cms::Exception("FFTJetBadConfig") << "invalid kernel scales" << std::endl;
0176 
0177   if (convolverMinBin || convolverMaxBin)
0178     if (convolverMinBin >= convolverMaxBin || convolverMaxBin > energyFlow->nEta())
0179       throw cms::Exception("FFTJetBadConfig") << "invalid convolver bin range" << std::endl;
0180 
0181   // FFT assumes that the grid extent in eta is 2*Pi. Adjust the
0182   // kernel scale in eta to compensate.
0183   kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
0184 
0185   // Build the FFT engine
0186   engine = std::make_unique<MyFFTEngine>(energyFlow->nEta(), energyFlow->nPhi());
0187 
0188   // 2d kernel
0189   kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
0190       new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
0191 
0192   // 2d convolver
0193   convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real>>(new fftjet::FrequencyKernelConvolver<Real, Complex>(
0194       engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
0195 }
0196 
0197 FFTJetPileupProcessor::~FFTJetPileupProcessor() {}
0198 
0199 // ------------ method called to produce the data  ------------
0200 void FFTJetPileupProcessor::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0201   loadInputCollection(iEvent);
0202   discretizeEnergyFlow();
0203 
0204   // Determine the average Et density for this event.
0205   // Needs to be done here, before mixing in another grid.
0206   const fftjet::Grid2d<Real>& g(*energyFlow);
0207   const double densityBeforeMixing = g.sum() / pileupEtaPhiArea;
0208 
0209   // Mix an extra grid (if requested)
0210   double densityAfterMixing = -1.0;
0211   if (!externalGridFiles.empty()) {
0212     mixExtraGrid();
0213     densityAfterMixing = g.sum() / pileupEtaPhiArea;
0214   }
0215 
0216   // Various useful variables
0217   const unsigned nScales = filterScales->size();
0218   const double* scales = &(*filterScales)[0];
0219   Real* convData = const_cast<Real*>(convolvedFlow->data());
0220   Real* sortData = convData + convolverMinBin * convolvedFlow->nPhi();
0221   const unsigned dataLen = convolverMaxBin ? (convolverMaxBin - convolverMinBin) * convolvedFlow->nPhi()
0222                                            : convolvedFlow->nEta() * convolvedFlow->nPhi();
0223 
0224   // Load flattenning factors from DB (if requested)
0225   if (loadFlatteningFactorsFromDB)
0226     loadFlatteningFactors(iSetup);
0227 
0228   // Go over all scales and perform the convolutions
0229   convolver->setEventData(g.data(), g.nEta(), g.nPhi());
0230   for (unsigned iscale = 0; iscale < nScales; ++iscale) {
0231     // Perform the convolution
0232     convolver->convolveWithKernel(scales[iscale], convData, convolvedFlow->nEta(), convolvedFlow->nPhi());
0233 
0234     // Apply the flattening factors
0235     if (!etaFlatteningFactors.empty())
0236       convolvedFlow->scaleData(&etaFlatteningFactors[0], etaFlatteningFactors.size());
0237 
0238     // Sort the convolved data
0239     std::sort(sortData, sortData + dataLen);
0240 
0241     // Determine the percentile points
0242     for (unsigned iper = 0; iper < nPercentiles; ++iper) {
0243       // Map percentile 0 into point number 0,
0244       // 1 into point number dataLen-1
0245       const double q = (iper + 0.5) / nPercentiles;
0246       const double dindex = q * (dataLen - 1U);
0247       const unsigned ilow = static_cast<unsigned>(std::floor(dindex));
0248       const double percentile =
0249           fftjet::lin_interpolate_1d(ilow, ilow + 1U, sortData[ilow], sortData[ilow + 1U], dindex);
0250 
0251       // Store the calculated percentile
0252       percentileData[iscale * nPercentiles + iper] = percentile;
0253     }
0254   }
0255 
0256   // Convert percentile data into a more convenient storable object
0257   // and put it into the event record
0258   iEvent.put(std::make_unique<reco::DiscretizedEnergyFlow>(
0259                  &percentileData[0], "FFTJetPileupProcessor", -0.5, nScales - 0.5, 0.0, nScales, nPercentiles),
0260              outputLabel);
0261 
0262   iEvent.put(std::make_unique<std::pair<double, double>>(densityBeforeMixing, densityAfterMixing), outputLabel);
0263 }
0264 
0265 void FFTJetPileupProcessor::mixExtraGrid() {
0266   const unsigned nFiles = externalGridFiles.size();
0267   if (currentFileNum > nFiles) {
0268     // This is the first time this function is called
0269     currentFileNum = 0;
0270     gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
0271     if (!gridStream.is_open())
0272       throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0273                                                  " failed to open external grid file "
0274                                               << externalGridFiles[currentFileNum] << std::endl;
0275   }
0276 
0277   const fftjet::Grid2d<float>* g = nullptr;
0278   const unsigned maxFail = 100U;
0279   unsigned nEnergyRejected = 0;
0280 
0281   while (!g) {
0282     g = fftjet::Grid2d<float>::read(gridStream);
0283 
0284     // If we can't read the grid, we need to switch to another file
0285     for (unsigned ntries = 0; ntries < nFiles && g == nullptr; ++ntries) {
0286       gridStream.close();
0287       currentFileNum = (currentFileNum + 1U) % nFiles;
0288       gridStream.open(externalGridFiles[currentFileNum].c_str(), std::ios_base::in | std::ios_base::binary);
0289       if (!gridStream.is_open())
0290         throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0291                                                    " failed to open external grid file "
0292                                                 << externalGridFiles[currentFileNum] << std::endl;
0293       g = fftjet::Grid2d<float>::read(gridStream);
0294     }
0295 
0296     if (g)
0297       if (g->sum() > externalGridMaxEnergy) {
0298         delete g;
0299         g = nullptr;
0300         if (++nEnergyRejected >= maxFail)
0301           throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0302                                                      " too many grids in a row ("
0303                                                   << nEnergyRejected << ") failed the maximum energy cut" << std::endl;
0304       }
0305   }
0306 
0307   if (g) {
0308     add_Grid2d_data(energyFlow.get(), *g);
0309     delete g;
0310   } else {
0311     // Too bad, no useful file found
0312     throw cms::Exception("FFTJetBadConfig") << "ERROR in FFTJetPileupProcessor::mixExtraGrid():"
0313                                                " no valid grid records found"
0314                                             << std::endl;
0315   }
0316 }
0317 
0318 void FFTJetPileupProcessor::loadFlatteningFactors(const edm::EventSetup& iSetup) {
0319   edm::ESHandle<FFTJetLookupTableSequence> h;
0320   StaticFFTJetLookupTableSequenceLoader::instance().load(iSetup, flatteningTableRecord, h);
0321   std::shared_ptr<npstat::StorableMultivariateFunctor> f = (*h)[flatteningTableCategory][flatteningTableName];
0322 
0323   // Fill out the table of flattening factors as a function of eta
0324   const unsigned nEta = energyFlow->nEta();
0325   etaFlatteningFactors.clear();
0326   etaFlatteningFactors.reserve(nEta);
0327   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
0328     const double eta = energyFlow->etaBinCenter(ieta);
0329     etaFlatteningFactors.push_back((*f)(&eta, 1U));
0330   }
0331 }
0332 
0333 //define this as a plug-in
0334 DEFINE_FWK_MODULE(FFTJetPileupProcessor);