JETAN/fastjet/fastjet/ClusterSequenceAreaBase.hh

   1 //STARTHEADER
   2 // $Id: ClusterSequenceAreaBase.hh 1503 2009-04-06 11:32:52Z salam $
   3 //
   4 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
   5 //
   6 //----------------------------------------------------------------------
   7 // This file is part of FastJet.
   8 //
   9 //  FastJet is free software; you can redistribute it and/or modify
  10 //  it under the terms of the GNU General Public License as published by
  11 //  the Free Software Foundation; either version 2 of the License, or
  12 //  (at your option) any later version.
  13 //
  14 //  The algorithms that underlie FastJet have required considerable
  15 //  development and are described in hep-ph/0512210. If you use
  16 //  FastJet as part of work towards a scientific publication, please
  17 //  include a citation to the FastJet paper.
  18 //
  19 //  FastJet is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 //  GNU General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU General Public License
  25 //  along with FastJet; if not, write to the Free Software
  26 //  Foundation, Inc.:
  27 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  28 //----------------------------------------------------------------------
  29 //ENDHEADER
  30
  31 #ifndef __FASTJET_CLUSTERSEQUENCEAREABASE_HH__
  32 #define __FASTJET_CLUSTERSEQUENCEAREABASE_HH__
  33
  34 #include "fastjet/ClusterSequence.hh"
  35 #include "fastjet/internal/LimitedWarning.hh"
  36 #include "fastjet/RangeDefinition.hh"
  37
  38 namespace fastjet {
  39
  40 /// base class that sets interface for extensions of ClusterSequence
  41 /// that provide information about the area of each jet;
  42 ///
  43 /// the virtual functions here all return 0, since no area determination
  44 /// is implemented.
  45 class ClusterSequenceAreaBase : public ClusterSequence {
  46 public:
  47
  48   /// a constructor which just carries out the construction of the
  49   /// parent class
  50   template<class L> ClusterSequenceAreaBase
  51          (const std::vector<L> & pseudojets,
  52           const JetDefinition & jet_def,
  53           const bool & writeout_combinations = false) :
  54            ClusterSequence(pseudojets, jet_def, writeout_combinations) {}
  55
  56
  57   /// default constructor
  58   ClusterSequenceAreaBase() {}
  59
  60
  61   /// destructor
  62   virtual ~ClusterSequenceAreaBase() {}
  63
  64
  65   /// return the area associated with the given jet; this base class
  66   /// returns 0.
  67   virtual double area       (const PseudoJet & jet) const {return 0.0;}
  68
  69   /// return the error (uncertainty) associated with the determination
  70   /// of the area of this jet; this base class returns 0.
  71   virtual double area_error (const PseudoJet & jet) const {return 0.0;}
  72
  73   /// return a PseudoJet whose 4-vector is defined by the following integral
  74   ///
  75   ///       \int drap d\phi PseudoJet("rap,phi,pt=one") *
  76   ///                           * Theta("rap,phi inside jet boundary")
  77   ///
  78   /// where PseudoJet("rap,phi,pt=one") is a 4-vector with the given
  79   /// rapidity (rap), azimuth (phi) and pt=1, while Theta("rap,phi
  80   /// inside jet boundary") is a function that is 1 when rap,phi
  81   /// define a direction inside the jet boundary and 0 otherwise.
  82   ///
  83   /// This base class returns a null 4-vector.
  84   virtual PseudoJet area_4vector(const PseudoJet & jet) const {
  85     return PseudoJet(0.0,0.0,0.0,0.0);}
  86
  87   /// true if a jet is made exclusively of ghosts
  88   ///
  89   /// NB: most area classes do not give any explicit ghost jets, but
  90   /// some do, and they should replace this function with their own
  91   /// version.
  92   virtual bool is_pure_ghost(const PseudoJet & jet) const {
  93     return false;
  94   }
  95
  96   /// returns true if ghosts are explicitly included within
  97   /// jets for this ClusterSequence;
  98   ///
  99   /// Derived classes that do include explicit ghosts should provide
 100   /// an alternative version of this routine and set it properly.
 101   virtual bool has_explicit_ghosts() const {
 102     return false;
 103   }
 104
 105   /// return the total area, within range, that is free of jets, in
 106   /// general based on the inclusive jets
 107   virtual double empty_area(const RangeDefinition & range) const;
 108
 109   /// return the total area, within range, that is free of jets, based
 110   /// on the supplied all_jets
 111   double empty_area_from_jets(const std::vector<PseudoJet> & all_jets,
 112                               const RangeDefinition & range) const;
 113
 114   /// return something similar to the number of pure ghost jets
 115   /// in the given range in an active area case.
 116   /// For the local implementation we return empty_area/(0.55 pi R^2),
 117   /// based on measured properties of ghost jets with kt and cam. Note
 118   /// that the number returned is a double.
 119   virtual double n_empty_jets(const RangeDefinition & range) const {
 120     double R = jet_def().R();
 121     return empty_area(range)/(0.55*pi*R*R);
 122   }
 123
 124   /// the median of (pt/area) for jets contained within range,
 125   /// making use also of the info on n_empty_jets
 126   double median_pt_per_unit_area(const RangeDefinition & range) const;
 127
 128   /// the median of (pt/area_4vector) for jets contained within
 129   /// making use also of the info on n_empty_jets
 130   double median_pt_per_unit_area_4vector(const RangeDefinition & range) const;
 131
 132   /// the function that does the work for median_pt_per_unit_area and
 133   /// median_pt_per_unit_area_4vector:
 134   /// - something_is_area_4vect = false -> use plain area
 135   /// - something_is_area_4vect = true  -> use 4-vector area
 136   double median_pt_per_unit_something(
 137                     const RangeDefinition & range, bool use_area_4vector) const;
 138
 139   /// using jets withing range (and with 4-vector areas if
 140   /// use_area_4vector), calculate the median pt/area, as well as an
 141   /// "error" (uncertainty), which is defined as the 1-sigma
 142   /// half-width of the distribution of pt/A, obtained by looking for
 143   /// the point below which we have (1-0.6827)/2 of the jets
 144   /// (including empty jets).
 145   ///
 146   /// The subtraction for a jet with uncorrected pt pt^U and area A is
 147   ///
 148   ///   pt^S = pt^U - median*A +- sigma*sqrt(A)
 149   ///
 150   /// where the error is only that associated with the fluctuations
 151   /// in the noise and not that associated with the noise having
 152   /// caused changes in the hard-particle content of the jet.
 153   ///
 154   /// NB: subtraction may also be done with 4-vector area of course,
 155   /// and this is recommended for jets with larger values of R, as
 156   /// long as rho has also been determined with a 4-vector area;
 157   /// using a scalar area causes one to neglect terms of relative
 158   /// order $R^2/8$ in the jet $p_t$.
 159   virtual void get_median_rho_and_sigma(const RangeDefinition & range,
 160                                         bool use_area_4vector,
 161                                         double & median, double & sigma,
 162                                         double & mean_area) const;
 163
 164   /// a more advanced version of get_median_rho_and_sigma, which allows
 165   /// one to use any "view" of the event containing all jets (so that,
 166   /// e.g. one might use Cam on a different resolution scale without
 167   /// have to rerun the algorithm).
 168   ///
 169   /// By default it will assume that "all" are not inclusive jets,
 170   /// so that in dealing with empty area it has to calculate
 171   /// the number of empty jets based on the empty area and the
 172   /// the observed <area> of jets rather than a surmised area
 173   ///
 174   /// Note that for small effective radii, this can cause problems
 175   /// because the harder jets get an area >> <ghost-jet-area>
 176   /// and so the estimate comes out all wrong. In these situations
 177   /// it is highly advisable to use an area with explicit ghosts, since
 178   /// then the "empty" jets are actually visible.
 179   virtual void get_median_rho_and_sigma(const std::vector<PseudoJet> & all_jets,
 180                                         const RangeDefinition & range,
 181                                         bool use_area_4vector,
 182                                         double & median, double & sigma,
 183                                         double & mean_area,
 184                                         bool all_are_inclusive = false) const;
 185
 186   /// same as the full version of get_median_rho_and_error, but without
 187   /// access to the mean_area
 188   virtual void get_median_rho_and_sigma(const RangeDefinition & range,
 189                                 bool use_area_4vector,
 190                                 double & median, double & sigma) const {
 191     double mean_area;
 192     get_median_rho_and_sigma(range,  use_area_4vector,
 193                              median,  sigma, mean_area);
 194   }
 195
 196
 197   /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range "range".
 198   /// exclude_above allows one to exclude large values of pt/area from fit.
 199   ///               (if negative, the cut is discarded)
 200   /// use_area_4vector = true uses the 4vector areas.
 201   virtual void parabolic_pt_per_unit_area(double & a, double & b,
 202                                           const RangeDefinition & range,
 203                                           double exclude_above=-1.0,
 204                                           bool use_area_4vector=false) const;
 205
 206   /// return a vector of all subtracted jets, using area_4vector, given rho.
 207   /// Only inclusive_jets above ptmin are subtracted and returned.
 208   /// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
 209   /// i.e. not necessarily ordered in pt once subtracted
 210   std::vector<PseudoJet> subtracted_jets(const double rho,
 211                                          const double ptmin=0.0) const;
 212
 213   /// return a vector of subtracted jets, using area_4vector.
 214   /// Only inclusive_jets above ptmin are subtracted and returned.
 215   /// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
 216   /// i.e. not necessarily ordered in pt once subtracted
 217   std::vector<PseudoJet> subtracted_jets(const RangeDefinition & range,
 218                                          const double ptmin=0.0) const;
 219
 220   /// return a subtracted jet, using area_4vector, given rho
 221   PseudoJet subtracted_jet(const PseudoJet & jet,
 222                            const double rho) const;
 223
 224   /// return a subtracted jet, using area_4vector; note
 225   /// that this is potentially inefficient if repeatedly used for many
 226   /// different jets, because rho will be recalculated each time
 227   /// around.
 228   PseudoJet subtracted_jet(const PseudoJet & jet,
 229                            const RangeDefinition & range) const;
 230
 231   /// return the subtracted pt, given rho
 232   double subtracted_pt(const PseudoJet & jet,
 233                        const double rho,
 234                        bool use_area_4vector=false) const;
 235
 236   /// return the subtracted pt; note that this is
 237   /// potentially inefficient if repeatedly used for many different
 238   /// jets, because rho will be recalculated each time around.
 239   double subtracted_pt(const PseudoJet & jet,
 240                        const RangeDefinition & range,
 241                        bool use_area_4vector=false) const;
 242
 243
 244 private:
 245   /// handle warning messages
 246   static LimitedWarning _warnings;
 247
 248   /// check the jet algorithm is suitable (and if not issue a warning)
 249   void _check_jet_alg_good_for_median() const;
 250
 251 };
 252
 253
 254
 255 } // fastjet namespace
 256
 257 #endif // __FASTJET_CLUSTERSEQUENCEAREABASE_HH__