FFTJetProducers/plugins/FFTJetProducer.h

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoJets/FFTJetProducers
0004 // Class:      FFTJetProducer
0005 //
0006 /**\class FFTJetProducer FFTJetProducer.h RecoJets/FFTJetProducers/plugins/FFTJetProducer.h
0007
0008  Description: makes jets using FFTJet clustering tree
0009
0010  Implementation:
0011      If you want to change the jet algorithm functionality (for example,
0012      by providing your own jet membership function), derive from this
0013      class and override the appropriate parser method (for example,
0014      parse_jetMembershipFunction). At the end of your own parser, don't
0015      forget to call the parser of the base class in order to get the default
0016      behavior when your special configuration is not provided (this is known
0017      as the "chain-of-responsibility" design pattern). If you also need to
0018      override "beginJob" and/or "produce" methods, the first thing to do in
0019      your method is to call the corresponding method of this base.
0020 */
0021 //
0022 // Original Author:  Igor Volobouev
0023 //         Created:  Sun Jun 20 14:32:36 CDT 2010
0024 //
0025 //
0026
0027 #ifndef RecoJets_FFTJetProducers_FFTJetProducer_h
0028 #define RecoJets_FFTJetProducers_FFTJetProducer_h
0029
0030 #include <memory>
0031
0032 // FFTJet headers
0033 #include "fftjet/AbsRecombinationAlg.hh"
0034 #include "fftjet/AbsVectorRecombinationAlg.hh"
0035 #include "fftjet/SparseClusteringTree.hh"
0036
0037 // Additional FFTJet headers
0038 #include "fftjet/VectorRecombinationAlgFactory.hh"
0039 #include "fftjet/RecombinationAlgFactory.hh"
0040
0041 #include "DataFormats/JetReco/interface/FFTCaloJetCollection.h"
0042 #include "DataFormats/JetReco/interface/FFTGenJetCollection.h"
0043 #include "DataFormats/JetReco/interface/FFTPFJetCollection.h"
0044 #include "DataFormats/JetReco/interface/FFTJPTJetCollection.h"
0045 #include "DataFormats/JetReco/interface/FFTBasicJetCollection.h"
0046 #include "DataFormats/JetReco/interface/FFTTrackJetCollection.h"
0047 #include "DataFormats/JetReco/interface/FFTJetProducerSummary.h"
0048 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0049
0050 #include "RecoJets/FFTJetAlgorithms/interface/clusteringTreeConverters.h"
0051 #include "RecoJets/FFTJetAlgorithms/interface/jetConverters.h"
0052 #include "RecoJets/FFTJetAlgorithms/interface/matchOneToOne.h"
0053 #include "RecoJets/FFTJetAlgorithms/interface/JetToPeakDistance.h"
0054 #include "RecoJets/FFTJetAlgorithms/interface/adjustForPileup.h"
0055
0056 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0057
0058 // framework include files
0059 #include "FWCore/Framework/interface/Frameworkfwd.h"
0060 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0061 #include "FWCore/Framework/interface/Event.h"
0062
0063 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0064
0065 // local FFTJet-related definitions
0066 #include "RecoJets/FFTJetAlgorithms/interface/fftjetTypedefs.h"
0067 #include "RecoJets/FFTJetAlgorithms/interface/AbsPileupCalculator.h"
0068 #include "RecoJets/FFTJetProducers/interface/FFTJetInterface.h"
0069
0070 namespace fftjetcms {
0071   class DiscretizedEnergyFlow;
0072 }
0073
0074 class CaloGeometry;
0075 class CaloGeometryRecord;
0076 class HcalTopology;
0077 class HcalRecNumberingRecord;
0078
0079 //
0080 // class declaration
0081 //
0082 class FFTJetProducer : public fftjetcms::FFTJetInterface {
0083 public:
0084   typedef fftjet::RecombinedJet<fftjetcms::VectorLike> RecoFFTJet;
0085   typedef fftjet::SparseClusteringTree<fftjet::Peak, long> SparseTree;
0086
0087   // Masks for the status bits. Do not add anything
0088   // here -- higher bits (starting with 0x1000) will be
0089   // used to indicate jet correction levels applied.
0090   enum StatusBits {
0091     RESOLUTION = 0xff,
0092     CONSTITUENTS_RESUMMED = 0x100,
0093     PILEUP_CALCULATED = 0x200,
0094     PILEUP_SUBTRACTED_4VEC = 0x400,
0095     PILEUP_SUBTRACTED_PT = 0x800
0096   };
0097
0098   enum Resolution { FIXED = 0, MAXIMALLY_STABLE, GLOBALLY_ADAPTIVE, LOCALLY_ADAPTIVE, FROM_GENJETS };
0099
0100   explicit FFTJetProducer(const edm::ParameterSet&);
0101   // Explicitly disable other ways to construct this object
0102   FFTJetProducer() = delete;
0103   FFTJetProducer(const FFTJetProducer&) = delete;
0104   FFTJetProducer& operator=(const FFTJetProducer&) = delete;
0105
0106   ~FFTJetProducer() override;
0107
0108   // Parser for the resolution enum
0109   static Resolution parse_resolution(const std::string& name);
0110
0111 protected:
0112   // Functions which should be overriden from the base
0113   void beginStream(edm::StreamID) override;
0114   void produce(edm::Event&, const edm::EventSetup&) override;
0115
0116   // The following functions can be overriden by derived classes
0117   // in order to adjust jet reconstruction algorithm behavior.
0118
0119   // Override the following method in order to implement
0120   // your own precluster selection strategy
0121   virtual void selectPreclusters(const SparseTree& tree,
0122                                  const fftjet::Functor1<bool, fftjet::Peak>& peakSelector,
0123                                  std::vector<fftjet::Peak>* preclusters);
0124
0125   // Precluster maker from GenJets (useful in calibration)
0126   virtual void genJetPreclusters(const SparseTree& tree,
0127                                  edm::Event&,
0128                                  const edm::EventSetup&,
0129                                  const fftjet::Functor1<bool, fftjet::Peak>& peakSelector,
0130                                  std::vector<fftjet::Peak>* preclusters);
0131
0132   // Override the following method (which by default does not do
0133   // anything) in order to implement your own process-dependent
0134   // assignment of membership functions to preclusters. This method
0135   // will be called once per event, just before the main algorithm.
0136   virtual void assignMembershipFunctions(std::vector<fftjet::Peak>* preclusters);
0137
0138   // Parser for the peak selector
0139   virtual std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak> > parse_peakSelector(const edm::ParameterSet&);
0140
0141   // Parser for the default jet membership function
0142   virtual std::unique_ptr<fftjet::ScaleSpaceKernel> parse_jetMembershipFunction(const edm::ParameterSet&);
0143
0144   // Parser for the background membership function
0145   virtual std::unique_ptr<fftjetcms::AbsBgFunctor> parse_bgMembershipFunction(const edm::ParameterSet&);
0146
0147   // Calculator for the recombination scale
0148   virtual std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > parse_recoScaleCalcPeak(const edm::ParameterSet&);
0149
0150   // Calculator for the recombination scale ratio
0151   virtual std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > parse_recoScaleRatioCalcPeak(
0152       const edm::ParameterSet&);
0153
0154   // Calculator for the membership function factor
0155   virtual std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > parse_memberFactorCalcPeak(const edm::ParameterSet&);
0156
0157   // Similar calculators for the iterative algorithm
0158   virtual std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > parse_recoScaleCalcJet(const edm::ParameterSet&);
0159
0160   virtual std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > parse_recoScaleRatioCalcJet(const edm::ParameterSet&);
0161
0162   virtual std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > parse_memberFactorCalcJet(const edm::ParameterSet&);
0163
0164   // Calculator of the distance between jets which is used to make
0165   // the decision about convergence of the iterative algorithm
0166   virtual std::unique_ptr<fftjet::Functor2<double, RecoFFTJet, RecoFFTJet> > parse_jetDistanceCalc(
0167       const edm::ParameterSet&);
0168
0169   // Pile-up density calculator
0170   virtual std::unique_ptr<fftjetcms::AbsPileupCalculator> parse_pileupDensityCalc(const edm::ParameterSet& ps);
0171
0172   // The following function performs most of the precluster selection
0173   // work in this module. You might want to reuse it if only a slight
0174   // modification of the "selectPreclusters" method is desired.
0175   void selectTreeNodes(const SparseTree& tree,
0176                        const fftjet::Functor1<bool, fftjet::Peak>& peakSelect,
0177                        std::vector<SparseTree::NodeId>* nodes);
0178
0179 private:
0180   typedef fftjet::AbsVectorRecombinationAlg<fftjetcms::VectorLike, fftjetcms::BgData> RecoAlg;
0181   typedef fftjet::AbsRecombinationAlg<fftjetcms::Real, fftjetcms::VectorLike, fftjetcms::BgData> GridAlg;
0182
0183   // Useful local utilities
0184   template <class Real>
0185   void loadSparseTreeData(const edm::Event&,
0186                           const edm::EDGetTokenT<reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > >& tok);
0187
0188   void removeFakePreclusters();
0189
0190   // The following methods do most of the work.
0191   // The following function tells us if the grid was rebuilt.
0192   bool loadEnergyFlow(const edm::Event& iEvent, std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& flow);
0193   void buildGridAlg();
0194   void prepareRecombinationScales();
0195   bool checkConvergence(const std::vector<RecoFFTJet>& previousIterResult, std::vector<RecoFFTJet>& thisIterResult);
0196   void determineGriddedConstituents();
0197   void determineVectorConstituents();
0198   void saveResults(edm::Event& iEvent, const edm::EventSetup&, unsigned nPreclustersFound);
0199
0200   template <typename Jet>
0201   void writeJets(edm::Event& iEvent, const edm::EventSetup&);
0202
0203   template <typename Jet>
0204   void makeProduces(const std::string& alias, const std::string& tag);
0205
0206   // The following function scans the pile-up density
0207   // and fills the pile-up grid. Can be overriden if
0208   // necessary.
0209   virtual void determinePileupDensityFromConfig(const edm::Event& iEvent,
0210                                                 std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& density);
0211
0212   // Similar function for getting pile-up shape from the database
0213   virtual void determinePileupDensityFromDB(const edm::Event& iEvent,
0214                                             const edm::EventSetup& iSetup,
0215                                             std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& density);
0216
0217   // The following function builds the pile-up estimate
0218   // for each jet
0219   void determinePileup();
0220
0221   // The following function returns the number of iterations
0222   // performed. If this number equals to or less than "maxIterations"
0223   // then the iterations have converged. If the number larger than
0224   // "maxIterations" then the iterations failed to converge (note,
0225   // however, that only "maxIterations" iterations would still be
0226   // performed).
0227   unsigned iterateJetReconstruction();
0228
0229   // A function to set jet status bits
0230   static void setJetStatusBit(RecoFFTJet* jet, int mask, bool value);
0231
0232   //
0233   // ----------member data ---------------------------
0234   //
0235
0236   // Local copy of the module configuration
0237   const edm::ParameterSet myConfiguration;
0238
0239   // Label for the tree produced by FFTJetPatRecoProducer
0240   const edm::InputTag treeLabel;
0241
0242   // Are we going to use energy flow discretization grid as input
0243   // to jet reconstruction?
0244   const bool useGriddedAlgorithm;
0245
0246   // Are we going to rebuild the energy flow discretization grid
0247   // or to reuse the grid made by FFTJetPatRecoProducer?
0248   const bool reuseExistingGrid;
0249
0250   // Are we iterating?
0251   const unsigned maxIterations;
0252
0253   // Parameters which affect iteration convergence
0254   const unsigned nJetsRequiredToConverge;
0255   const double convergenceDistance;
0256
0257   // Are we building assignments of collection members to jets?
0258   const bool assignConstituents;
0259
0260   // Are we resumming the constituents to determine jet 4-vectors?
0261   // This might make sense if FFTJet is used in the crisp, gridded
0262   // mode to determine jet areas, and vector recombination is desired.
0263   const bool resumConstituents;
0264
0265   // Are we going to subtract the pile-up? Note that
0266   // pile-up subtraction does not modify eta and phi moments.
0267   const bool calculatePileup;
0268   const bool subtractPileup;
0269   const bool subtractPileupAs4Vec;
0270
0271   // Label for the pile-up energy flow. Must be specified
0272   // if the pile-up is subtracted.
0273   const edm::InputTag pileupLabel;
0274
0275   // Scale for the peak selection (if the scale is fixed)
0276   const double fixedScale;
0277
0278   // Minimum and maximum scale for searching stable configurations
0279   const double minStableScale;
0280   const double maxStableScale;
0281
0282   // Stability "alpha"
0283   const double stabilityAlpha;
0284
0285   // Not sure at this point how to treat noise... For now, this is
0286   // just a single configurable number...
0287   const double noiseLevel;
0288
0289   // Number of clusters requested (if the scale is adaptive)
0290   const unsigned nClustersRequested;
0291
0292   // Maximum eta for the grid-based algorithm
0293   const double gridScanMaxEta;
0294
0295   // Parameters related to the recombination algorithm
0296   const std::string recombinationAlgorithm;
0297   const bool isCrisp;
0298   const double unlikelyBgWeight;
0299   const double recombinationDataCutoff;
0300
0301   // Label for the genJets used as seeds for jets
0302   const edm::InputTag genJetsLabel;
0303
0304   // Maximum number of preclusters to use as jet seeds.
0305   // This does not take into account the preclusters
0306   // for which the value of the membership factor is 0.
0307   const unsigned maxInitialPreclusters;
0308
0309   // Resolution. The corresponding parameter value
0310   // should be one of "fixed", "maximallyStable",
0311   // "globallyAdaptive", "locallyAdaptive", or "fromGenJets".
0312   Resolution resolution;
0313
0314   // Parameters related to the pileup shape stored
0315   // in the database
0316   std::string pileupTableRecord;
0317   std::string pileupTableName;
0318   std::string pileupTableCategory;
0319   bool loadPileupFromDB;
0320
0321   // Scales used
0322   std::unique_ptr<std::vector<double> > iniScales;
0323
0324   // The sparse clustering tree
0325   SparseTree sparseTree;
0326
0327   // Peak selector for the peaks already in the tree
0328   std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak> > peakSelector;
0329
0330   // Recombination algorithms and related quantities
0331   std::unique_ptr<RecoAlg> recoAlg;
0332   std::unique_ptr<GridAlg> gridAlg;
0333   std::unique_ptr<fftjet::ScaleSpaceKernel> jetMembershipFunction;
0334   std::unique_ptr<fftjetcms::AbsBgFunctor> bgMembershipFunction;
0335
0336   // Calculator for the recombination scale
0337   std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > recoScaleCalcPeak;
0338
0339   // Calculator for the recombination scale ratio
0340   std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > recoScaleRatioCalcPeak;
0341
0342   // Calculator for the membership function factor
0343   std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > memberFactorCalcPeak;
0344
0345   // Similar calculators for the iterative algorithm
0346   std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > recoScaleCalcJet;
0347   std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > recoScaleRatioCalcJet;
0348   std::unique_ptr<fftjet::Functor1<double, RecoFFTJet> > memberFactorCalcJet;
0349
0350   // Calculator for the jet distance used to estimate convergence
0351   // of the iterative algorithm
0352   std::unique_ptr<fftjet::Functor2<double, RecoFFTJet, RecoFFTJet> > jetDistanceCalc;
0353
0354   // Vector of selected tree nodes
0355   std::vector<SparseTree::NodeId> nodes;
0356
0357   // Vector of selected preclusters
0358   std::vector<fftjet::Peak> preclusters;
0359
0360   // Vector of reconstructed jets (we will refill it in every event)
0361   std::vector<RecoFFTJet> recoJets;
0362
0363   // Cluster occupancy calculated as a function of level number
0364   std::vector<unsigned> occupancy;
0365
0366   // The thresholds obtained by the LOCALLY_ADAPTIVE method
0367   std::vector<double> thresholds;
0368
0369   // Minimum, maximum and used level calculated by some algorithms
0370   unsigned minLevel, maxLevel, usedLevel;
0371
0372   // Unclustered/unused energy produced during recombination
0373   fftjetcms::VectorLike unclustered;
0374   double unused;
0375
0376   // Quantities defined below are used in the iterative mode only
0377   std::vector<fftjet::Peak> iterPreclusters;
0378   std::vector<RecoFFTJet> iterJets;
0379   unsigned iterationsPerformed;
0380
0381   // Vectors of constituents
0382   std::vector<std::vector<reco::CandidatePtr> > constituents;
0383
0384   // Vector of pile-up. We will subtract it from the
0385   // 4-vectors of reconstructed jets.
0386   std::vector<fftjetcms::VectorLike> pileup;
0387
0388   // The pile-up transverse energy density discretization grid.
0389   // Note that this is _density_, not energy. To get energy,
0390   // multiply by cell area.
0391   std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> > pileupEnergyFlow;
0392
0393   // The functor that calculates the pile-up density
0394   std::unique_ptr<fftjetcms::AbsPileupCalculator> pileupDensityCalc;
0395
0396   // Memory buffers related to pile-up subtraction
0397   std::vector<fftjet::AbsKernel2d*> memFcns2dVec;
0398   std::vector<double> doubleBuf;
0399   std::vector<unsigned> cellCountsVec;
0400
0401   // Tokens for data access
0402   edm::EDGetTokenT<reco::PattRecoTree<double, reco::PattRecoPeak<double> > > input_recotree_token_d_;
0403   edm::EDGetTokenT<reco::PattRecoTree<float, reco::PattRecoPeak<float> > > input_recotree_token_f_;
0404   edm::EDGetTokenT<std::vector<reco::FFTAnyJet<reco::GenJet> > > input_genjet_token_;
0405   edm::EDGetTokenT<reco::DiscretizedEnergyFlow> input_energyflow_token_;
0406   edm::EDGetTokenT<reco::FFTJetPileupSummary> input_pusummary_token_;
0407
0408   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometry_token_;
0409   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> topology_token_;
0410 };
0411
0412 #endif  // RecoJets_FFTJetProducers_FFTJetProducer_h