FFTJetProducers/plugins/FFTJetPatRecoProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FFTJetProducers
0004 // Class:      FFTJetPatRecoProducer
0005 //
0006 /**\class FFTJetPatRecoProducer FFTJetPatRecoProducer.cc RecoJets/FFTJetProducer/plugins/FFTJetPatRecoProducer.cc
0007
0008  Description: Runs FFTJet pattern recognition stage and saves the results
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Tue Jun 15 12:45:45 CDT 2010
0016 //
0017 //
0018
0019 #include <fstream>
0020 #include <memory>
0021
0022 // FFTJet headers
0023 #include "fftjet/ProximityClusteringTree.hh"
0024 #include "fftjet/ClusteringSequencer.hh"
0025 #include "fftjet/ClusteringTreeSparsifier.hh"
0026 #include "fftjet/FrequencyKernelConvolver.hh"
0027 #include "fftjet/FrequencySequentialConvolver.hh"
0028 #include "fftjet/DiscreteGauss1d.hh"
0029 #include "fftjet/DiscreteGauss2d.hh"
0030
0031 // framework include files
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include "DataFormats/Common/interface/View.h"
0037 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0038 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0039 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0040 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
0041 #include "DataFormats/JetReco/interface/TrackJetCollection.h"
0042 #include "DataFormats/Candidate/interface/Candidate.h"
0043
0044 // Energy flow object
0045 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0046
0047 // parameter parser header
0048 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0049
0050 // functions which manipulate storable trees
0051 #include "RecoJets/FFTJetAlgorithms/interface/clusteringTreeConverters.h"
0052
0053 // functions which manipulate energy discretization grids
0054 #include "RecoJets/FFTJetAlgorithms/interface/gridConverters.h"
0055
0056 // useful utilities collected in the second base
0057 #include "RecoJets/FFTJetProducers/interface/FFTJetInterface.h"
0058
0059 using namespace fftjetcms;
0060
0061 //
0062 // class declaration
0063 //
0064 class FFTJetPatRecoProducer : public FFTJetInterface {
0065 public:
0066   explicit FFTJetPatRecoProducer(const edm::ParameterSet&);
0067   FFTJetPatRecoProducer() = delete;
0068   FFTJetPatRecoProducer(const FFTJetPatRecoProducer&) = delete;
0069   FFTJetPatRecoProducer& operator=(const FFTJetPatRecoProducer&) = delete;
0070   ~FFTJetPatRecoProducer() override;
0071
0072 protected:
0073   // Useful local typedefs
0074   typedef fftjet::ProximityClusteringTree<fftjet::Peak, long> ClusteringTree;
0075   typedef fftjet::SparseClusteringTree<fftjet::Peak, long> SparseTree;
0076   typedef fftjet::ClusteringSequencer<Real> Sequencer;
0077   typedef fftjet::ClusteringTreeSparsifier<fftjet::Peak, long> Sparsifier;
0078
0079   // methods
0080   void produce(edm::Event&, const edm::EventSetup&) override;
0081
0082   void buildKernelConvolver(const edm::ParameterSet&);
0083   fftjet::PeakFinder buildPeakFinder(const edm::ParameterSet&);
0084
0085   template <class Real>
0086   void buildSparseProduct(edm::Event&) const;
0087
0088   template <class Real>
0089   void buildDenseProduct(edm::Event&) const;
0090
0091   // The complete clustering tree
0092   ClusteringTree* clusteringTree;
0093
0094   // Basically, we need to create FFTJet objects
0095   // ClusteringSequencer and ClusteringTreeSparsifier
0096   // which will subsequently perform most of the work
0097   std::unique_ptr<Sequencer> sequencer;
0098   std::unique_ptr<Sparsifier> sparsifier;
0099
0100   // The FFT engine(s)
0101   std::unique_ptr<MyFFTEngine> engine;
0102   std::unique_ptr<MyFFTEngine> anotherEngine;
0103
0104   // The pattern recognition kernel(s)
0105   std::unique_ptr<fftjet::AbsFrequencyKernel> kernel2d;
0106   std::unique_ptr<fftjet::AbsFrequencyKernel1d> etaKernel;
0107   std::unique_ptr<fftjet::AbsFrequencyKernel1d> phiKernel;
0108
0109   // The kernel convolver
0110   std::unique_ptr<fftjet::AbsConvolverBase<Real> > convolver;
0111
0112   // The peak selector for the clustering tree
0113   std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak> > peakSelector;
0114
0115   // Distance calculator for the clustering tree
0116   std::unique_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak> > distanceCalc;
0117
0118   // The sparse clustering tree
0119   SparseTree sparseTree;
0120
0121   // The following parameters will define the behavior
0122   // of the algorithm wrt insertion of the complete event
0123   // into the clustering tree
0124   const double completeEventDataCutoff;
0125
0126   // Are we going to make clustering trees?
0127   const bool makeClusteringTree;
0128
0129   // Are we going to verify the data conversion for double precision
0130   // storage?
0131   const bool verifyDataConversion;
0132
0133   // Are we going to store the discretization grid?
0134   const bool storeDiscretizationGrid;
0135
0136   // Sparsify the clustering tree?
0137   const bool sparsify;
0138
0139 private:
0140   // Members needed for storing grids externally
0141   std::ofstream externalGridStream;
0142   bool storeGridsExternally;
0143   fftjet::Grid2d<float>* extGrid;
0144 };
0145
0146 //
0147 // constructors and destructor
0148 //
0149 FFTJetPatRecoProducer::FFTJetPatRecoProducer(const edm::ParameterSet& ps)
0150     : FFTJetInterface(ps),
0151       clusteringTree(nullptr),
0152       completeEventDataCutoff(ps.getParameter<double>("completeEventDataCutoff")),
0153       makeClusteringTree(ps.getParameter<bool>("makeClusteringTree")),
0154       verifyDataConversion(ps.getUntrackedParameter<bool>("verifyDataConversion", false)),
0155       storeDiscretizationGrid(ps.getParameter<bool>("storeDiscretizationGrid")),
0156       sparsify(ps.getParameter<bool>("sparsify")),
0157       extGrid(nullptr) {
0158   // register your products
0159   if (makeClusteringTree) {
0160     if (storeInSinglePrecision())
0161       produces<reco::PattRecoTree<float, reco::PattRecoPeak<float> > >(outputLabel);
0162     else
0163       produces<reco::PattRecoTree<double, reco::PattRecoPeak<double> > >(outputLabel);
0164   }
0165   if (storeDiscretizationGrid)
0166     produces<reco::DiscretizedEnergyFlow>(outputLabel);
0167
0168   // Check if we want to write the grids into an external file
0169   const std::string externalGridFile(ps.getParameter<std::string>("externalGridFile"));
0170   storeGridsExternally = !externalGridFile.empty();
0171   if (storeGridsExternally) {
0172     externalGridStream.open(externalGridFile.c_str(), std::ios_base::out | std::ios_base::binary);
0173     if (!externalGridStream.is_open())
0174       throw cms::Exception("FFTJetBadConfig")
0175           << "FFTJetPatRecoProducer failed to open file " << externalGridFile << std::endl;
0176   }
0177
0178   if (!makeClusteringTree && !storeDiscretizationGrid && !storeGridsExternally) {
0179     throw cms::Exception("FFTJetBadConfig")
0180         << "FFTJetPatRecoProducer is not configured to produce anything" << std::endl;
0181   }
0182
0183   // now do whatever other initialization is needed
0184
0185   // Build the discretization grid
0186   energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0187   checkConfig(energyFlow, "invalid discretization grid");
0188
0189   // Build the FFT engine(s), pattern recognition kernel(s),
0190   // and the kernel convolver
0191   buildKernelConvolver(ps);
0192
0193   // Build the peak selector
0194   peakSelector = fftjet_PeakSelector_parser(ps.getParameter<edm::ParameterSet>("PeakSelectorConfiguration"));
0195   checkConfig(peakSelector, "invalid peak selector");
0196
0197   // Build the initial set of pattern recognition scales
0198   std::unique_ptr<std::vector<double> > iniScales =
0199       fftjet_ScaleSet_parser(ps.getParameter<edm::ParameterSet>("InitialScales"));
0200   checkConfig(iniScales, "invalid set of scales");
0201
0202   // Do we want to use the adaptive clustering tree algorithm?
0203   const unsigned maxAdaptiveScales = ps.getParameter<unsigned>("maxAdaptiveScales");
0204   const double minAdaptiveRatioLog = ps.getParameter<double>("minAdaptiveRatioLog");
0205   if (minAdaptiveRatioLog <= 0.0)
0206     throw cms::Exception("FFTJetBadConfig") << "bad adaptive ratio logarithm limit" << std::endl;
0207
0208   // Make sure that all standard scales are larger than the
0209   // complete event scale
0210   if (getEventScale() > 0.0) {
0211     const double cs = getEventScale();
0212     const unsigned nscales = iniScales->size();
0213     for (unsigned i = 0; i < nscales; ++i)
0214       if (cs >= (*iniScales)[i])
0215         throw cms::Exception("FFTJetBadConfig") << "incompatible scale for complete event" << std::endl;
0216   }
0217
0218   // At this point we are ready to build the clustering sequencer
0219   sequencer = std::make_unique<Sequencer>(
0220       convolver.get(), peakSelector.get(), buildPeakFinder(ps), *iniScales, maxAdaptiveScales, minAdaptiveRatioLog);
0221
0222   // Build the clustering tree sparsifier
0223   const edm::ParameterSet& SparsifierConfiguration(ps.getParameter<edm::ParameterSet>("SparsifierConfiguration"));
0224   sparsifier = fftjet_ClusteringTreeSparsifier_parser(SparsifierConfiguration);
0225   checkConfig(sparsifier, "invalid sparsifier parameters");
0226
0227   // Build the distance calculator for the clustering tree
0228   const edm::ParameterSet& TreeDistanceCalculator(ps.getParameter<edm::ParameterSet>("TreeDistanceCalculator"));
0229   distanceCalc = fftjet_DistanceCalculator_parser(TreeDistanceCalculator);
0230   checkConfig(distanceCalc, "invalid tree distance calculator");
0231
0232   // Build the clustering tree itself
0233   clusteringTree = new ClusteringTree(distanceCalc.get());
0234 }
0235
0236 void FFTJetPatRecoProducer::buildKernelConvolver(const edm::ParameterSet& ps) {
0237   // Check the parameter named "etaDependentScaleFactors". If the vector
0238   // of scales is empty we will use 2d kernel, otherwise use 1d kernels
0239   const std::vector<double> etaDependentScaleFactors(ps.getParameter<std::vector<double> >("etaDependentScaleFactors"));
0240
0241   // Make sure that the number of scale factors provided is correct
0242   const bool use2dKernel = etaDependentScaleFactors.empty();
0243   if (!use2dKernel)
0244     if (etaDependentScaleFactors.size() != energyFlow->nEta())
0245       throw cms::Exception("FFTJetBadConfig") << "invalid number of eta-dependent scale factors" << std::endl;
0246
0247   // Get the eta and phi scales for the kernel(s)
0248   double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
0249   const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
0250   if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
0251     throw cms::Exception("FFTJetBadConfig") << "invalid kernel scale" << std::endl;
0252
0253   // FFT assumes that the grid extent in eta is 2*Pi. Adjust the
0254   // kernel scale in eta to compensate.
0255   kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
0256
0257   // Are we going to try to fix the efficiency near detector edges?
0258   const bool fixEfficiency = ps.getParameter<bool>("fixEfficiency");
0259
0260   // Minimum and maximum eta bin for the convolver
0261   unsigned convolverMinBin = 0, convolverMaxBin = 0;
0262   if (fixEfficiency || !use2dKernel) {
0263     convolverMinBin = ps.getParameter<unsigned>("convolverMinBin");
0264     convolverMaxBin = ps.getParameter<unsigned>("convolverMaxBin");
0265   }
0266
0267   if (use2dKernel) {
0268     // Build the FFT engine
0269     engine = std::make_unique<MyFFTEngine>(energyFlow->nEta(), energyFlow->nPhi());
0270
0271     // 2d kernel
0272     kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
0273         new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
0274
0275     // 2d convolver
0276     convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(new fftjet::FrequencyKernelConvolver<Real, Complex>(
0277         engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
0278   } else {
0279     // Need separate FFT engines for eta and phi
0280     engine = std::make_unique<MyFFTEngine>(1, energyFlow->nEta());
0281     anotherEngine = std::make_unique<MyFFTEngine>(1, energyFlow->nPhi());
0282
0283     // 1d kernels
0284     etaKernel =
0285         std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelEtaScale, energyFlow->nEta()));
0286
0287     phiKernel =
0288         std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelPhiScale, energyFlow->nPhi()));
0289
0290     // Sequential convolver
0291     convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(
0292         new fftjet::FrequencySequentialConvolver<Real, Complex>(engine.get(),
0293                                                                 anotherEngine.get(),
0294                                                                 etaKernel.get(),
0295                                                                 phiKernel.get(),
0296                                                                 etaDependentScaleFactors,
0297                                                                 convolverMinBin,
0298                                                                 convolverMaxBin,
0299                                                                 fixEfficiency));
0300   }
0301 }
0302
0303 fftjet::PeakFinder FFTJetPatRecoProducer::buildPeakFinder(const edm::ParameterSet& ps) {
0304   const double peakFinderMaxEta = ps.getParameter<double>("peakFinderMaxEta");
0305   if (peakFinderMaxEta <= 0.0)
0306     throw cms::Exception("FFTJetBadConfig") << "invalid peak finder eta cut" << std::endl;
0307   const double maxMagnitude = ps.getParameter<double>("peakFinderMaxMagnitude");
0308   int minBin = energyFlow->getEtaBin(-peakFinderMaxEta);
0309   if (minBin < 0)
0310     minBin = 0;
0311   int maxBin = energyFlow->getEtaBin(peakFinderMaxEta) + 1;
0312   if (maxBin < 0)
0313     maxBin = 0;
0314   return fftjet::PeakFinder(maxMagnitude, true, minBin, maxBin);
0315 }
0316
0317 FFTJetPatRecoProducer::~FFTJetPatRecoProducer() {
0318   // do anything here that needs to be done at desctruction time
0319   // (e.g. close files, deallocate resources etc.)
0320   if (storeGridsExternally)
0321     externalGridStream.close();
0322   delete clusteringTree;
0323   delete extGrid;
0324 }
0325
0326 //
0327 // member functions
0328 //
0329 template <class Real>
0330 void FFTJetPatRecoProducer::buildSparseProduct(edm::Event& ev) const {
0331   typedef reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > StoredTree;
0332
0333   auto tree = std::make_unique<StoredTree>();
0334
0335   sparsePeakTreeToStorable(sparseTree, sequencer->maxAdaptiveScales(), tree.get());
0336
0337   // Check that we can restore the tree
0338   if (verifyDataConversion && !storeInSinglePrecision()) {
0339     SparseTree check;
0340     const std::vector<double>& scalesUsed(sequencer->getInitialScales());
0341     sparsePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
0342     if (sparseTree != check)
0343       throw cms::Exception("FFTJetInterface") << "Data conversion failed for sparse clustering tree" << std::endl;
0344   }
0345
0346   ev.put(std::move(tree), outputLabel);
0347 }
0348
0349 template <class Real>
0350 void FFTJetPatRecoProducer::buildDenseProduct(edm::Event& ev) const {
0351   typedef reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > StoredTree;
0352
0353   auto tree = std::make_unique<StoredTree>();
0354
0355   densePeakTreeToStorable(*clusteringTree, sequencer->maxAdaptiveScales(), tree.get());
0356
0357   // Check that we can restore the tree
0358   if (verifyDataConversion && !storeInSinglePrecision()) {
0359     ClusteringTree check(distanceCalc.get());
0360     const std::vector<double>& scalesUsed(sequencer->getInitialScales());
0361     densePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
0362     if (*clusteringTree != check)
0363       throw cms::Exception("FFTJetInterface") << "Data conversion failed for dense clustering tree" << std::endl;
0364   }
0365
0366   ev.put(std::move(tree), outputLabel);
0367 }
0368
0369 // ------------ method called to produce the data  ------------
0370 void FFTJetPatRecoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0371   loadInputCollection(iEvent);
0372   discretizeEnergyFlow();
0373
0374   if (makeClusteringTree) {
0375     sequencer->run(*energyFlow, clusteringTree);
0376     if (getEventScale() > 0.0)
0377       sequencer->insertCompleteEvent(getEventScale(), *energyFlow, clusteringTree, completeEventDataCutoff);
0378
0379     if (sparsify) {
0380       sparsifier->sparsify(*clusteringTree, &sparseTree);
0381
0382       // Do not call the "sortNodes" method of the sparse tree here.
0383       // Currently, the nodes are sorted by daughter number.
0384       // This is the way we want it in storage because the stored
0385       // tree does not include daughter ordering info explicitly.
0386
0387       if (storeInSinglePrecision())
0388         buildSparseProduct<float>(iEvent);
0389       else
0390         buildSparseProduct<double>(iEvent);
0391     } else {
0392       if (storeInSinglePrecision())
0393         buildDenseProduct<float>(iEvent);
0394       else
0395         buildDenseProduct<double>(iEvent);
0396     }
0397   }
0398
0399   if (storeDiscretizationGrid) {
0400     const fftjet::Grid2d<Real>& g(*energyFlow);
0401
0402     auto flow = std::make_unique<reco::DiscretizedEnergyFlow>(
0403         g.data(), g.title(), g.etaMin(), g.etaMax(), g.phiBin0Edge(), g.nEta(), g.nPhi());
0404
0405     if (verifyDataConversion) {
0406       fftjet::Grid2d<Real> check(
0407           flow->nEtaBins(), flow->etaMin(), flow->etaMax(), flow->nPhiBins(), flow->phiBin0Edge(), flow->title());
0408       check.blockSet(flow->data(), flow->nEtaBins(), flow->nPhiBins());
0409       assert(g == check);
0410     }
0411
0412     iEvent.put(std::move(flow), outputLabel);
0413   }
0414
0415   if (storeGridsExternally) {
0416     if (extGrid)
0417       copy_Grid2d_data(extGrid, *energyFlow);
0418     else
0419       extGrid = convert_Grid2d_to_float(*energyFlow);
0420     if (!extGrid->write(externalGridStream)) {
0421       throw cms::Exception("FFTJetPatRecoProducer::produce")
0422           << "Failed to write grid data into an external file" << std::endl;
0423     }
0424   }
0425 }
0426
0427 //define this as a plug-in
0428 DEFINE_FWK_MODULE(FFTJetPatRecoProducer);