fastjet/contrib/JetCleanser.hh

0001 // JetCleanser Package
0002 // Questions/Comments? dkrohn@physics.harvard.edu mattlow@uchicago.edu schwartz@physics.harvard.edu liantaow@uchicago.edu
0003 //
0004 // Copyright (c) 2013
0005 // David Krohn, Matthew Low, Matthew Schwartz, and Lian-Tao Wang
0006 //
0007 // $Id$
0008 //
0009 //----------------------------------------------------------------------
0010 // This file is part of FastJet contrib.
0011 //
0012 // It is free software; you can redistribute it and/or modify it under
0013 // the terms of the GNU General Public License as published by the
0014 // Free Software Foundation; either version 2 of the License, or (at
0015 // your option) any later version.
0016 //
0017 // It is distributed in the hope that it will be useful, but WITHOUT
0018 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0019 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
0020 // License for more details.
0021 //
0022 // You should have received a copy of the GNU General Public License
0023 // along with this code. If not, see <http://www.gnu.org/licenses/>.
0024 //----------------------------------------------------------------------
0025
0026 #ifndef __FASTJET_CONTRIB_JETCLEANSER_HH__
0027 #define __FASTJET_CONTRIB_JETCLEANSER_HH__
0028
0029 #include <fastjet/internal/base.hh>
0030
0031 #include "fastjet/ClusterSequence.hh"
0032 #include "fastjet/Error.hh"
0033 #include "fastjet/JetDefinition.hh"
0034 //#include "fastjet/FunctionOfPseudoJet.hh"
0035
0036 #include <map>
0037 #include <sstream>
0038 #include <string>
0039
0040 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0041
0042 namespace contrib{
0043
0044 //------------------------------------------------------------------------
0045 /// \class JetCleanser
0046 /// This class implements jet cleansing which is a substructure technique
0047 /// designed to correct for the effects of pileup.
0048 ///
0049 /// Jet cleansing takes as input either:
0050 ///  (A) "calorimeter" or "particle" four-vectors which are charged and neutral
0051 ///       particles all grouped together, charged tracks coming from the primary
0052 ///       interaction, and charged tracks coming from pileup.
0053 ///  (B) neutral particles, charged tracks coming from the primary
0054 ///       interaction, and charged tracks coming from pileup.
0055 ///
0056 /// The output is a jet for which the constituents themselves have been
0057 /// corrected so that the output can be used to calculate both
0058 /// kinematic quantities and shape variables.  The rescale subjets are also
0059 /// available via the pieces() function.
0060 ///
0061 class JetCleanser {
0062
0063 public:
0064   enum cleansing_mode {
0065     jvf_cleansing,      /// gamma0 = gamma1
0066     linear_cleansing,   /// approximate gamma0 as constant
0067     gaussian_cleansing  /// approximate gamma's as truncated gaussians and maximize likelihood
0068   };
0069
0070   enum input_mode {
0071     input_nc_together,  /// input is all charged and neutrals together (as particles or calo cells),
0072                         ///          charged primary interaction particles, and charged pileup particles
0073     input_nc_separate   /// input is all neutral particles, charged primary interaction particles,
0074                         ///          and charged pileup particles
0075   };
0076
0077 public:
0078   // constructors
0079   JetCleanser(JetDefinition subjet_def, cleansing_mode cmode, input_mode imode);
0080   JetCleanser(double rsub, cleansing_mode cmode, input_mode imode);
0081
0082   // destructor
0083   ~JetCleanser(){}
0084
0085   // settings
0086   void SetGroomingParameters(double fcut, int nsjmin);
0087   inline void SetTrimming(double fcut) { SetGroomingParameters(fcut,0); }
0088   inline void SetFiltering(int nsj) { SetGroomingParameters(1.0,nsj); }
0089   void SetLinearParameters(double g0_mean=0.67);
0090   void SetGaussianParameters(double g0_mean=0.67, double g1_mean=0.67, double g0_width=0.15, double g1_width=0.25);
0091
0092   // standard usage
0093   std::string description() const;
0094   PseudoJet operator()(const PseudoJet & jet,                             // For use with input_nc_together mode,
0095                        const std::vector<PseudoJet> & tracks_lv,          // it takes a plain jet, charged tracks from
0096                        const std::vector<PseudoJet> & tracks_pu) const;   // the leading vertex, and charged tracks
0097                                                                           // from pileup.
0098   PseudoJet operator()(const std::vector<PseudoJet> & neutrals_all,       // For use with input_nc_separate mode,
0099                        const std::vector<PseudoJet> & tracks_lv,          // it takes all neutral particles, charged
0100                        const std::vector<PseudoJet> & tracks_pu) const;   // tracks from the leading vertex, and charged
0101                                                                           // tracks from pileup.
0102   void   _RunTests();
0103
0104
0105 private:
0106   // Modification to satisfy C++11 (thanks to Gavin Salam)
0107   //static const double jc_zero = 1.0e-6;
0108   static const double jc_zero;
0109
0110   double _rsub;
0111   double _fcut;
0112   double _nsjmin;
0113   JetDefinition _subjet_def;
0114
0115   cleansing_mode _cleansing_mode;
0116   input_mode _input_mode;
0117
0118   double _linear_gamma0_mean;
0119   double _gaussian_gamma0_mean;
0120   double _gaussian_gamma0_width;
0121   double _gaussian_gamma1_mean;
0122   double _gaussian_gamma1_width;
0123
0124   // helper functions
0125   void   _SetDefaults();
0126   void   _CheckRescalingValues(double & pt_all, const double & ptc_lv, const double & ptc_pu) const;
0127   double _GetSubjetRescaling_nctogether(double pt_all, double ptc_lv, double ptc_pu) const;
0128   double _GetSubjetRescaling_ncseparate(double ptn_all, double ptc_lv, double ptc_pu) const;
0129   double _GaussianGetMinimizedGamma0(double pt_all, double ptc_lv, double ptc_pu) const;
0130   double _GaussianGetGamma1(double gamma0, double pt_all, double ptc_lv, double ptc_pu) const;
0131   double _GaussianFunction(double x, void * params) const;
0132   void   _RunTestRescaling(double pt_all, double ptc_lv, double ptc_pu) const;
0133 };
0134
0135
0136 // helper function
0137 std::vector<PseudoJet> RescalePseudoJetVector(const std::vector<PseudoJet> & jets, const double s_factor);
0138
0139 // helper function
0140 PseudoJet RescalePseudoJetConstituents(const PseudoJet & jet, const double s_factor);
0141
0142 // helper function
0143 std::vector< std::vector<PseudoJet> > ClusterSets(const JetDefinition & jet_def,
0144                                                   const std::vector<PseudoJet> & cluster_set,
0145                                                   const std::vector< std::vector<PseudoJet> > & follow_sets,
0146                                                   const double &ptmin=0.0);
0147
0148
0149
0150 //------------------------------------------------------------------------
0151 /// \class FollowSetGhostInfo
0152 /// This class keeps tracks of which set of PseudoJets and the index where
0153 /// the associated ghost PseudoJet came from.
0154 ///
0155 class FollowSetGhostInfo : public PseudoJet::UserInfoBase {
0156   public:
0157     inline FollowSetGhostInfo(int set_id, int ind_id) { _set_id = set_id; _ind_id = ind_id; }
0158     inline int GetSetId() { return _set_id; }
0159     inline int GetIndId() { return _ind_id; }
0160
0161   private:
0162     int _set_id;
0163     int _ind_id;
0164 };
0165
0166 } // namespace contrib
0167
0168 FASTJET_END_NAMESPACE
0169
0170 #endif  // __FASTJET_CONTRIB_JETCLEANSER_HH__