FastJet  3.0.6
JetMedianBackgroundEstimator.hh
1 #ifndef __FASTJET_BACKGROUND_ESTIMATOR_HH__
2 #define __FASTJET_BACKGROUND_ESTIMATOR_HH__
3 
4 //STARTHEADER
5 // $Id: JetMedianBackgroundEstimator.hh 2689 2011-11-14 14:51:06Z soyez $
6 //
7 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8 //
9 //----------------------------------------------------------------------
10 // This file is part of FastJet.
11 //
12 // FastJet is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU General Public License as published by
14 // the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // The algorithms that underlie FastJet have required considerable
18 // development and are described in hep-ph/0512210. If you use
19 // FastJet as part of work towards a scientific publication, please
20 // include a citation to the FastJet paper.
21 //
22 // FastJet is distributed in the hope that it will be useful,
23 // but WITHOUT ANY WARRANTY; without even the implied warranty of
24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 // GNU General Public License for more details.
26 //
27 // You should have received a copy of the GNU General Public License
28 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
29 //----------------------------------------------------------------------
30 //ENDHEADER
31 
32 #include <fastjet/ClusterSequenceAreaBase.hh>
33 #include <fastjet/AreaDefinition.hh>
34 #include <fastjet/FunctionOfPseudoJet.hh>
35 #include <fastjet/Selector.hh>
36 #include <fastjet/tools/BackgroundEstimatorBase.hh>
37 #include <iostream>
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 
42 /// @ingroup tools_background
43 /// \class JetMedianBackgroundEstimator
44 ///
45 /// Class to estimate the pt density of the background per unit area,
46 /// using the median of the distribution of pt/area from jets that
47 /// pass some selection criterion.
48 ///
49 /// Events are passed either in the form of the event particles (in
50 /// which they're clustered by the class), a ClusterSequenceArea (in
51 /// which case the jets used are those returned by "inclusive_jets()")
52 /// or directly as a set of jets.
53 ///
54 /// The selection criterion is typically a geometrical one (e.g. all
55 /// jets with |y|<2) sometimes supplemented with some kinematical
56 /// restriction (e.g. exclusion of the two hardest jets). It is passed
57 /// to the class through a Selector.
58 ///
59 /// Beware:
60 /// by default, to correctly handle partially empty events, the
61 /// class attempts to calculate an "empty area", based
62 /// (schematically) on
63 ///
64 /// range.total_area() - sum_{jets_in_range} jets.area()
65 ///
66 /// For ranges with small areas, this can be inaccurate (particularly
67 /// relevant in dense events where empty_area should be zero and ends
68 /// up not being zero).
69 ///
70 /// This calculation of empty area can be avoided if a
71 /// ClusterSequenceArea class with explicit ghosts
72 /// (ActiveAreaExplicitGhosts) is used. This is _recommended_
73 /// unless speed requirements cause you to use Voronoi areas. For
74 /// speedy background estimation you could also consider using
75 /// GridMedianBackgroundEstimator.
76 ///
77 ///
79 public:
80  /// @name constructors and destructors
81  //\{
82  //----------------------------------------------------------------
83  /// Constructor that sets the rho range as well as the jet
84  /// definition and area definition to be used to cluster the
85  /// particles. Prior to the estimation of rho, one has to provide
86  /// the particles to cluster using set_particles(...)
87  ///
88  /// \param rho_range the Selector specifying which jets will be considered
89  /// \param jet_def the jet definition to use for the clustering
90  /// \param area_def the area definition to use for the clustering
91  JetMedianBackgroundEstimator(const Selector &rho_range,
92  const JetDefinition &jet_def,
93  const AreaDefinition &area_def);
94 
95  /// ctor from a ClusterSequenceAreaBase with area
96  ///
97  /// \param rho_range the Selector specifying which jets will be considered
98  /// \param csa the ClusterSequenceArea to use
99  ///
100  /// Pre-conditions:
101  /// - one should be able to estimate the "empty area" (i.e. the area
102  /// not occupied by jets). This is feasible if at least one of the following
103  /// conditions is satisfied:
104  /// ( i) the ClusterSequence has explicit ghosts
105  /// (ii) the range has a computable area.
106  /// - the jet algorithm must be suited for median computation
107  /// (otherwise a warning will be issues)
108  ///
109  /// Note that selectors with e.g. hardest-jets exclusion do not have
110  /// a well-defined area. For this reasons, it is STRONGLY advised to
111  /// use an area with explicit ghosts.
112  JetMedianBackgroundEstimator(const Selector &rho_range,
113  const ClusterSequenceAreaBase &csa);
114 
115 
116  /// Default constructor that optionally sets the rho range. The
117  /// configuration must be done later calling
118  /// set_cluster_sequence(...) or set_jets(...).
119  ///
120  /// \param rho_range the Selector specifying which jets will be considered
121  ///
122  JetMedianBackgroundEstimator(const Selector &rho_range = SelectorIdentity())
123  : _rho_range(rho_range), _jet_def(JetDefinition()) { reset(); }
124 
125 
126  /// default dtor
128 
129  //\}
130 
131 
132  /// @name setting a new event
133  //\{
134  //----------------------------------------------------------------
135 
136  /// tell the background estimator that it has a new event, composed
137  /// of the specified particles.
138  virtual void set_particles(const std::vector<PseudoJet> & particles);
139 
140  /// (re)set the cluster sequence (with area support) to be used by
141  /// future calls to rho() etc.
142  ///
143  /// \param csa the cluster sequence area
144  ///
145  /// Pre-conditions:
146  /// - one should be able to estimate the "empty area" (i.e. the area
147  /// not occupied by jets). This is feasible if at least one of the following
148  /// conditions is satisfied:
149  /// ( i) the ClusterSequence has explicit ghosts
150  /// (ii) the range selected has a computable area.
151  /// - the jet algorithm must be suited for median computation
152  /// (otherwise a warning will be issues)
153  ///
154  /// Note that selectors with e.g. hardest-jets exclusion do not have
155  /// a well-defined area. For this reasons, it is STRONGLY advised to
156  /// use an area with explicit ghosts.
157  void set_cluster_sequence(const ClusterSequenceAreaBase & csa);
158 
159  /// (re)set the jets (which must have area support) to be used by future
160  /// calls to rho() etc.; for the conditions that must be satisfied
161  /// by the jets, see the Constructor that takes jets.
162  void set_jets(const std::vector<PseudoJet> &jets);
163 
164  /// (re)set the selector to be used for future calls to rho() etc.
165  void set_selector(const Selector & rho_range_selector) {
166  _rho_range = rho_range_selector;
167  _uptodate = false;
168  }
169 
170  //\}
171 
172 
173  /// @name retrieving fundamental information
174  //\{
175  //----------------------------------------------------------------
176 
177  /// get rho, the median background density per unit area
178  double rho() const;
179 
180  /// get sigma, the background fluctuations per unit area
181  double sigma() const;
182 
183  /// get rho, the median background density per unit area, locally at
184  /// the position of a given jet.
185  ///
186  /// If the Selector associated with the range takes a reference jet
187  /// (i.e. is relocatable), then for subsequent operations the
188  /// Selector has that jet set as its reference.
189  double rho(const PseudoJet & jet);
190 
191  /// get sigma, the background fluctuations per unit area,
192  /// locally at the position of a given jet.
193  ///
194  /// If the Selector associated with the range takes a reference jet
195  /// (i.e. is relocatable), then for subsequent operations the
196  /// Selector has that jet set as its reference.
197  double sigma(const PseudoJet &jet);
198 
199  /// returns true if this background estimator has support for
200  /// determination of sigma
201  virtual bool has_sigma() {return true;}
202 
203  //\}
204 
205  /// @name retrieving additional useful information
206  //\{
207  //----------------------------------------------------------------
208  /// Returns the mean area of the jets used to actually compute the
209  /// background properties in the last call of rho() or sigma()
210  double mean_area() const{
211  _recompute_if_needed();
212  return _mean_area;
213  }
214 
215  /// returns the number of jets used to actually compute the
216  /// background properties in the last call of rho() or sigma()
217  unsigned int n_jets_used() const{
218  _recompute_if_needed();
219  return _n_jets_used;
220  }
221 
222  /// Returns the estimate of the area (within the range defined by
223  /// the selector) that is not occupied by jets. The value is that
224  /// for the last call of rho() or sigma()
225  ///
226  /// The answer is defined to be zero if the area calculation
227  /// involved explicit ghosts; if the area calculation was an active
228  /// area, then use is made of the active area's internal list of
229  /// pure ghost jets (taking those that pass the selector); otherwise
230  /// it is based on the difference between the selector's total area
231  /// and the area of the jets that pass the selector.
232  ///
233  /// The result here is just the cached result of the corresponding
234  /// call to the ClusterSequenceAreaBase function.
235  double empty_area() const{
236  _recompute_if_needed();
237  return _empty_area;
238  }
239 
240  /// Returns the number of empty jets used when computing the
241  /// background properties. The value is that for the last call of
242  /// rho() or sigma().
243  ///
244  /// If the area has explicit ghosts the result is zero; for active
245  /// areas it is the number of internal pure ghost jets that pass the
246  /// selector; otherwise it is deduced from the empty area, divided by
247  /// \f$ 0.55 \pi R^2 \f$ (the average pure-ghost-jet area).
248  ///
249  /// The result here is just the cached result of the corresponding
250  /// call to the ClusterSequenceAreaBase function.
251  double n_empty_jets() const{
252  _recompute_if_needed();
253  return _n_empty_jets;
254  }
255 
256  //}
257 
258 
259  /// @name configuring behaviour
260  //\{
261  //----------------------------------------------------------------
262 
263  /// Resets the class to its default state, including the choice to
264  /// use 4-vector areas.
265  ///
266  void reset();
267 
268  /// By default when calculating pt/Area for a jet, it is the
269  /// transverse component of the 4-vector area that is used in the ratiof \f$p_t/A\f$.
270  /// Calling this function with a "false" argument causes the scalar area to
271  /// be used instead.
272  ///
273  /// While the difference between the two choices is usually small,
274  /// for high-precision work it is usually the 4-vector area that is
275  /// to be preferred.
276  ///
277  /// \param use_it whether one uses the 4-vector area or not (true by default)
278  void set_use_area_4vector(bool use_it = true){
279  _use_area_4vector = use_it;
280  _uptodate = false;
281  }
282 
283  /// check if the estimator uses the 4-vector area or the scalar area
284  bool use_area_4vector() const{ return _use_area_4vector;}
285 
286  /// The FastJet v2.X sigma calculation had a small spurious offset
287  /// in the limit of a small number of jets. This is fixed by default
288  /// in versions 3 upwards. The old behaviour can be obtained with a
289  /// call to this function.
290  void set_provide_fj2_sigma(bool provide_fj2_sigma = true) {
291  _provide_fj2_sigma = provide_fj2_sigma;
292  _uptodate = false;
293  }
294 
295  /// Set a pointer to a class that calculates the quantity whose
296  /// median will be calculated; if the pointer is null then pt/area
297  /// is used (as occurs also if this function is not called).
298  ///
299  /// Note that this is still <i>preliminary</i> in FastJet 3.0 and
300  /// that backward compatibility is not guaranteed in future releases
301  /// of FastJet
302  void set_jet_density_class(const FunctionOfPseudoJet<double> * jet_density_class);
303 
304  /// return the pointer to the jet density class
306  return _jet_density_class;
307  }
308 
309  /// Set a pointer to a class that calculates the rescaling factor as
310  /// a function of the jet (position). Note that the rescaling factor
311  /// is used both in the determination of the "global" rho (the pt/A
312  /// of each jet is divided by this factor) and when asking for a
313  /// local rho (the result is multiplied by this factor).
314  ///
315  /// The BackgroundRescalingYPolynomial class can be used to get a
316  /// rescaling that depends just on rapidity.
317  virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) {
318  BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
319  _uptodate = false;
320  }
321 
322  //\}
323 
324  /// @name description
325  //\{
326  //----------------------------------------------------------------
327 
328  /// returns a textual description of the background estimator
329  std::string description() const;
330 
331  //\}
332 
333 
334 private:
335 
336  /// do the actual job
337  void _compute() const;
338 
339  /// check if the properties need to be recomputed
340  /// and do so if needed
341  void _recompute_if_needed() const {
342  if (!_uptodate) _compute();
343  _uptodate = true;
344  }
345 
346  /// for estimation using a selector that takes a reference jet
347  /// (i.e. a selector that can be relocated) this function allows one
348  /// to set its position.
349  ///
350  /// Note that this HAS to be called before any attempt to compute
351  /// the background properties. The call is, however, performed
352  /// automatically by the functions rho(jet) and sigma(jet).
353  void _recompute_if_needed(const PseudoJet &jet);
354 
355  /// check that the underlying structure is still alive
356  /// throw an error otherwise
357  void _check_csa_alive() const;
358 
359  /// check that the algorithm used for the clustering is adapted for
360  /// background estimation (i.e. either kt or C/A)
361  /// Issue a warning otherwise
362  void _check_jet_alg_good_for_median() const;
363 
364  // the basic parameters of this class (passed through the variou ctors)
365  Selector _rho_range; ///< range to compute the background in
366  JetDefinition _jet_def; ///< the jet def to use for teh clustering
367  AreaDefinition _area_def; ///< the area def to use for teh clustering
368  std::vector<PseudoJet> _included_jets; ///< jets to be used
369 
370  // the tunable aprameters of the class
371  bool _use_area_4vector;
372  bool _provide_fj2_sigma;
373  const FunctionOfPseudoJet<double> * _jet_density_class;
374  //SharedPtr<BackgroundRescalingBase> _rescaling_class_sharedptr;
375 
376  // the actual results of the computation
377  mutable double _rho; ///< background estimated density per unit area
378  mutable double _sigma; ///< background estimated fluctuations
379  mutable double _mean_area; ///< mean area of the jets used to estimate the background
380  mutable unsigned int _n_jets_used; ///< number of jets used to estimate the background
381  mutable double _n_empty_jets; ///< number of empty (pure-ghost) jets
382  mutable double _empty_area; ///< the empty (pure-ghost/unclustered) area!
383 
384  // internal variables
385  SharedPtr<PseudoJetStructureBase> _csi; ///< allows to check if _csa is still valid
386  PseudoJet _current_reference; ///< current reference jet
387  mutable bool _uptodate; ///< true when the background computation is up-to-date
388 
389  /// handle warning messages
390  static LimitedWarning _warnings;
391  static LimitedWarning _warnings_zero_area;
392  static LimitedWarning _warnings_preliminary;
393 };
394 
395 
396 
397 
398 //----------------------------------------------------------------------
399 /// @ingroup tools_background
400 /// \class BackgroundJetPtDensity
401 /// Class that implements pt/area_4vector.perp() for background estimation
402 /// <i>(this is a preliminary class)</i>.
404 public:
405  virtual double result(const PseudoJet & jet) const {
406  return jet.perp() / jet.area_4vector().perp();
407  }
408  virtual std::string description() const {return "BackgroundJetPtDensity";}
409 };
410 
411 
412 //----------------------------------------------------------------------
413 /// @ingroup tools_background
414 /// \class BackgroundJetScalarPtDensity
415 /// Class that implements (scalar pt sum of jet)/(scalar area of jet)
416 /// for background estimation <i>(this is a preliminary class)</i>.
417 ///
418 /// Optionally it can return a quantity based on the sum of pt^n,
419 /// e.g. for use in subtracting fragementation function moments.
421 public:
422  /// Default constructor provides background estimation with scalar pt sum
423  BackgroundJetScalarPtDensity() : _pt_power(1) {}
424 
425  /// Constructor to provide background estimation based on
426  /// \f$ sum_{i\in jet} p_{ti}^{n} \f$
427  BackgroundJetScalarPtDensity(double n) : _pt_power(n) {}
428 
429  virtual double result(const PseudoJet & jet) const;
430 
431  virtual std::string description() const {return "BackgroundScalarJetPtDensity";}
432 
433 private:
434  double _pt_power;
435 };
436 
437 //----------------------------------------------------------------------
438 /// @ingroup tools_background
439 /// \class BackgroundJetPtMDensity
440 /// Class that implements
441 /// \f$ \frac{1}{A} \sum_{i \in jet} (\sqrt{p_{ti}^2+m^2} - p_{ti}) \f$
442 /// for background estimation <i>(this is a preliminary class)</i>.
443 ///
444 ///
445 /// This is useful for correcting jet masses in cases where the event
446 /// involves massive particles.
448 public:
449  virtual double result(const PseudoJet & jet) const {
450  std::vector<PseudoJet> constituents = jet.constituents();
451  double scalar_ptm = 0;
452  for (unsigned i = 0; i < constituents.size(); i++) {
453  scalar_ptm += constituents[i].mperp() - constituents[i].perp();
454  }
455  return scalar_ptm / jet.area();
456  }
457 
458  virtual std::string description() const {return "BackgroundPtMDensity";}
459 };
460 
461 
462 
463 FASTJET_END_NAMESPACE
464 
465 #endif // __BACKGROUND_ESTIMATOR_HH__
466 
BackgroundJetScalarPtDensity(double n)
Constructor to provide background estimation based on .
Class to estimate the pt density of the background per unit area, using the median of the distributio...
virtual double result(const PseudoJet &jet) const
the action of the function this has to be overloaded in derived classes
BackgroundJetScalarPtDensity()
Default constructor provides background estimation with scalar pt sum.
double n_empty_jets() const
Returns the number of empty jets used when computing the background properties.
Class that implements pt/area_4vector.perp() for background estimation (this is a preliminary class)...
double empty_area() const
Returns the estimate of the area (within the range defined by the selector) that is not occupied by j...
void set_selector(const Selector &rho_range_selector)
(re)set the selector to be used for future calls to rho() etc.
bool use_area_4vector() const
check if the estimator uses the 4-vector area or the scalar area
virtual std::string description() const
returns a description of the function (an empty string by default)
Class that implements for background estimation (this is a preliminary class).
virtual double area() const
return the jet (scalar) area.
Definition: PseudoJet.cc:703
JetMedianBackgroundEstimator(const Selector &rho_range=SelectorIdentity())
Default constructor that optionally sets the rho range.
class that holds a generic area definition
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
unsigned int n_jets_used() const
returns the number of jets used to actually compute the background properties in the last call of rho...
void set_provide_fj2_sigma(bool provide_fj2_sigma=true)
The FastJet v2.X sigma calculation had a small spurious offset in the limit of a small number of jets...
virtual double result(const PseudoJet &jet) const
the action of the function this has to be overloaded in derived classes
base class that sets interface for extensions of ClusterSequence that provide information about the a...
const FunctionOfPseudoJet< double > * jet_density_class() const
return the pointer to the jet density class
Abstract base class that provides the basic interface for classes that estimate levels of background ...
virtual bool has_sigma()
returns true if this background estimator has support for determination of sigma
void set_use_area_4vector(bool use_it=true)
By default when calculating pt/Area for a jet, it is the transverse component of the 4-vector area th...
Class that implements (scalar pt sum of jet)/(scalar area of jet) for background estimation (this is ...
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:114
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
virtual void set_rescaling_class(const FunctionOfPseudoJet< double > *rescaling_class_in)
Set a pointer to a class that calculates the rescaling factor as a function of the jet (position)...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:566
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65
virtual std::string description() const
returns a description of the function (an empty string by default)
virtual std::string description() const
returns a description of the function (an empty string by default)
class that is intended to hold a full definition of the jet clusterer
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:718
double mean_area() const
Returns the mean area of the jets used to actually compute the background properties in the last call...