FastJet  3.0.6
JetMedianBackgroundEstimator.cc
1 //STARTHEADER
2 // $Id: JetMedianBackgroundEstimator.cc 2689 2011-11-14 14:51:06Z soyez $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
30 #include <fastjet/ClusterSequenceArea.hh>
31 #include <fastjet/ClusterSequenceStructure.hh>
32 #include <iostream>
33 
34 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
35 
36 using namespace std;
37 
38 double BackgroundJetScalarPtDensity::result(const PseudoJet & jet) const {
39  std::vector<PseudoJet> constituents = jet.constituents();
40  double scalar_pt = 0;
41  for (unsigned i = 0; i < constituents.size(); i++) {
42  scalar_pt += pow(constituents[i].perp(), _pt_power);
43  }
44  return scalar_pt / jet.area();
45 }
46 
47 
48 //----------------------------------------------------------------------
49 double BackgroundRescalingYPolynomial::result(const PseudoJet & jet) const {
50  double y = jet.rap();
51  double y2 = y*y;
52  double rescaling = _a0 + _a1*y + _a2*y2 + _a3*y2*y + _a4*y2*y2;
53  return rescaling;
54 }
55 
56 /// allow for warnings
57 LimitedWarning JetMedianBackgroundEstimator::_warnings;
58 LimitedWarning JetMedianBackgroundEstimator::_warnings_zero_area;
59 LimitedWarning JetMedianBackgroundEstimator::_warnings_preliminary;
60 
61 
62 //---------------------------------------------------------------------
63 // class JetMedianBackgroundEstimator
64 // Class to estimate the density of the background per unit area
65 //---------------------------------------------------------------------
66 
67 //----------------------------------------------------------------------
68 // ctors and dtors
69 //----------------------------------------------------------------------
70 // ctor that allows to set only the particles later on
71 JetMedianBackgroundEstimator::JetMedianBackgroundEstimator(const Selector &rho_range,
72  const JetDefinition &jet_def,
73  const AreaDefinition &area_def)
74  : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def) {
75 
76  // initialise things decently
77  reset();
78 
79  // make a few checks
80  _check_jet_alg_good_for_median();
81 }
82 
83 
84 //----------------------------------------------------------------------
85 // ctor from a cluster sequence
86 // - csa the ClusterSequenceArea to use
87 // - rho_range the Selector specifying which jets will be considered
89  : _rho_range(rho_range), _jet_def(JetDefinition()){
90 
91  // initialise things properly
92  reset();
93 
94  // tell the BGE about the cluster sequence
96 }
97 
98 
99 //----------------------------------------------------------------------
100 // setting a new event
101 //----------------------------------------------------------------------
102 // tell the background estimator that it has a new event, composed
103 // of the specified particles.
104 void JetMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
105  // make sure that we have been provided a genuine jet definition
106  if (_jet_def.jet_algorithm() == undefined_jet_algorithm)
107  throw Error("JetMedianBackgroundEstimator::set_particles can only be called if you set the jet (and area) definition explicitly through the class constructor");
108 
109  // initialise things decently (including setting uptodate to false!)
110  //reset();
111  _uptodate=false;
112 
113  // cluster the particles
114  //
115  // One may argue that it is better to cache the particles and only
116  // do the clustering later but clustering the particles now has 2
117  // practical advantages:
118  // - it allows to une only '_included_jets' in all that follows
119  // - it avoids adding another flag to ensure particles are
120  // clustered only once
121  ClusterSequenceArea *csa = new ClusterSequenceArea(particles, _jet_def, _area_def);
122  _included_jets = csa->inclusive_jets();
123 
124  // store the CS for later on
125  _csi = csa->structure_shared_ptr();
127 }
128 
129 //----------------------------------------------------------------------
130 // (re)set the cluster sequence (with area support) to be used by
131 // future calls to rho() etc.
132 //
133 // \param csa the cluster sequence area
134 //
135 // Pre-conditions:
136 // - one should be able to estimate the "empty area" (i.e. the area
137 // not occupied by jets). This is feasible if at least one of the following
138 // conditions is satisfied:
139 // ( i) the ClusterSequence has explicit ghosts
140 // (ii) the range selected has a computable area.
141 // - the jet algorithm must be suited for median computation
142 // (otherwise a warning will be issues)
143 //
144 // Note that selectors with e.g. hardest-jets exclusion do not have
145 // a well-defined area. For this reasons, it is STRONGLY advised to
146 // use an area with explicit ghosts.
148  _csi = csa.structure_shared_ptr();
149 
150  // sanity checks
151  //---------------
152  // (i) check the alg is appropriate
153  _check_jet_alg_good_for_median();
154 
155  // (ii) check that, if there are no explicit ghosts, the selector has a finite area
156  if ((!csa.has_explicit_ghosts()) && (!_rho_range.has_finite_area())){
157  throw Error("JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
158  }
159 
160  // get the initial list of jets
161  _included_jets = csa.inclusive_jets();
162 
163  _uptodate = false;
164 }
165 
166 
167 //----------------------------------------------------------------------
168 // (re)set the jets (which must have area support) to be used by future
169 // calls to rho() etc.; for the conditions that must be satisfied
170 // by the jets, see the Constructor that takes jets.
171 void JetMedianBackgroundEstimator::set_jets(const vector<PseudoJet> &jets) {
172 
173  if (! jets.size())
174  throw Error("JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: At least one jet is needed to compute the background properties");
175 
176  // sanity checks
177  //---------------
178  // (o) check that there is an underlying CS shared by all the jets
179  if (! (jets[0].has_associated_cluster_sequence()) && (jets[0].has_area()))
180  throw Error("JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
181 
182  _csi = jets[0].structure_shared_ptr();
183  ClusterSequenceStructure * csi = dynamic_cast<ClusterSequenceStructure*>(_csi());
184  const ClusterSequenceAreaBase * csab = csi->validated_csab();
185 
186  for (unsigned int i=1;i<jets.size(); i++){
187  if (! jets[i].has_associated_cluster_sequence()) // area automatic if the next test succeeds
188  throw Error("JetMedianBackgroundEstimator::set_jets(...): the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
189 
190  if (jets[i].structure_shared_ptr().get() != _csi.get())
191  throw Error("JetMedianBackgroundEstimator::set_jets(...): all the jets used to estimate the background properties must share the same ClusterSequence");
192  }
193 
194  // (i) check the alg is appropriate
195  _check_jet_alg_good_for_median();
196 
197  // (ii) check that, if there are no explicit ghosts, the selector has a finite area
198  if ((!csab->has_explicit_ghosts()) && (!_rho_range.has_finite_area())){
199  throw Error("JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
200  }
201 
202 
203  // get the initial list of jets
204  _included_jets = jets;
205 
206  // ensure recalculation of quantities that need it
207  _uptodate = false;
208 }
209 
210 
211 //----------------------------------------------------------------------
212 // retrieving fundamental information
213 //----------------------------------------------------------------
214 
215 // get rho, the median background density per unit area
217  if (_rho_range.takes_reference())
218  throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
219  _recompute_if_needed();
220  return _rho;
221 }
222 
223 // get sigma, the background fluctuations per unit area
225  if (_rho_range.takes_reference())
226  throw Error("The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
227  _recompute_if_needed();
228  return _sigma;
229 }
230 
231 // get rho, the median background density per unit area, locally at
232 // the position of a given jet.
233 //
234 // If the Selector associated with the range takes a reference jet
235 // (i.e. is relocatable), then for subsequent operations the
236 // Selector has that jet set as its reference.
238  _recompute_if_needed(jet);
239  double our_rho = _rho;
240  if (_rescaling_class != 0) {
241  our_rho *= (*_rescaling_class)(jet);
242  }
243  return our_rho;
244 }
245 
246 // get sigma, the background fluctuations per unit area,
247 // locally at the position of a given jet.
248 //
249 // If the Selector associated with the range takes a reference jet
250 // (i.e. is relocatable), then for subsequent operations the
251 // Selector has that jet set as its reference.
253  _recompute_if_needed(jet);
254  double our_sigma = _sigma;
255  if (_rescaling_class != 0) {
256  our_sigma *= (*_rescaling_class)(jet);
257  }
258  return our_sigma;
259 }
260 
261 
262 //----------------------------------------------------------------------
263 // configuring behaviour
264 //----------------------------------------------------------------------
265 // reset to default values
266 //
267 // set the variou options to their default values
269  // set the remaining default parameters
270  set_use_area_4vector(); // true by default
271  set_provide_fj2_sigma(false);
272 
273  // reset the computed values
274  _rho = _sigma = 0.0;
275  _n_jets_used = _n_empty_jets = 0;
276  _empty_area = _mean_area = 0.0;
277 
278  _included_jets.clear();
279 
280  _jet_density_class = 0; // null pointer
281  _rescaling_class = 0; // null pointer
282 
283  _uptodate = false;
284 }
285 
286 
287 // Set a pointer to a class that calculates the quantity whose
288 // median will be calculated; if the pointer is null then pt/area
289 // is used (as occurs also if this function is not called).
291  _warnings_preliminary.warn("JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.0. Their interface may differ in future releases (without guaranteeing backward compatibility).");
292  _jet_density_class = jet_density_class_in;
293  _uptodate = false;
294 }
295 
296 
297 
298 //----------------------------------------------------------------------
299 // description
300 //----------------------------------------------------------------------
302  ostringstream desc;
303  desc << "JetMedianBackgroundEstimator, using " << _jet_def.description()
304  << " with " << _area_def.description()
305  << " and selecting jets with " << _rho_range.description();
306  return desc.str();
307 }
308 
309 
310 
311 //----------------------------------------------------------------------
312 // computation of the background properties
313 //----------------------------------------------------------------------
314 // for estimation using a relocatable selector (i.e. local range)
315 // this allows to set its position. Note that this HAS to be called
316 // before any attempt to compute the background properties
317 void JetMedianBackgroundEstimator::_recompute_if_needed(const PseudoJet &jet){
318  // if the range is relocatable, handles its relocation
319  if (_rho_range.takes_reference()){
320  // check that the reference is not the same as the previous one
321  // (would avoid an unnecessary recomputation)
322  if (jet == _current_reference) return;
323 
324  // relocate the range and make sure things get recomputed the next
325  // time one tries to get some information
326  _rho_range.set_reference(jet);
327  _uptodate=false;
328  }
329 
330  _recompute_if_needed();
331 }
332 
333 // do the actual job
334 void JetMedianBackgroundEstimator::_compute() const {
335  // check if the clustersequence is still valid
336  _check_csa_alive();
337 
338  // fill the vector of pt/area (or the quantity from the jet density class)
339  // - in the range
340  vector<double> vector_for_median;
341  double total_area = 0.0;
342  _n_jets_used = 0;
343 
344  // apply the selector to the included jets
345  vector<PseudoJet> selected_jets = _rho_range(_included_jets);
346 
347  // compute the pt/area for the selected jets
348  for (unsigned i = 0; i < selected_jets.size(); i++) {
349  const PseudoJet & current_jet = selected_jets[i];
350 
351  double this_area = (_use_area_4vector) ? current_jet.area_4vector().perp() : current_jet.area();
352 
353  if (this_area>0){
354  double median_input;
355  if (_jet_density_class == 0) {
356  median_input = current_jet.perp()/this_area;
357  } else {
358  median_input = (*_jet_density_class)(current_jet);
359  }
360  if (_rescaling_class != 0) {
361  median_input /= (*_rescaling_class)(current_jet);
362  }
363  vector_for_median.push_back(median_input);
364  total_area += this_area;
365  _n_jets_used++;
366  } else {
367  _warnings_zero_area.warn("JetMedianBackgroundEstimator::_compute(...): discarded jet with zero area. Zero-area jets may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A).");
368  }
369 
370  }
371 
372  // there is nothing inside our region, so answer will always be zero
373  if (vector_for_median.size() == 0) {
374  _rho = 0.0;
375  _sigma = 0.0;
376  _mean_area = 0.0;
377  return;
378  }
379 
380  // determine the number of empty jets
381  const ClusterSequenceAreaBase * csab = (dynamic_cast<ClusterSequenceStructure*>(_csi()))->validated_csab();
382  if (csab->has_explicit_ghosts()) {
383  _empty_area = 0.0;
384  _n_empty_jets = 0;
385  } else {
386  _empty_area = csab->empty_area(_rho_range);
387  _n_empty_jets = csab->n_empty_jets(_rho_range);
388  }
389 
390  double total_njets = _n_jets_used + _n_empty_jets;
391  total_area += _empty_area;
392 
393  double stand_dev;
394  _median_and_stddev(vector_for_median, _n_empty_jets, _rho, stand_dev,
395  _provide_fj2_sigma);
396 
397  // process and store the results (_rho was already stored above)
398  _mean_area = total_area / total_njets;
399  _sigma = stand_dev * sqrt(_mean_area);
400 
401  // record that the computation has been performed
402  _uptodate = true;
403 }
404 
405 
406 
407 // check that the underlying structure is still alive;
408 // throw an error otherwise
409 void JetMedianBackgroundEstimator::_check_csa_alive() const{
410  ClusterSequenceStructure* csa = dynamic_cast<ClusterSequenceStructure*>(_csi());
411  if (csa == 0) {
412  throw Error("JetMedianBackgroundEstimator: there is no cluster sequence associated with the JetMedianBackgroundEstimator");
413  }
414  if (! dynamic_cast<ClusterSequenceStructure*>(_csi())->has_associated_cluster_sequence())
415  throw Error("JetMedianBackgroundEstimator: modifications are no longer possible as the underlying ClusterSequence has gone out of scope");
416 }
417 
418 
419 // check that the algorithm used for the clustering is suitable for
420 // background estimation (i.e. either kt or C/A).
421 // Issue a warning otherwise
422 void JetMedianBackgroundEstimator::_check_jet_alg_good_for_median() const{
423  const JetDefinition * jet_def = &_jet_def;
424 
425  // if no explicit jet def has been provided, fall back on the
426  // cluster sequence
427  if (_jet_def.jet_algorithm() == undefined_jet_algorithm){
428  const ClusterSequence * cs = dynamic_cast<ClusterSequenceStructure*>(_csi())->validated_cs();
429  jet_def = &(cs->jet_def());
430  }
431 
432  if (jet_def->jet_algorithm() != kt_algorithm
433  && jet_def->jet_algorithm() != cambridge_algorithm
435  _warnings.warn("JetMedianBackgroundEstimator: jet_def being used may not be suitable for estimating diffuse backgrounds (good alternatives are kt, cam)");
436  }
437 }
438 
439 
440 
441 
442 FASTJET_END_NAMESPACE
443 
444 
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:121
a version of cambridge with a special distance measure for particles whose pt is < extra_param() ...
std::string description() const
returns a textual description of the background estimator
deals with clustering
General class for user to obtain ClusterSequence with additional area information.
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined ...
void warn(const std::string &warning)
outputs a warning to standard error (or the user&#39;s default warning stream if set) ...
std::vector< PseudoJet > inclusive_jets(const double &ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
JetAlgorithm jet_algorithm() const
return information about the definition...
bool takes_reference() const
returns true if this can be applied jet by jet
Definition: Selector.hh:268
virtual double n_empty_jets(const Selector &selector) const
return something similar to the number of pure ghost jets in the given selector&#39;s range in an active ...
virtual double area() const
return the jet (scalar) area.
Definition: PseudoJet.cc:703
the longitudinally invariant kt algorithm
JetMedianBackgroundEstimator(const Selector &rho_range, const JetDefinition &jet_def, const AreaDefinition &area_def)
Constructor that sets the rho range as well as the jet definition and area definition to be used to c...
class that holds a generic area definition
double sigma() const
get sigma, the background fluctuations per unit area
virtual bool has_explicit_ghosts() const
returns true if ghosts are explicitly included within jets for this ClusterSequence; ...
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
Contains any information related to the clustering that should be directly accessible to PseudoJet...
virtual double empty_area(const Selector &selector) const
return the total area, corresponding to the given Selector, that is free of jets, in general based on...
const JetDefinition & jet_def() const
return a reference to the jet definition
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
void set_jet_density_class(const FunctionOfPseudoJet< double > *jet_density_class)
Set a pointer to a class that calculates the quantity whose median will be calculated; if the pointer...
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...
void _median_and_stddev(const std::vector< double > &quantity_vector, double n_empty_jets, double &median, double &stand_dev_if_gaussian, bool do_fj2_calculation=false) const
given a quantity in a vector (e.g.
std::string description() const
return a textual description of the current jet definition
base class that sets interface for extensions of ClusterSequence that provide information about the a...
std::string description() const
return a description of the current area definition
virtual void set_particles(const std::vector< PseudoJet > &particles)
tell the background estimator that it has a new event, composed of the specified particles.
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
void reset()
Resets the class to its default state, including the choice to use 4-vector areas.
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...
void set_jets(const std::vector< PseudoJet > &jets)
(re)set the jets (which must have area support) to be used by future calls to rho() etc...
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
Definition: Selector.hh:231
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:218
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
Definition: Selector.hh:273
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
double rho() const
get rho, the median background density per unit area
virtual const ClusterSequenceAreaBase * validated_csab() const
if the jet has valid area information then return a pointer to the associated ClusterSequenceAreaBase...
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
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
void set_cluster_sequence(const ClusterSequenceAreaBase &csa)
(re)set the cluster sequence (with area support) to be used by future calls to rho() etc...
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