Fork me on GitHub

source: svn/trunk/external/fastjet/tools/JetMedianBackgroundEstimator.hh@ 1204

Last change on this file since 1204 was 999, checked in by Pavel Demin, 12 years ago

add fastjet/tools

File size: 17.6 KB
Line 
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
39FASTJET_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///
78class JetMedianBackgroundEstimator : public BackgroundEstimatorBase {
79public:
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
127 ~JetMedianBackgroundEstimator(){}
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
305 const FunctionOfPseudoJet<double> * jet_density_class() const{
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
334private:
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>.
403class BackgroundJetPtDensity : public FunctionOfPseudoJet<double> {
404public:
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.
420class BackgroundJetScalarPtDensity : public FunctionOfPseudoJet<double> {
421public:
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
433private:
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.
447class BackgroundJetPtMDensity : public FunctionOfPseudoJet<double> {
448public:
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
463FASTJET_END_NAMESPACE
464
465#endif // __BACKGROUND_ESTIMATOR_HH__
466
Note: See TracBrowser for help on using the repository browser.