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 | |
---|
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 |
---|
88 | JetMedianBackgroundEstimator::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. |
---|
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(); |
---|
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. |
---|
147 | void 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. |
---|
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 |
---|
216 | double 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 |
---|
224 | double 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. |
---|
237 | double 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. |
---|
252 | double 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 |
---|
268 | void 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). |
---|
290 | void 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 | //---------------------------------------------------------------------- |
---|
301 | string 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 |
---|
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 | |
---|