Fork me on GitHub

source: svn/trunk/external/fastjet/tools/JetMedianBackgroundEstimator.cc@ 1368

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

add fastjet/tools

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 16.8 KB
Line 
1//STARTHEADER
2// $Id: JetMedianBackgroundEstimator.cc 999 2013-03-04 11:48:06Z pavel $
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
34FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
35
36using namespace std;
37
38double 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//----------------------------------------------------------------------
49double 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
57LimitedWarning JetMedianBackgroundEstimator::_warnings;
58LimitedWarning JetMedianBackgroundEstimator::_warnings_zero_area;
59LimitedWarning 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
71JetMedianBackgroundEstimator::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
88JetMedianBackgroundEstimator::JetMedianBackgroundEstimator( const Selector &rho_range, const ClusterSequenceAreaBase &csa)
89 : _rho_range(rho_range), _jet_def(JetDefinition()){
90
91 // initialise things properly
92 reset();
93
94 // tell the BGE about the cluster sequence
95 set_cluster_sequence(csa);
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.
104void 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();
126 csa->delete_self_when_unused();
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.
147void JetMedianBackgroundEstimator::set_cluster_sequence(const ClusterSequenceAreaBase & csa) {
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.
171void 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
216double JetMedianBackgroundEstimator::rho() const {
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
224double JetMedianBackgroundEstimator::sigma() const {
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.
237double JetMedianBackgroundEstimator::rho(const PseudoJet & jet) {
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.
252double JetMedianBackgroundEstimator::sigma(const PseudoJet &jet) {
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
268void JetMedianBackgroundEstimator::reset(){
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).
290void JetMedianBackgroundEstimator::set_jet_density_class(const FunctionOfPseudoJet<double> * jet_density_class_in) {
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//----------------------------------------------------------------------
301string JetMedianBackgroundEstimator::description() const {
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
317void 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
334void 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
409void 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
422void 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
442FASTJET_END_NAMESPACE
443
444
Note: See TracBrowser for help on using the repository browser.