FastJet 3.4.1
BackgroundEstimatorBase.hh
1#ifndef __FASTJET_BACKGROUND_ESTIMATOR_BASE_HH__
2#define __FASTJET_BACKGROUND_ESTIMATOR_BASE_HH__
3
4//FJSTARTHEADER
5// $Id$
6//
7// Copyright (c) 2005-2023, 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/FunctionOfPseudoJet.hh"
36#include "fastjet/Selector.hh"
37#include "fastjet/Error.hh"
38#include <iostream>
39
40FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41
42
43/// @ingroup tools_background
44/// @name helpers to handle the result of the background estimation
45//\{
46///
47/// /// a class that holds the result of the calculation
48///
49/// By default it provides access to the main background properties:
50/// rho, rho_m, sigma and sigma_m. If background estimators derived
51/// from the base class want to store more information, this can be
52/// done using the "Extra" information.
54public:
55 /// ctor wo initialisation
57 : _rho(0.0), _sigma(0.0), _rho_m(0.0), _sigma_m(0.0),
58 _has_sigma(false), _has_rho_m(false),
59 _mean_area(0.0){}
60
61
62 /// @name for accessing information about the background
63 ///@{
64
65 /// background density per unit area
66 double rho() const {return _rho;}
67
68 /// background fluctuations per unit square-root area
69 /// must be multipled by sqrt(area) to get fluctuations for a region
70 /// of a given area.
71 double sigma() const {return _sigma;}
72
73 /// true if this background estimate has a determination of sigma
74 bool has_sigma() {return true;}
75
76 /// purely longitudinal (particle-mass-induced)
77 /// component of the background density per unit area
78 double rho_m() const {return _rho_m;}
79
80 /// fluctuations in the purely longitudinal (particle-mass-induced)
81 /// component of the background density per unit square-root area
82 double sigma_m() const {return _sigma_m;}
83
84 /// true if this background estimate has a determination of rho_m.
85 /// Support for sigma_m is automatic if one has sigma and rho_m support.
86 bool has_rho_m() const {return _has_rho_m;}
87
88 /// mean area of the patches used to compute the background properties
89 double mean_area() const {return _mean_area;}
90
91 /// base class for extra information
92 class Extras {
93 public:
94 // dummy ctor
95 Extras(){};
96
97 // dummy virtual dtor
98 // makes it polymorphic to allow for dynamic_cast
99 virtual ~Extras(){};
100 };
101
102 /// returns true if the background estimate has extra info
103 bool has_extras() const{
104 return _extras.get();
105 }
106
107 /// returns true if the background estimate has extra info
108 /// compatible with the provided template type
109 template<typename T>
110 bool has_extras() const{
111 return _extras.get() && dynamic_cast<const T *>(_extras.get());
112 }
113
114 /// returns a reference to the extra information associated with a
115 /// given BackgroundEstimator. It assumes that the extra
116 /// information is reachable with class name
117 /// BackgroundEstimator::Extras
118 template<typename BackgroundEstimator>
119 const typename BackgroundEstimator::Extras & extras() const{
120 return dynamic_cast<const typename BackgroundEstimator::Extras &>(* _extras.get());
121 }
122
123 ///@}
124
125
126 /// @name for setting information about the background (internal FJ use)
127 ///@{
128
129 /// reset to default
130 void reset(){
131 _rho = _sigma = _rho_m = _sigma_m = _mean_area = 0.0;
132 _has_sigma = _has_rho_m = false;
133 _extras.reset();
134 }
135 void set_rho(double rho_in) {_rho = rho_in;}
136 void set_sigma(double sigma_in) {_sigma = sigma_in;}
137 void set_has_sigma(bool has_sigma_in) {_has_sigma = has_sigma_in;}
138 void set_rho_m(double rho_m_in) {_rho_m = rho_m_in;}
139 void set_sigma_m(double sigma_m_in) {_sigma_m = sigma_m_in;}
140 void set_has_rho_m(bool has_rho_m_in) {_has_rho_m = has_rho_m_in;}
141 void set_mean_area(double mean_area_in) {_mean_area = mean_area_in;}
142
143 /// apply a rescaling factor (to rho, rho_m, sigma, sigma_m)
144 void apply_rescaling_factor(double rescaling_factor){
145 _rho *= rescaling_factor;
146 _sigma *= rescaling_factor;
147 _rho_m *= rescaling_factor;
148 _sigma_m *= rescaling_factor;
149 }
150
151 /// sets the extra info based on the provided pointer
152 ///
153 /// When calling this method, the BackgroundEstimate class takes
154 /// ownership of the pointer (and is responsible for deleting it)
155 void set_extras(Extras *extras_in) {
156 _extras.reset(extras_in);
157 }
158 ///@}
159
160
161protected:
162 double _rho; ///< background estimated density per unit area
163 double _sigma; ///< background estimated fluctuations
164 double _rho_m; ///< "mass" background estimated density per unit area
165 double _sigma_m; ///< "mass" background estimated fluctuations
166 bool _has_sigma; ///< true if this estimate has a determination of sigma
167 bool _has_rho_m; ///< true if this estimate has a determination of rho_m
168 double _mean_area; ///< mean area of the patches used to compute the bkg properties
169
170
171 SharedPtr<Extras> _extras;
172
173};
174
175
176/// @ingroup tools_background
177/// \class BackgroundEstimatorBase
178///
179/// Abstract base class that provides the basic interface for classes
180/// that estimate levels of background radiation in hadron and
181/// heavy-ion collider events.
182///
184public:
185 /// @name constructors and destructors
186 //\{
187 //----------------------------------------------------------------
188 BackgroundEstimatorBase() : _rescaling_class(0) {
189 _set_cache_unavailable();
190 }
191
192#ifdef FASTJET_HAVE_THREAD_SAFETY
193 /// because of the internal atomic variale, we need to explicitly
194 /// implement a copy ctor
196#endif
197
198 /// a default virtual destructor that does nothing
200 //\}
201
202 /// @name setting a new event
203 //\{
204 //----------------------------------------------------------------
205
206 /// tell the background estimator that it has a new event, composed
207 /// of the specified particles.
208 virtual void set_particles(const std::vector<PseudoJet> & particles) = 0;
209
210 /// an alternative call that takes a random number generator seed
211 /// (typically a vector of length 2) to ensure reproducibility of
212 /// background estimators that rely on random numbers (specifically
213 /// JetMedianBackgroundEstimator with ghosted areas)
214 virtual void set_particles_with_seed(const std::vector<PseudoJet> & particles, const std::vector<int> & /*seed*/) {
215 set_particles(particles);
216 }
217
218 //\}
219
220 /// return a pointer to a copy of this BGE; the user is responsible
221 /// for eventually deleting the resulting object.
222 virtual BackgroundEstimatorBase * copy() const = 0;
223
224 /// @name retrieving fundamental information
225 //\{
226 //----------------------------------------------------------------
227 /// get the full set of background properties
228 virtual BackgroundEstimate estimate() const = 0;
229
230 /// get the full set of background properties for a given reference jet
231 virtual BackgroundEstimate estimate(const PseudoJet &jet) const = 0;
232
233 /// get rho, the background density per unit area
234 virtual double rho() const = 0;
235
236 /// get sigma, the background fluctuations per unit square-root area;
237 /// must be multipled by sqrt(area) to get fluctuations for a region
238 /// of a given area.
239 virtual double sigma() const {
240 throw Error("sigma() not supported for this Background Estimator");
241 }
242
243 /// get rho, the background density per unit area, locally at the
244 /// position of a given jet. Note that this is not const, because a
245 /// user may then wish to query other aspects of the background that
246 /// could depend on the position of the jet last used for a rho(jet)
247 /// determination.
248 virtual double rho(const PseudoJet & jet) = 0;
249
250 /// get sigma, the background fluctuations per unit area, locally at
251 /// the position of a given jet. As for rho(jet), it is non-const.
252 virtual double sigma(const PseudoJet & /*jet*/) {
253 throw Error("sigma(jet) not supported for this Background Estimator");
254 }
255
256 /// returns true if this background estimator has support for
257 /// determination of sigma
258 virtual bool has_sigma() const {return false;}
259
260 //----------------------------------------------------------------
261 // now do the same thing for rho_m and sigma_m
262
263 /// returns rho_m, the purely longitudinal, particle-mass-induced
264 /// component of the background density per unit area
265 virtual double rho_m() const{
266 throw Error("rho_m() not supported for this Background Estimator");
267 }
268
269 /// returns sigma_m, a measure of the fluctuations in the purely
270 /// longitudinal, particle-mass-induced component of the background
271 /// density per unit square-root area; must be multipled by sqrt(area) to get
272 /// fluctuations for a region of a given area.
273 virtual double sigma_m() const {
274 throw Error("sigma_m() not supported for this Background Estimator");
275 }
276
277 /// Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
278 virtual double rho_m(const PseudoJet & /*jet*/){
279 throw Error("rho_m(jet) not supported for this Background Estimator");
280 }
281
282 /// Returns sigma_m locally at the jet position. As for rho(jet), it is non-const.
283 virtual double sigma_m(const PseudoJet & /*jet*/) {
284 throw Error("sigma_m(jet) not supported for this Background Estimator");
285 }
286
287 /// Returns true if this background estimator has support for
288 /// determination of rho_m.
289 ///
290 /// Note that support for sigma_m is automatic is one has sigma and
291 /// rho_m support.
292 virtual bool has_rho_m() const {return false;}
293 //\}
294
295
296 /// @name configuring the behaviour
297 //\{
298 //----------------------------------------------------------------
299
300 /// Set a pointer to a class that calculates the rescaling factor as
301 /// a function of the jet (position). Note that the rescaling factor
302 /// is used both in the determination of the "global" rho (the pt/A
303 /// of each jet is divided by this factor) and when asking for a
304 /// local rho (the result is multiplied by this factor).
305 ///
306 /// The BackgroundRescalingYPolynomial class can be used to get a
307 /// rescaling that depends just on rapidity.
308 ///
309 /// There is currently no support for different rescaling classes
310 /// for rho and rho_m determinations.
311 virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) { _rescaling_class = rescaling_class_in; }
312
313 /// return the pointer to the jet density class
315 return _rescaling_class;
316 }
317
318 //\}
319
320 /// @name description
321 //\{
322 //----------------------------------------------------------------
323
324 /// returns a textual description of the background estimator
325 virtual std::string description() const = 0;
326
327 //\}
328
329
330protected:
331 /// @name helpers for derived classes
332 ///
333 /// Note that these helpers are related to median-based estimation
334 /// of the background, so there is no guarantee that they will
335 /// remain in this base class in the long term
336 //\{
337 //----------------------------------------------------------------
338
339 /// given a quantity in a vector (e.g. pt_over_area) and knowledge
340 /// about the number of empty jets, calculate the median and
341 /// stand_dev_if_gaussian (roughly from the 16th percentile)
342 ///
343 /// If do_fj2_calculation is set to true then this performs FastJet
344 /// 2.X estimation of the standard deviation, which has a spurious
345 /// offset in the limit of a small number of jets.
346 void _median_and_stddev(const std::vector<double> & quantity_vector,
347 double n_empty_jets,
348 double & median,
349 double & stand_dev_if_gaussian,
350 bool do_fj2_calculation = false
351 ) const;
352
353 /// computes a percentile of a given _sorted_ vector
354 /// \param sorted_quantity_vector the vector contains the data sample
355 /// \param percentile the percentile (defined between 0 and 1) to compute
356 /// \param nempty an additional number of 0's
357 /// (considered at the beginning of
358 /// the quantity vector)
359 /// \param do_fj2_calculation carry out the calculation as it
360 /// was done in fj2 (suffers from "edge effects")
361 double _percentile(const std::vector<double> & sorted_quantity_vector,
362 const double percentile,
363 const double nempty=0.0,
364 const bool do_fj2_calculation = false) const;
365
366 //\}
367
368 const FunctionOfPseudoJet<double> * _rescaling_class;
369 static LimitedWarning _warnings_empty_area;
370
371
372 // cached actual results of the computation
373 mutable bool _cache_available;
374 mutable BackgroundEstimate _cached_estimate; ///< all the info about what is computed
375 //\}
376
377 /// @name helpers for thread safety
378 ///
379 /// Note that these helpers are related to median-based estimation
380 /// of the background, so there is no guarantee that they will
381 /// remain in this base class in the long term
382 //\{
383#ifdef FASTJET_HAVE_THREAD_SAFETY
384 // true when one is currently writing to cache (i.e. when the spin lock is set)
385 mutable std::atomic<bool> _writing_to_cache;
386
387 void _set_cache_unavailable(){
388 _cache_available = false;
389 _writing_to_cache = false;
390 }
391
392 // // allows us to lock things down before caching basic (patches) info
393 // std::mutex _jets_caching_mutex;
394#else
395 void _set_cache_unavailable(){
396 _cache_available = false;
397 }
398#endif
399
400 void _lock_if_needed() const;
401 void _unlock_if_needed() const;
402
403 //\}
404
405};
406
407
408
409//----------------------------------------------------------------------
410/// @ingroup tools_background
411/// A background rescaling that is a simple polynomial in y
413public:
414 /// construct a background rescaling polynomial of the form
415 /// a0 + a1*y + a2*y^2 + a3*y^3 + a4*y^4
416 ///
417 /// The following values give a reasonable reproduction of the
418 /// Pythia8 tune 4C background shape for pp collisions at
419 /// sqrt(s)=7TeV:
420 ///
421 /// - a0 = 1.157
422 /// - a1 = 0
423 /// - a2 = -0.0266
424 /// - a3 = 0
425 /// - a4 = 0.000048
426 ///
428 double a1=0,
429 double a2=0,
430 double a3=0,
431 double a4=0) : _a0(a0), _a1(a1), _a2(a2), _a3(a3), _a4(a4) {}
432
433 /// return the rescaling factor associated with this jet
434 virtual double result(const PseudoJet & jet) const;
435private:
436 double _a0, _a1, _a2, _a3, _a4;
437};
438
439
440
441
442
443FASTJET_END_NAMESPACE
444
445#endif // __BACKGROUND_ESTIMATOR_BASE_HH__
446
447// //--------------------------------------------------
448// // deprecated
449// class JetMedianBGE{
450// BackgroundEstimateDefinition();
451//
452// ....;
453// }
454// //--------------------------------------------------
455//
456// class BackgroundEstimateDefinition{
457// //const EventBackground get_event_background(particles, <seed>) const;
458//
459//
460// //--------------------------------------------------
461// // DEPRECATED
462// void set_particles() {
463//
464// _worker = ...;
465// double rho(const PseudoJet &/*jet*/) const{ _worker.rho();}
466//
467// //--------------------------------------------------
468//
469// EventBackground(Worker?) _cached_event_background;
470// };
471//
472// class EventBackground{
473// EventBackground(particles, BackgroundEstimateDefinition);
474//
475//
476// class EventBackgroundWorker{
477//
478// ...
479// };
480//
481//
482// BackgroundEstimate estimate() const;
483// BackgroundEstimate estimate(jet) const;
484//
485// // do we want this:
486// double rho();
487// double rho(jet);
488// //?
489//
490//
491// mutable BackgroundEstimate _estimate;
492//
493//
494// SharedPtr<EventBackgroundWorker> _event_background_worker;
495// }
496//
497// class BackgroundEstimate{
498//
499// double rho();
500//
501// SharedPtr<EventBackgroundWorker> _event_background_worker;
502//
503// private:
504// _rho;
505// _sigma;
506// _rho_m;
507// _sigma_m;
508// _has_rho_m;
509// _has_sigma;
510//
511// // info specific to JMBGE: _reference_jet, mean_area, n_jets_used, n_empty_jets, _empty_area
512// // all need to go in the estimate in general
513// // info specific to GMBGE: RectangularGrid, _mean_area (can go either in the the def, the eventBG or the estimate
514//
515// };
base class for extra information
/// a class that holds the result of the calculation
void set_extras(Extras *extras_in)
sets the extra info based on the provided pointer
double _mean_area
mean area of the patches used to compute the bkg properties
bool has_rho_m() const
true if this background estimate has a determination of rho_m.
bool has_extras() const
returns true if the background estimate has extra info compatible with the provided template type
double _sigma_m
"mass" background estimated fluctuations
bool _has_rho_m
true if this estimate has a determination of rho_m
double sigma_m() const
fluctuations in the purely longitudinal (particle-mass-induced) component of the background density p...
double _rho_m
"mass" background estimated density per unit area
void apply_rescaling_factor(double rescaling_factor)
apply a rescaling factor (to rho, rho_m, sigma, sigma_m)
double rho() const
background density per unit area
const BackgroundEstimator::Extras & extras() const
returns a reference to the extra information associated with a given BackgroundEstimator.
bool has_sigma()
true if this background estimate has a determination of sigma
double sigma() const
background fluctuations per unit square-root area must be multipled by sqrt(area) to get fluctuations...
bool has_extras() const
returns true if the background estimate has extra info
double rho_m() const
purely longitudinal (particle-mass-induced) component of the background density per unit area
double _sigma
background estimated fluctuations
double _rho
background estimated density per unit area
BackgroundEstimate()
ctor wo initialisation
double mean_area() const
mean area of the patches used to compute the background properties
bool _has_sigma
true if this estimate has a determination of sigma
Abstract base class that provides the basic interface for classes that estimate levels of background ...
virtual double sigma() const
get sigma, the background fluctuations per unit square-root area; must be multipled by sqrt(area) to ...
virtual std::string description() const =0
returns a textual description of the background estimator
const FunctionOfPseudoJet< double > * rescaling_class() const
return the pointer to the jet density class
virtual double sigma_m() const
returns sigma_m, a measure of the fluctuations in the purely longitudinal, particle-mass-induced comp...
virtual BackgroundEstimate estimate() const =0
get the full set of background properties
virtual double rho() const =0
get rho, the background density per unit area
virtual double sigma_m(const PseudoJet &)
Returns sigma_m locally at the jet position. As for rho(jet), it is non-const.
virtual double sigma(const PseudoJet &)
get sigma, the background fluctuations per unit area, locally at the position of a given jet.
virtual ~BackgroundEstimatorBase()
a default virtual destructor that does nothing
virtual double rho_m() const
returns rho_m, the purely longitudinal, particle-mass-induced component of the background density per...
virtual double rho_m(const PseudoJet &)
Returns rho_m locally at the jet position. As for rho(jet), it is non-const.
virtual bool has_rho_m() const
Returns true if this background estimator has support for determination of rho_m.
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).
BackgroundEstimate _cached_estimate
all the info about what is computed
virtual double rho(const PseudoJet &jet)=0
get rho, the background density per unit area, locally at the position of a given jet.
virtual void set_particles_with_seed(const std::vector< PseudoJet > &particles, const std::vector< int > &)
an alternative call that takes a random number generator seed (typically a vector of length 2) to ens...
virtual BackgroundEstimate estimate(const PseudoJet &jet) const =0
get the full set of background properties for a given reference jet
virtual BackgroundEstimatorBase * copy() const =0
return a pointer to a copy of this BGE; the user is responsible for eventually deleting the resulting...
virtual void set_particles(const std::vector< PseudoJet > &particles)=0
tell the background estimator that it has a new event, composed of the specified particles.
virtual bool has_sigma() const
returns true if this background estimator has support for determination of sigma
A background rescaling that is a simple polynomial in y.
BackgroundRescalingYPolynomial(double a0=1, double a1=0, double a2=0, double a3=0, double a4=0)
construct a background rescaling polynomial of the form a0 + a1*y + a2*y^2 + a3*y^3 + a4*y^4
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:52
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition: PseudoJet.hh:68
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition: SharedPtr.hh:341