FFTJetProducers/plugins/FFTJetEFlowSmoother.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FFTJetProducers
0004 // Class:      FFTJetEFlowSmoother
0005 //
0006 /**\class FFTJetEFlowSmoother FFTJetEFlowSmoother.cc RecoJets/FFTJetProducers/plugins/FFTJetEFlowSmoother.cc
0007
0008  Description: Runs FFTJet filtering code for multiple scales and saves the results
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Thu Jun  2 18:49:49 CDT 2011
0016 //
0017 //
0018
0019 #include <cmath>
0020 #include <memory>
0021
0022 // FFTJet headers
0023 #include "fftjet/FrequencyKernelConvolver.hh"
0024 #include "fftjet/DiscreteGauss2d.hh"
0025 #include "fftjet/EquidistantSequence.hh"
0026 #include "fftjet/interpolate.hh"
0027 #include "fftjet/FrequencyKernelConvolver.hh"
0028 #include "fftjet/FrequencySequentialConvolver.hh"
0029 #include "fftjet/DiscreteGauss1d.hh"
0030 #include "fftjet/DiscreteGauss2d.hh"
0031
0032 // framework include files
0033 #include "FWCore/Framework/interface/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036
0037 #include "DataFormats/Common/interface/View.h"
0038 #include <TH3F.h>
0039
0040 // parameter parser header
0041 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0042
0043 // useful utilities collected in the second base
0044 #include "RecoJets/FFTJetProducers/interface/FFTJetInterface.h"
0045
0046 using namespace fftjetcms;
0047
0048 //
0049 // class declaration
0050 //
0051 class FFTJetEFlowSmoother : public FFTJetInterface {
0052 public:
0053   explicit FFTJetEFlowSmoother(const edm::ParameterSet&);
0054   FFTJetEFlowSmoother() = delete;
0055   FFTJetEFlowSmoother(const FFTJetEFlowSmoother&) = delete;
0056   FFTJetEFlowSmoother& operator=(const FFTJetEFlowSmoother&) = delete;
0057   ~FFTJetEFlowSmoother() override;
0058
0059 protected:
0060   // methods
0061   void produce(edm::Event&, const edm::EventSetup&) override;
0062
0063 private:
0064   void buildKernelConvolver(const edm::ParameterSet&);
0065
0066   // Storage for convolved energy flow
0067   std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> > convolvedFlow;
0068
0069   // Filtering scales
0070   std::unique_ptr<std::vector<double> > iniScales;
0071
0072   // The FFT engine(s)
0073   std::unique_ptr<MyFFTEngine> engine;
0074   std::unique_ptr<MyFFTEngine> anotherEngine;
0075
0076   // The pattern recognition kernel(s)
0077   std::unique_ptr<fftjet::AbsFrequencyKernel> kernel2d;
0078   std::unique_ptr<fftjet::AbsFrequencyKernel1d> etaKernel;
0079   std::unique_ptr<fftjet::AbsFrequencyKernel1d> phiKernel;
0080
0081   // The kernel convolver
0082   std::unique_ptr<fftjet::AbsConvolverBase<Real> > convolver;
0083
0084   // The scale power to use for scaling Et
0085   double scalePower;
0086
0087   // Overall factor to multiply Et with
0088   double etConversionFactor;
0089 };
0090
0091 //
0092 // constructors and destructor
0093 //
0094 FFTJetEFlowSmoother::FFTJetEFlowSmoother(const edm::ParameterSet& ps)
0095     : FFTJetInterface(ps),
0096       scalePower(ps.getParameter<double>("scalePower")),
0097       etConversionFactor(ps.getParameter<double>("etConversionFactor")) {
0098   // Build the discretization grid
0099   energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0100   checkConfig(energyFlow, "invalid discretization grid");
0101
0102   // Copy of the grid which will be used for convolutions
0103   convolvedFlow = std::make_unique<fftjet::Grid2d<fftjetcms::Real> >(*energyFlow);
0104
0105   // Build the initial set of pattern recognition scales
0106   iniScales = fftjet_ScaleSet_parser(ps.getParameter<edm::ParameterSet>("InitialScales"));
0107   checkConfig(iniScales, "invalid set of scales");
0108
0109   // Build the FFT engine(s), pattern recognition kernel(s),
0110   // and the kernel convolver
0111   buildKernelConvolver(ps);
0112
0113   // Build the set of pattern recognition scales
0114   produces<TH3F>(outputLabel);
0115 }
0116
0117 void FFTJetEFlowSmoother::buildKernelConvolver(const edm::ParameterSet& ps) {
0118   // Check the parameter named "etaDependentScaleFactors". If the vector
0119   // of scales is empty we will use 2d kernel, otherwise use 1d kernels
0120   const std::vector<double> etaDependentScaleFactors(ps.getParameter<std::vector<double> >("etaDependentScaleFactors"));
0121
0122   // Make sure that the number of scale factors provided is correct
0123   const bool use2dKernel = etaDependentScaleFactors.empty();
0124   if (!use2dKernel)
0125     if (etaDependentScaleFactors.size() != energyFlow->nEta())
0126       throw cms::Exception("FFTJetBadConfig") << "invalid number of eta-dependent scale factors" << std::endl;
0127
0128   // Get the eta and phi scales for the kernel(s)
0129   double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
0130   const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
0131   if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
0132     throw cms::Exception("FFTJetBadConfig") << "invalid kernel scale" << std::endl;
0133
0134   // FFT assumes that the grid extent in eta is 2*Pi. Adjust the
0135   // kernel scale in eta to compensate.
0136   kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
0137
0138   // Are we going to try to fix the efficiency near detector edges?
0139   const bool fixEfficiency = ps.getParameter<bool>("fixEfficiency");
0140
0141   // Minimum and maximum eta bin for the convolver
0142   unsigned convolverMinBin = 0, convolverMaxBin = 0;
0143   if (fixEfficiency || !use2dKernel) {
0144     convolverMinBin = ps.getParameter<unsigned>("convolverMinBin");
0145     convolverMaxBin = ps.getParameter<unsigned>("convolverMaxBin");
0146   }
0147
0148   if (use2dKernel) {
0149     // Build the FFT engine
0150     engine = std::make_unique<MyFFTEngine>(energyFlow->nEta(), energyFlow->nPhi());
0151
0152     // 2d kernel
0153     kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
0154         new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
0155
0156     // 2d convolver
0157     convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(new fftjet::FrequencyKernelConvolver<Real, Complex>(
0158         engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
0159   } else {
0160     // Need separate FFT engines for eta and phi
0161     engine = std::make_unique<MyFFTEngine>(1, energyFlow->nEta());
0162     anotherEngine = std::make_unique<MyFFTEngine>(1, energyFlow->nPhi());
0163
0164     // 1d kernels
0165     etaKernel =
0166         std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelEtaScale, energyFlow->nEta()));
0167
0168     phiKernel =
0169         std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelPhiScale, energyFlow->nPhi()));
0170
0171     // Sequential convolver
0172     convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(
0173         new fftjet::FrequencySequentialConvolver<Real, Complex>(engine.get(),
0174                                                                 anotherEngine.get(),
0175                                                                 etaKernel.get(),
0176                                                                 phiKernel.get(),
0177                                                                 etaDependentScaleFactors,
0178                                                                 convolverMinBin,
0179                                                                 convolverMaxBin,
0180                                                                 fixEfficiency));
0181   }
0182 }
0183
0184 FFTJetEFlowSmoother::~FFTJetEFlowSmoother() {}
0185
0186 // ------------ method called to produce the data  ------------
0187 void FFTJetEFlowSmoother::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0188   loadInputCollection(iEvent);
0189   discretizeEnergyFlow();
0190
0191   // Various useful variables
0192   const fftjet::Grid2d<Real>& g(*energyFlow);
0193   const unsigned nScales = iniScales->size();
0194   const double* scales = &(*iniScales)[0];
0195   Real* convData = const_cast<Real*>(convolvedFlow->data());
0196   const unsigned nEta = g.nEta();
0197   const unsigned nPhi = g.nPhi();
0198   const double bin0edge = g.phiBin0Edge();
0199
0200   // We will fill the following histo
0201   auto pTable = std::make_unique<TH3F>("FFTJetEFlowSmoother",
0202                                        "FFTJetEFlowSmoother",
0203                                        nScales + 1U,
0204                                        -1.5,
0205                                        nScales - 0.5,
0206                                        nEta,
0207                                        g.etaMin(),
0208                                        g.etaMax(),
0209                                        nPhi,
0210                                        bin0edge,
0211                                        bin0edge + 2.0 * M_PI);
0212   TH3F* h = pTable.get();
0213   h->SetDirectory(nullptr);
0214   h->GetXaxis()->SetTitle("Scale");
0215   h->GetYaxis()->SetTitle("Eta");
0216   h->GetZaxis()->SetTitle("Phi");
0217
0218   // Fill the original thing
0219   double factor = etConversionFactor * pow(getEventScale(), scalePower);
0220   for (unsigned ieta = 0; ieta < nEta; ++ieta)
0221     for (unsigned iphi = 0; iphi < nPhi; ++iphi)
0222       h->SetBinContent(1U, ieta + 1U, iphi + 1U, factor * g.binValue(ieta, iphi));
0223
0224   // Go over all scales and perform the convolutions
0225   convolver->setEventData(g.data(), nEta, nPhi);
0226   for (unsigned iscale = 0; iscale < nScales; ++iscale) {
0227     factor = etConversionFactor * pow(scales[iscale], scalePower);
0228
0229     // Perform the convolution
0230     convolver->convolveWithKernel(scales[iscale], convData, convolvedFlow->nEta(), convolvedFlow->nPhi());
0231
0232     // Fill the output histo
0233     for (unsigned ieta = 0; ieta < nEta; ++ieta) {
0234       const Real* etaData = convData + ieta * nPhi;
0235       for (unsigned iphi = 0; iphi < nPhi; ++iphi)
0236         h->SetBinContent(iscale + 2U, ieta + 1U, iphi + 1U, factor * etaData[iphi]);
0237     }
0238   }
0239
0240   iEvent.put(std::move(pTable), outputLabel);
0241 }
0242
0243 //define this as a plug-in
0244 DEFINE_FWK_MODULE(FFTJetEFlowSmoother);