Fork me on GitHub

source: git/external/fastjet/ClusterSequenceAreaBase.cc@ 78d1846

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 78d1846 was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4
  • Property mode set to 100644
File size: 18.3 KB
Line 
1
2//FJSTARTHEADER
3// $Id: ClusterSequenceAreaBase.cc 4079 2016-03-09 12:20:03Z soyez $
4//
5// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
6//
7//----------------------------------------------------------------------
8// This file is part of FastJet.
9//
10// FastJet is free software; you can redistribute it and/or modify
11// it under the terms of the GNU General Public License as published by
12// the Free Software Foundation; either version 2 of the License, or
13// (at your option) any later version.
14//
15// The algorithms that underlie FastJet have required considerable
16// development. They are described in the original FastJet paper,
17// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
18// FastJet as part of work towards a scientific publication, please
19// quote the version you use and include a citation to the manual and
20// optionally also to hep-ph/0512210.
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//FJENDHEADER
31
32
33
34
35#include "fastjet/ClusterSequenceAreaBase.hh"
36#include <algorithm>
37
38FASTJET_BEGIN_NAMESPACE
39
40using namespace std;
41
42
43/// allow for warnings
44LimitedWarning ClusterSequenceAreaBase::_warnings;
45LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
46LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;
47
48//----------------------------------------------------------------------
49/// return the total area, within the selector's range, that is free
50/// of jets.
51///
52/// Calculate this as (range area) - \sum_{i in range} A_i
53///
54/// for ClusterSequences with explicit ghosts, assume that there will
55/// never be any empty area, i.e. it is always filled in by pure
56/// ghosts jets. This holds for seq.rec. algorithms
57double ClusterSequenceAreaBase::empty_area(const Selector & selector) const {
58
59 if (has_explicit_ghosts()) {return 0.0;}
60 else { return empty_area_from_jets(inclusive_jets(0.0), selector);}
61
62}
63
64//----------------------------------------------------------------------
65/// return the total area, within range, that is free of jets.
66///
67/// Calculate this as (range area) - \sum_{i in range} A_i
68///
69double ClusterSequenceAreaBase::empty_area_from_jets(
70 const std::vector<PseudoJet> & all_jets,
71 const Selector & selector) const {
72 _check_selector_good_for_median(selector);
73
74 double empty = selector.area();
75 for (unsigned i = 0; i < all_jets.size(); i++) {
76 if (selector.pass(all_jets[i])) empty -= area(all_jets[i]);
77 }
78 return empty;
79}
80
81// this is deprecated but used by other deprecated methods. So we hide
82// the implementation in a protected method so that (i) it can still
83// be used internally (without generating a compile-time warning when
84// building FastJet) and the interface can be marked as deprecated.
85// This can disappear once all the public interfaces have disappeared.
86double ClusterSequenceAreaBase::median_pt_per_unit_area(const Selector & selector) const {
87 return _median_pt_per_unit_area(selector);
88}
89
90// the hidden implementation
91double ClusterSequenceAreaBase::_median_pt_per_unit_area(const Selector & selector) const {
92 return _median_pt_per_unit_something(selector,false);
93}
94
95
96
97// this is deprecated but used by other deprecated methods. So we hide
98// the implementation in a protected method so that (i) it can still
99// be used internally (without generating a compile-time warning when
100// building FastJet) and the interface can be marked as deprecated.
101// This can disappear once all the public interfaces have disappeared.
102double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const Selector & selector) const {
103 return _median_pt_per_unit_area_4vector(selector);
104}
105
106// the deprecated interface
107double ClusterSequenceAreaBase::_median_pt_per_unit_area_4vector(const Selector & selector) const {
108 return _median_pt_per_unit_something(selector,true);
109}
110
111
112//----------------------------------------------------------------------
113// this is deprecated but used by other deprecated methods. So we hide
114// the implementation in a protected method so that (i) it can still
115// be used internally (without generating a compile-time warning when
116// building FastJet) and the interface can be marked as deprecated.
117// This can disappear once all the public interfaces have disappeared.
118double ClusterSequenceAreaBase::median_pt_per_unit_something(
119 const Selector & selector, bool use_area_4vector) const {
120 return _median_pt_per_unit_something(selector, use_area_4vector);
121}
122
123// the median of (pt/area) for jets contained within range, counting
124// the empty area as if it were made up of a collection of empty
125// jets each of area (0.55 * pi R^2).
126double ClusterSequenceAreaBase::_median_pt_per_unit_something(
127 const Selector & selector, bool use_area_4vector) const {
128 double median, sigma, mean_area;
129 _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
130 return median;
131}
132
133
134//----------------------------------------------------------------------
135/// fits a form pt_per_unit_area(y) = a + b*y^2 for jets in range.
136/// exclude_above allows one to exclude large values of pt/area from
137/// fit. use_area_4vector = true uses the 4vector areas.
138void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
139 double & a, double & b, const Selector & selector,
140 double exclude_above, bool use_area_4vector) const {
141 return _parabolic_pt_per_unit_area(a, b, selector, exclude_above, use_area_4vector);
142}
143
144void ClusterSequenceAreaBase::_parabolic_pt_per_unit_area(
145 double & a, double & b, const Selector & selector,
146 double exclude_above, bool use_area_4vector) const {
147 // sanity check on the selector: we require a finite area and that
148 // it applies jet by jet (see BackgroundEstimator for more advanced
149 // usage)
150 _check_selector_good_for_median(selector);
151
152 int n=0;
153 int n_excluded = 0;
154 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
155
156 vector<PseudoJet> incl_jets = inclusive_jets();
157
158 for (unsigned i = 0; i < incl_jets.size(); i++) {
159 if (selector.pass(incl_jets[i])) {
160 double this_area;
161 if ( use_area_4vector ) {
162 this_area = area_4vector(incl_jets[i]).perp();
163 } else {
164 this_area = area(incl_jets[i]);
165 }
166 double f = incl_jets[i].perp()/this_area;
167 if (exclude_above <= 0.0 || f < exclude_above) {
168 double x = incl_jets[i].rap(); double x2 = x*x;
169 mean_f += f;
170 mean_x2 += x2;
171 mean_x4 += x2*x2;
172 mean_fx2 += f*x2;
173 n++;
174 } else {
175 n_excluded++;
176 }
177 }
178 }
179
180 if (n <= 1) {
181 // meaningful results require at least two jets inside the
182 // area -- mind you if there are empty jets we should be in
183 // any case doing something special...
184 a = 0.0;
185 b = 0.0;
186 } else {
187 mean_f /= n;
188 mean_x2 /= n;
189 mean_x4 /= n;
190 mean_fx2 /= n;
191
192 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
193 a = mean_f - b*mean_x2;
194 }
195 //cerr << "n_excluded = "<< n_excluded << endl;
196}
197
198
199//----------------------------------------------------------------------
200void ClusterSequenceAreaBase::get_median_rho_and_sigma(
201 const Selector & selector, bool use_area_4vector,
202 double & median, double & sigma, double & mean_area) const {
203 _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
204}
205
206void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
207 const Selector & selector, bool use_area_4vector,
208 double & median, double & sigma, double & mean_area) const {
209
210 vector<PseudoJet> incl_jets = inclusive_jets();
211 _get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
212 median, sigma, mean_area, true);
213}
214
215void ClusterSequenceAreaBase::get_median_rho_and_sigma(
216 const vector<PseudoJet> & all_jets,
217 const Selector & selector, bool use_area_4vector,
218 double & median, double & sigma, double & mean_area,
219 bool all_are_incl) const {
220 _get_median_rho_and_sigma(all_jets, selector, use_area_4vector,
221 median, sigma, mean_area, all_are_incl);
222}
223
224void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
225 const vector<PseudoJet> & all_jets,
226 const Selector & selector, bool use_area_4vector,
227 double & median, double & sigma, double & mean_area,
228 bool all_are_incl) const {
229
230 _check_jet_alg_good_for_median();
231
232 // sanity check on the selector: we require a finite area and that
233 // it applies jet by jet (see BackgroundEstimator for more advanced
234 // usage)
235 _check_selector_good_for_median(selector);
236
237 vector<double> pt_over_areas;
238 double total_area = 0.0;
239 double total_njets = 0;
240
241 for (unsigned i = 0; i < all_jets.size(); i++) {
242 if (selector.pass(all_jets[i])) {
243 double this_area;
244 if (use_area_4vector) {
245 this_area = area_4vector(all_jets[i]).perp();
246 } else {
247 this_area = area(all_jets[i]);
248 }
249
250 if (this_area>0) {
251 pt_over_areas.push_back(all_jets[i].perp()/this_area);
252 } else {
253 _warnings_zero_area.warn("ClusterSequenceAreaBase::get_median_rho_and_sigma(...): 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).");
254 }
255
256 total_area += this_area;
257 total_njets += 1.0;
258 }
259 }
260
261 // there is nothing inside our region, so answer will always be zero
262 if (pt_over_areas.size() == 0) {
263 median = 0.0;
264 sigma = 0.0;
265 mean_area = 0.0;
266 return;
267 }
268
269 // get median (pt/area) [this is the "old" median definition. It considers
270 // only the "real" jets in calculating the median, i.e. excluding the
271 // only-ghost ones; it will be supplemented with more info below]
272 sort(pt_over_areas.begin(), pt_over_areas.end());
273
274 // now get the median & error, accounting for empty jets
275 // define the fractions of distribution at median, median-1sigma
276 double posn[2] = {0.5, (1.0-0.6827)/2.0};
277 double res[2];
278
279 double n_empty, empty_a;
280 if (has_explicit_ghosts()) {
281 // NB: the following lines of code are potentially incorrect in cases
282 // where there are unclustered particles (empty_area would do a better job,
283 // at least for active areas). This is not an issue with kt or C/A, or other
284 // algorithms that cluster all particles (and the median estimation should in
285 // any case only be done with kt or C/A!)
286 empty_a = 0.0;
287 n_empty = 0;
288 } else if (all_are_incl) {
289 // the default case
290 empty_a = empty_area(selector);
291 n_empty = n_empty_jets(selector);
292 } else {
293 // this one is intended to be used when e.g. one runs C/A, then looks at its
294 // exclusive jets in order to get an effective smaller R value, and passes those
295 // to this routine.
296 empty_a = empty_area_from_jets(all_jets, selector);
297 mean_area = total_area / total_njets; // temporary value
298 n_empty = empty_a / mean_area;
299 }
300 //cout << "*** tot_area = " << total_area << ", empty_a = " << empty_a << endl;
301 //cout << "*** n_empty = " << n_empty << ", ntotal = " << total_njets << endl;
302 total_njets += n_empty;
303 total_area += empty_a;
304
305 // we need an int (rather than an unsigned int) with the size of the
306 // pt_over_areas array, because we'll often be doing subtraction of
307 // -1, negating it, etc. All of these operations go crazy with unsigned ints.
308 int pt_over_areas_size = pt_over_areas.size();
309 if (n_empty < -pt_over_areas_size/4.0)
310 _warnings_empty_area.warn("ClusterSequenceAreaBase::get_median_rho_and_sigma(...): the estimated empty area is suspiciously large and negative and may lead to an over-estimation of rho. This may be due to (i) a rare statistical fluctuation or (ii) too small a range used to estimate the background properties.");
311
312 for (int i = 0; i < 2; i++) {
313 double nj_median_pos =
314 (pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
315 double nj_median_ratio;
316 if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
317 int int_nj_median = int(nj_median_pos);
318
319 // avoid potential overflow issues
320 if (int_nj_median+1 > pt_over_areas_size-1){
321 int_nj_median = pt_over_areas_size-2;
322 nj_median_pos = pt_over_areas_size-1;
323 }
324
325 nj_median_ratio =
326 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
327 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
328 } else {
329 nj_median_ratio = 0.0;
330 }
331 res[i] = nj_median_ratio;
332 }
333 median = res[0];
334 double error = res[0] - res[1];
335 mean_area = total_area / total_njets;
336 sigma = error * sqrt(mean_area);
337}
338
339
340// return a vector of all subtracted jets, using area_4vector, given rho.
341// Only inclusive_jets above ptmin are subtracted and returned.
342// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
343// i.e. not necessarily ordered in pt once subtracted
344vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
345 const double ptmin)
346 const {
347 return _subtracted_jets(rho,ptmin);
348}
349
350vector<PseudoJet> ClusterSequenceAreaBase::_subtracted_jets(const double rho,
351 const double ptmin)
352 const {
353 vector<PseudoJet> sub_jets;
354 vector<PseudoJet> jets_local = sorted_by_pt(inclusive_jets(ptmin));
355 for (unsigned i=0; i<jets_local.size(); i++) {
356 PseudoJet sub_jet = _subtracted_jet(jets_local[i],rho);
357 sub_jets.push_back(sub_jet);
358 }
359 return sub_jets;
360}
361
362// return a vector of subtracted jets, using area_4vector.
363// Only inclusive_jets above ptmin are subtracted and returned.
364// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
365// i.e. not necessarily ordered in pt once subtracted
366vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
367 const Selector & selector,
368 const double ptmin)
369 const {
370 double rho = _median_pt_per_unit_area_4vector(selector);
371 return _subtracted_jets(rho,ptmin);
372}
373
374
375/// return a subtracted jet, using area_4vector, given rho
376PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
377 const double rho) const {
378 return _subtracted_jet(jet, rho);
379}
380
381PseudoJet ClusterSequenceAreaBase::_subtracted_jet(const PseudoJet & jet,
382 const double rho) const {
383 PseudoJet area4vect = area_4vector(jet);
384 PseudoJet sub_jet;
385 // sanity check
386 if (rho*area4vect.perp() < jet.perp() ) {
387 sub_jet = jet - rho*area4vect;
388 } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
389
390 // make sure the subtracted jet has the same index (cluster, user, csw)
391 // (i.e. "looks like") the original jet
392 sub_jet.set_cluster_hist_index(jet.cluster_hist_index());
393 sub_jet.set_user_index(jet.user_index());
394 // do not use CS::_set_structure_shared_ptr here, which should
395 // only be called to maintain the tally during construction
396 sub_jet.set_structure_shared_ptr(jet.structure_shared_ptr());
397 return sub_jet;
398}
399
400
401/// return a subtracted jet, using area_4vector; note that this is
402/// potentially inefficient if repeatedly used for many different
403/// jets, because rho will be recalculated each time around.
404PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
405 const Selector & selector) const {
406 return _subtracted_jet(jet, selector);
407}
408
409PseudoJet ClusterSequenceAreaBase::_subtracted_jet(const PseudoJet & jet,
410 const Selector & selector) const {
411 double rho = _median_pt_per_unit_area_4vector(selector);
412 PseudoJet sub_jet = _subtracted_jet(jet, rho);
413 return sub_jet;
414}
415
416
417/// return the subtracted pt, given rho
418double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
419 const double rho,
420 bool use_area_4vector) const {
421 return _subtracted_pt(jet, rho, use_area_4vector);
422}
423
424double ClusterSequenceAreaBase::_subtracted_pt(const PseudoJet & jet,
425 const double rho,
426 bool use_area_4vector) const {
427 if ( use_area_4vector ) {
428 PseudoJet sub_jet = _subtracted_jet(jet,rho);
429 return sub_jet.perp();
430 } else {
431 return jet.perp() - rho*area(jet);
432 }
433}
434
435
436/// return the subtracted pt; note that this is
437/// potentially inefficient if repeatedly used for many different
438/// jets, because rho will be recalculated each time around.
439double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
440 const Selector & selector,
441 bool use_area_4vector) const {
442 if ( use_area_4vector ) {
443 PseudoJet sub_jet = _subtracted_jet(jet,selector);
444 return sub_jet.perp();
445 } else {
446 double rho = _median_pt_per_unit_area(selector);
447 return _subtracted_pt(jet,rho,false);
448 }
449}
450
451// check the selector is suited for the computations i.e. applies jet
452// by jet and has a finite area
453void ClusterSequenceAreaBase::_check_selector_good_for_median(const Selector &selector) const{
454 // make sure the selector has a finite area
455 if ((! has_explicit_ghosts()) && (! selector.has_finite_area())){
456 throw Error("ClusterSequenceAreaBase: empty area can only be computed from selectors with a finite area");
457 }
458
459 // make sure the selector applies jet by jet
460 if (! selector.applies_jet_by_jet()){
461 throw Error("ClusterSequenceAreaBase: empty area can only be computed from selectors that apply jet by jet");
462 }
463}
464
465
466/// check the jet algorithm is suitable (and if not issue a warning)
467void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
468 if (jet_def().jet_algorithm() != kt_algorithm
469 && jet_def().jet_algorithm() != cambridge_algorithm
470 && jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
471 _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
472 }
473}
474
475
476
477FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.