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
434  && jet_def->jet_algorithm() != cambridge_for_passive_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
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'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 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...
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...
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