1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
//STARTHEADER
// $Id: BackgroundEstimatorBase.cc 2616 2011-09-30 18:03:40Z salam $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
#include "fastjet/tools/BackgroundEstimatorBase.hh"
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
LimitedWarning BackgroundEstimatorBase::_warnings_empty_area;
//----------------------------------------------------------------------
// given a quantity in a vector (e.g. pt_over_area) and knowledge
// about the number of empty jets, calculate the median and
// stand_dev_if_gaussian (roughly from the 16th percentile)
//
// If do_fj2_calculation is set to true then this performs FastJet
// 2.X estimation of the standard deviation, which has a spurious
// offset in the limit of a small number of jets.
void BackgroundEstimatorBase::_median_and_stddev(const vector<double> & quantity_vector,
double n_empty_jets,
double & median,
double & stand_dev_if_gaussian,
bool do_fj2_calculation) const {
// this check is redundant (the code below behaves sensibly even
// with a zero size), but serves as a reminder of what happens if
// the quantity vector is zero-sized
if (quantity_vector.size() == 0) {
median = 0;
stand_dev_if_gaussian = 0;
return;
}
vector<double> sorted_quantity_vector = quantity_vector;
sort(sorted_quantity_vector.begin(), sorted_quantity_vector.end());
// empty area can sometimes be negative; with small ranges this can
// become pathological, so warn the user
int n_jets_used = sorted_quantity_vector.size();
if (n_empty_jets < -n_jets_used/4.0)
_warnings_empty_area.warn("BackgroundEstimatorBase::_median_and_stddev(...): 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.");
// now get the median & error, accounting for empty jets;
// define the fractions of distribution at median, median-1sigma
double posn[2] = {0.5, (1.0-0.6827)/2.0};
double res[2];
for (int i = 0; i < 2; i++) {
res[i] = _percentile(sorted_quantity_vector, posn[i], n_empty_jets,
do_fj2_calculation);
}
median = res[0];
stand_dev_if_gaussian = res[0] - res[1];
}
//----------------------------------------------------------------------
// computes a percentile of a given _sorted_ vector of quantities
// - sorted_quantities the (sorted) vector contains the data sample
// - percentile the percentile (defined between 0 and 1) to compute
// - nempty an additional number of 0's
// (considered at the beginning of
// the quantity vector)
// - do_fj2_calculation carry out the calculation as it
// was done in fj2 (suffers from "edge effects")
double BackgroundEstimatorBase::_percentile(const vector<double> & sorted_quantities,
const double percentile,
const double nempty,
const bool do_fj2_calculation
) const {
assert(percentile >= 0.0 && percentile <= 1.0);
int quantities_size = sorted_quantities.size();
if (quantities_size == 0) return 0;
double total_njets = quantities_size + nempty;
double percentile_pos;
if (do_fj2_calculation) {
percentile_pos = (total_njets-1)*percentile - nempty;
} else {
percentile_pos = (total_njets)*percentile - nempty - 0.5;
}
double result;
if (percentile_pos >= 0 && quantities_size > 1) {
int int_percentile_pos = int(percentile_pos);
// avoid potential overflow issues
if (int_percentile_pos+1 > quantities_size-1){
int_percentile_pos = quantities_size-2;
percentile_pos = quantities_size-1;
}
result =
sorted_quantities[int_percentile_pos] * (int_percentile_pos+1-percentile_pos)
+ sorted_quantities[int_percentile_pos+1] * (percentile_pos - int_percentile_pos);
} else if (percentile_pos > -0.5 && quantities_size >= 1
&& !do_fj2_calculation) {
// in the LHS of this "bin", just keep a constant value (we could have
// interpolated to zero, but this might misbehave in cases where all jets
// are active, because it would go to zero too fast)
result = sorted_quantities[0];
} else {
result = 0.0;
}
return result;
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
|
