FastJet 3.0.2
JetMedianBackgroundEstimator.hh
00001 #ifndef __FASTJET_BACKGROUND_ESTIMATOR_HH__
00002 #define __FASTJET_BACKGROUND_ESTIMATOR_HH__
00003 
00004 //STARTHEADER
00005 // $Id: JetMedianBackgroundEstimator.hh 2689 2011-11-14 14:51:06Z soyez $
00006 //
00007 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00008 //
00009 //----------------------------------------------------------------------
00010 // This file is part of FastJet.
00011 //
00012 //  FastJet is free software; you can redistribute it and/or modify
00013 //  it under the terms of the GNU General Public License as published by
00014 //  the Free Software Foundation; either version 2 of the License, or
00015 //  (at your option) any later version.
00016 //
00017 //  The algorithms that underlie FastJet have required considerable
00018 //  development and are described in hep-ph/0512210. If you use
00019 //  FastJet as part of work towards a scientific publication, please
00020 //  include a citation to the FastJet paper.
00021 //
00022 //  FastJet is distributed in the hope that it will be useful,
00023 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00024 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00025 //  GNU General Public License for more details.
00026 //
00027 //  You should have received a copy of the GNU General Public License
00028 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00029 //----------------------------------------------------------------------
00030 //ENDHEADER
00031 
00032 #include <fastjet/ClusterSequenceAreaBase.hh>
00033 #include <fastjet/AreaDefinition.hh>
00034 #include <fastjet/FunctionOfPseudoJet.hh>
00035 #include <fastjet/Selector.hh>
00036 #include <fastjet/tools/BackgroundEstimatorBase.hh>
00037 #include <iostream>
00038 
00039 FASTJET_BEGIN_NAMESPACE     // defined in fastjet/internal/base.hh
00040 
00041 
00042 /// @ingroup tools_background
00043 /// \class JetMedianBackgroundEstimator
00044 ///
00045 /// Class to estimate the pt density of the background per unit area,
00046 /// using the median of the distribution of pt/area from jets that
00047 /// pass some selection criterion.
00048 ///
00049 /// Events are passed either in the form of the event particles (in
00050 /// which they're clustered by the class), a ClusterSequenceArea (in
00051 /// which case the jets used are those returned by "inclusive_jets()")
00052 /// or directly as a set of jets.
00053 ///
00054 /// The selection criterion is typically a geometrical one (e.g. all
00055 /// jets with |y|<2) sometimes supplemented with some kinematical
00056 /// restriction (e.g. exclusion of the two hardest jets). It is passed
00057 /// to the class through a Selector.
00058 ///
00059 /// Beware: 
00060 ///   by default, to correctly handle partially empty events, the
00061 ///   class attempts to calculate an "empty area", based
00062 ///   (schematically) on
00063 ///
00064 ///          range.total_area() - sum_{jets_in_range} jets.area()
00065 ///  
00066 ///   For ranges with small areas, this can be inaccurate (particularly 
00067 ///   relevant in dense events where empty_area should be zero and ends
00068 ///   up not being zero).
00069 ///
00070 ///   This calculation of empty area can be avoided if a
00071 ///   ClusterSequenceArea class with explicit ghosts
00072 ///   (ActiveAreaExplicitGhosts) is used.  This is _recommended_
00073 ///   unless speed requirements cause you to use Voronoi areas. For
00074 ///   speedy background estimation you could also consider using
00075 ///   GridMedianBackgroundEstimator.
00076 ///
00077 ///
00078 class JetMedianBackgroundEstimator : public BackgroundEstimatorBase {
00079 public:
00080   /// @name constructors and destructors
00081   //\{
00082   //----------------------------------------------------------------
00083   /// Constructor that sets the rho range as well as the jet
00084   /// definition and area definition to be used to cluster the
00085   /// particles. Prior to the estimation of rho, one has to provide
00086   /// the particles to cluster using set_particles(...)
00087   ///
00088   /// \param rho_range  the Selector specifying which jets will be considered
00089   /// \param jet_def    the jet definition to use for the clustering
00090   /// \param area_def   the area definition to use for the clustering
00091   JetMedianBackgroundEstimator(const Selector &rho_range,
00092                                const JetDefinition &jet_def,
00093                                const AreaDefinition &area_def);
00094 
00095   /// ctor from a ClusterSequenceAreaBase with area
00096   ///
00097   /// \param rho_range   the Selector specifying which jets will be considered
00098   /// \param csa         the ClusterSequenceArea to use
00099   ///
00100   /// Pre-conditions: 
00101   ///  - one should be able to estimate the "empty area" (i.e. the area
00102   ///    not occupied by jets). This is feasible if at least one of the following
00103   ///    conditions is satisfied:
00104   ///     ( i) the ClusterSequence has explicit ghosts
00105   ///     (ii) the range has a computable area.
00106   ///  - the jet algorithm must be suited for median computation
00107   ///    (otherwise a warning will be issues)
00108   ///
00109   /// Note that selectors with e.g. hardest-jets exclusion do not have
00110   /// a well-defined area. For this reasons, it is STRONGLY advised to
00111   /// use an area with explicit ghosts.
00112   JetMedianBackgroundEstimator(const Selector &rho_range,
00113                                const ClusterSequenceAreaBase &csa);
00114 
00115 
00116   /// Default constructor that optionally sets the rho range. The
00117   /// configuration must be done later calling
00118   /// set_cluster_sequence(...) or set_jets(...).
00119   ///
00120   /// \param rho_range   the Selector specifying which jets will be considered
00121   ///
00122   JetMedianBackgroundEstimator(const Selector &rho_range = SelectorIdentity())
00123     : _rho_range(rho_range), _jet_def(JetDefinition()) { reset(); }
00124   
00125 
00126   /// default dtor
00127   ~JetMedianBackgroundEstimator(){}
00128 
00129   //\}
00130 
00131 
00132   /// @name setting a new event
00133   //\{
00134   //----------------------------------------------------------------
00135 
00136   /// tell the background estimator that it has a new event, composed
00137   /// of the specified particles.
00138   virtual void set_particles(const std::vector<PseudoJet> & particles);
00139 
00140   /// (re)set the cluster sequence (with area support) to be used by
00141   /// future calls to rho() etc. 
00142   ///
00143   /// \param csa  the cluster sequence area
00144   ///
00145   /// Pre-conditions: 
00146   ///  - one should be able to estimate the "empty area" (i.e. the area
00147   ///    not occupied by jets). This is feasible if at least one of the following
00148   ///    conditions is satisfied:
00149   ///     ( i) the ClusterSequence has explicit ghosts
00150   ///     (ii) the range selected has a computable area.
00151   ///  - the jet algorithm must be suited for median computation
00152   ///    (otherwise a warning will be issues)
00153   ///
00154   /// Note that selectors with e.g. hardest-jets exclusion do not have
00155   /// a well-defined area. For this reasons, it is STRONGLY advised to
00156   /// use an area with explicit ghosts.
00157   void set_cluster_sequence(const ClusterSequenceAreaBase & csa);
00158 
00159   /// (re)set the jets (which must have area support) to be used by future
00160   /// calls to rho() etc.; for the conditions that must be satisfied
00161   /// by the jets, see the Constructor that takes jets.
00162   void set_jets(const std::vector<PseudoJet> &jets);
00163 
00164   /// (re)set the selector to be used for future calls to rho() etc.
00165   void set_selector(const Selector & rho_range_selector) {
00166     _rho_range = rho_range_selector;
00167     _uptodate = false;
00168   }
00169 
00170   //\}
00171 
00172 
00173   /// @name  retrieving fundamental information
00174   //\{
00175   //----------------------------------------------------------------
00176 
00177   /// get rho, the median background density per unit area
00178   double rho() const;
00179 
00180   /// get sigma, the background fluctuations per unit area
00181   double sigma() const;
00182 
00183   /// get rho, the median background density per unit area, locally at
00184   /// the position of a given jet.
00185   ///
00186   /// If the Selector associated with the range takes a reference jet
00187   /// (i.e. is relocatable), then for subsequent operations the
00188   /// Selector has that jet set as its reference.
00189   double rho(const PseudoJet & jet);
00190 
00191   /// get sigma, the background fluctuations per unit area,
00192   /// locally at the position of a given jet.
00193   ///
00194   /// If the Selector associated with the range takes a reference jet
00195   /// (i.e. is relocatable), then for subsequent operations the
00196   /// Selector has that jet set as its reference.
00197   double sigma(const PseudoJet &jet);
00198 
00199   /// returns true if this background estimator has support for
00200   /// determination of sigma
00201   virtual bool has_sigma() {return true;}
00202 
00203   //\}
00204   
00205   /// @name  retrieving additional useful information
00206   //\{
00207   //----------------------------------------------------------------
00208   /// Returns the mean area of the jets used to actually compute the
00209   /// background properties in the last call of rho() or sigma()
00210   double mean_area() const{
00211     _recompute_if_needed();
00212     return _mean_area;
00213   }
00214   
00215   /// returns the number of jets used to actually compute the
00216   /// background properties in the last call of rho() or sigma()
00217   unsigned int n_jets_used() const{
00218     _recompute_if_needed();
00219     return _n_jets_used;
00220   }
00221 
00222   /// Returns the estimate of the area (within the range defined by
00223   /// the selector) that is not occupied by jets. The value is that
00224   /// for the last call of rho() or sigma()
00225   ///
00226   /// The answer is defined to be zero if the area calculation
00227   /// involved explicit ghosts; if the area calculation was an active
00228   /// area, then use is made of the active area's internal list of
00229   /// pure ghost jets (taking those that pass the selector); otherwise
00230   /// it is based on the difference between the selector's total area
00231   /// and the area of the jets that pass the selector.
00232   ///
00233   /// The result here is just the cached result of the corresponding
00234   /// call to the ClusterSequenceAreaBase function.
00235   double empty_area() const{
00236     _recompute_if_needed();
00237     return _empty_area;
00238   }
00239 
00240   /// Returns the number of empty jets used when computing the
00241   /// background properties. The value is that for the last call of
00242   /// rho() or sigma().
00243   ///
00244   /// If the area has explicit ghosts the result is zero; for active
00245   /// areas it is the number of internal pure ghost jets that pass the
00246   /// selector; otherwise it is deduced from the empty area, divided by 
00247   /// \f$ 0.55 \pi R^2 \f$ (the average pure-ghost-jet area).
00248   ///
00249   /// The result here is just the cached result of the corresponding
00250   /// call to the ClusterSequenceAreaBase function.
00251   double n_empty_jets() const{
00252     _recompute_if_needed();
00253     return _n_empty_jets;
00254   }
00255 
00256   //}
00257 
00258 
00259   /// @name configuring behaviour
00260   //\{
00261   //----------------------------------------------------------------
00262 
00263   /// Resets the class to its default state, including the choice to
00264   /// use 4-vector areas.
00265   ///
00266   void reset();
00267 
00268   /// By default when calculating pt/Area for a jet, it is the
00269   /// transverse component of the 4-vector area that is used in the ratiof \f$p_t/A\f$. 
00270   /// Calling this function with a "false" argument causes the scalar area to
00271   /// be used instead. 
00272   ///
00273   /// While the difference between the two choices is usually small,
00274   /// for high-precision work it is usually the 4-vector area that is
00275   /// to be preferred.
00276   ///
00277   ///  \param use_it             whether one uses the 4-vector area or not (true by default)
00278   void set_use_area_4vector(bool use_it = true){
00279     _use_area_4vector = use_it;
00280     _uptodate = false;
00281   }  
00282 
00283   /// check if the estimator uses the 4-vector area or the scalar area
00284   bool use_area_4vector() const{ return _use_area_4vector;}
00285 
00286   /// The FastJet v2.X sigma calculation had a small spurious offset
00287   /// in the limit of a small number of jets. This is fixed by default
00288   /// in versions 3 upwards. The old behaviour can be obtained with a
00289   /// call to this function.
00290   void set_provide_fj2_sigma(bool provide_fj2_sigma = true) {
00291     _provide_fj2_sigma = provide_fj2_sigma;
00292     _uptodate = false;
00293   }
00294 
00295   /// Set a pointer to a class that calculates the quantity whose
00296   /// median will be calculated; if the pointer is null then pt/area
00297   /// is used (as occurs also if this function is not called).
00298   ///
00299   /// Note that this is still <i>preliminary</i> in FastJet 3.0 and
00300   /// that backward compatibility is not guaranteed in future releases
00301   /// of FastJet
00302   void set_jet_density_class(const FunctionOfPseudoJet<double> * jet_density_class);
00303 
00304   /// return the pointer to the jet density class
00305   const FunctionOfPseudoJet<double> *  jet_density_class() const{
00306     return _jet_density_class;
00307   }
00308 
00309   /// Set a pointer to a class that calculates the rescaling factor as
00310   /// a function of the jet (position). Note that the rescaling factor
00311   /// is used both in the determination of the "global" rho (the pt/A
00312   /// of each jet is divided by this factor) and when asking for a
00313   /// local rho (the result is multiplied by this factor).
00314   ///
00315   /// The BackgroundRescalingYPolynomial class can be used to get a
00316   /// rescaling that depends just on rapidity.
00317   virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) {
00318     BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
00319     _uptodate = false;
00320   }
00321 
00322   //\}
00323 
00324   /// @name description
00325   //\{
00326   //----------------------------------------------------------------
00327 
00328   /// returns a textual description of the background estimator
00329   std::string description() const;
00330 
00331   //\}
00332 
00333 
00334 private:
00335 
00336   /// do the actual job
00337   void _compute() const;
00338    
00339   /// check if the properties need to be recomputed 
00340   /// and do so if needed
00341   void _recompute_if_needed() const {
00342     if (!_uptodate) _compute();
00343     _uptodate = true;
00344   }
00345 
00346   /// for estimation using a selector that takes a reference jet
00347   /// (i.e. a selector that can be relocated) this function allows one
00348   /// to set its position.
00349   ///
00350   /// Note that this HAS to be called before any attempt to compute
00351   /// the background properties. The call is, however, performed
00352   /// automatically by the functions rho(jet) and sigma(jet).
00353   void _recompute_if_needed(const PseudoJet &jet);
00354 
00355   /// check that the underlying structure is still alive
00356   /// throw an error otherwise
00357   void _check_csa_alive() const;
00358 
00359   /// check that the algorithm used for the clustering is adapted for
00360   /// background estimation (i.e. either kt or C/A)
00361   /// Issue a warning otherwise
00362   void _check_jet_alg_good_for_median() const;
00363   
00364   // the basic parameters of this class (passed through the variou ctors)
00365   Selector _rho_range;                   ///< range to compute the background in
00366   JetDefinition _jet_def;                ///< the jet def to use for teh clustering
00367   AreaDefinition _area_def;              ///< the area def to use for teh clustering
00368   std::vector<PseudoJet> _included_jets; ///< jets to be used
00369   
00370   // the tunable aprameters of the class
00371   bool _use_area_4vector;
00372   bool _provide_fj2_sigma;
00373   const FunctionOfPseudoJet<double> * _jet_density_class;
00374   //SharedPtr<BackgroundRescalingBase> _rescaling_class_sharedptr;
00375   
00376   // the actual results of the computation
00377   mutable double _rho;               ///< background estimated density per unit area
00378   mutable double _sigma;             ///< background estimated fluctuations
00379   mutable double _mean_area;         ///< mean area of the jets used to estimate the background
00380   mutable unsigned int _n_jets_used; ///< number of jets used to estimate the background
00381   mutable double _n_empty_jets;      ///< number of empty (pure-ghost) jets
00382   mutable double _empty_area;        ///< the empty (pure-ghost/unclustered) area!
00383 
00384   // internal variables
00385   SharedPtr<PseudoJetStructureBase> _csi; ///< allows to check if _csa is still valid
00386   PseudoJet _current_reference;           ///< current reference jet
00387   mutable bool _uptodate;                 ///< true when the background computation is up-to-date
00388 
00389   /// handle warning messages
00390   static LimitedWarning _warnings;
00391   static LimitedWarning _warnings_zero_area;
00392   static LimitedWarning _warnings_preliminary;
00393 };
00394 
00395 
00396 
00397 
00398 //----------------------------------------------------------------------
00399 /// @ingroup tools_background
00400 /// \class BackgroundJetPtDensity
00401 /// Class that implements pt/area_4vector.perp() for background estimation
00402 /// <i>(this is a preliminary class)</i>.
00403 class BackgroundJetPtDensity : public FunctionOfPseudoJet<double> {
00404 public:
00405   virtual double result(const PseudoJet & jet) const {
00406     return jet.perp() / jet.area_4vector().perp();
00407   }
00408   virtual std::string description() const {return "BackgroundJetPtDensity";}
00409 };
00410 
00411 
00412 //----------------------------------------------------------------------
00413 /// @ingroup tools_background
00414 /// \class BackgroundJetScalarPtDensity
00415 /// Class that implements (scalar pt sum of jet)/(scalar area of jet)
00416 /// for background estimation <i>(this is a preliminary class)</i>.
00417 ///
00418 /// Optionally it can return a quantity based on the sum of pt^n,
00419 /// e.g. for use in subtracting fragementation function moments.
00420 class BackgroundJetScalarPtDensity : public FunctionOfPseudoJet<double> {
00421 public:
00422   /// Default constructor provides background estimation with scalar pt sum
00423   BackgroundJetScalarPtDensity() : _pt_power(1) {}
00424 
00425   /// Constructor to provide background estimation based on 
00426   /// \f$ sum_{i\in jet} p_{ti}^{n} \f$
00427   BackgroundJetScalarPtDensity(double n) : _pt_power(n) {}
00428 
00429   virtual double result(const PseudoJet & jet) const;
00430 
00431   virtual std::string description() const {return "BackgroundScalarJetPtDensity";}
00432 
00433 private:
00434   double _pt_power;
00435 };
00436 
00437 //----------------------------------------------------------------------
00438 /// @ingroup tools_background
00439 /// \class BackgroundJetPtMDensity
00440 /// Class that implements
00441 /// \f$  \frac{1}{A} \sum_{i \in jet} (\sqrt{p_{ti}^2+m^2} - p_{ti}) \f$
00442 /// for background estimation <i>(this is a preliminary class)</i>.
00443 /// 
00444 ///
00445 /// This is useful for correcting jet masses in cases where the event
00446 /// involves massive particles.
00447 class BackgroundJetPtMDensity : public FunctionOfPseudoJet<double> {
00448 public:
00449   virtual double result(const PseudoJet & jet) const {
00450     std::vector<PseudoJet> constituents = jet.constituents();
00451     double scalar_ptm = 0;
00452     for (unsigned i = 0; i < constituents.size(); i++) {
00453       scalar_ptm += constituents[i].mperp() - constituents[i].perp();
00454     }
00455     return scalar_ptm / jet.area();
00456   }
00457 
00458   virtual std::string description() const {return "BackgroundPtMDensity";}
00459 };
00460 
00461 
00462 
00463 FASTJET_END_NAMESPACE
00464 
00465 #endif  // __BACKGROUND_ESTIMATOR_HH__
00466 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends