FastJet 3.0.2
|
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