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