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