Fork me on GitHub

source: git/external/fastjet/tools/GridMedianBackgroundEstimator.cc@ 150ff64

Last change on this file since 150ff64 was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 12.1 KB
Line 
1//FJSTARTHEADER
2// $Id: GridMedianBackgroundEstimator.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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. They are described in the original FastJet paper,
16// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17// FastJet as part of work towards a scientific publication, please
18// quote the version you use and include a citation to the manual and
19// optionally also to hep-ph/0512210.
20//
21// FastJet is distributed in the hope that it will be useful,
22// but WITHOUT ANY WARRANTY; without even the implied warranty of
23// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24// GNU General Public License for more details.
25//
26// You should have received a copy of the GNU General Public License
27// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28//----------------------------------------------------------------------
29//FJENDHEADER
30
31
32#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
33using namespace std;
34
35FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
36
37
38//----------------------------------------------------------------------
39// setting a new event
40//----------------------------------------------------------------------
41// tell the background estimator that it has a new event, composed
42// of the specified particles.
43void GridMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
44 vector<double> scalar_pt(n_tiles(), 0.0);
45
46#ifdef FASTJET_GMBGE_USEFJGRID
47 assert(all_tiles_equal_area());
48 //assert(n_good_tiles() == n_tiles()); // not needed now that we have an implementation
49#endif
50
51 // check if we need to compute only rho or both rho and rho_m
52 if (_enable_rho_m){
53 // both rho and rho_m
54 //
55 // this requires a few other variables
56 vector<double> scalar_dt(n_tiles(), 0.0);
57 double pt, dt;
58 for (unsigned i = 0; i < particles.size(); i++) {
59 int j = tile_index(particles[i]);
60 if (j >= 0){
61 pt = particles[i].pt();
62 dt = particles[i].mt() - pt;
63 if (_rescaling_class == 0){
64 scalar_pt[j] += pt;
65 scalar_dt[j] += dt;
66 } else {
67 double r = (*_rescaling_class)(particles[i]);
68 scalar_pt[j] += pt/r;
69 scalar_dt[j] += dt/r;
70 }
71 }
72 }
73 // sort things for _percentile
74 sort(scalar_dt.begin(), scalar_dt.end());
75
76 // compute rho_m and sigma_m (see comment below for the
77 // normaliosation of sigma)
78 double p50 = _percentile(scalar_dt, 0.5);
79 _rho_m = p50 / mean_tile_area();
80 _sigma_m = (p50-_percentile(scalar_dt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area());
81 } else {
82 // only rho
83 //fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0);
84 for (unsigned i = 0; i < particles.size(); i++) {
85 int j = tile_index(particles[i]);
86 if (j >= 0){
87 if (_rescaling_class == 0){
88 scalar_pt[j] += particles[i].pt();
89 } else {
90 scalar_pt[j] += particles[i].pt()/(*_rescaling_class)(particles[i]);
91 }
92 }
93 }
94 }
95
96 // if there are some "bad" tiles, then we need to exclude them from
97 // the calculation of the median. We'll do this by condensing the
98 // scalar_pt vector down to just the values for the tiles that are
99 // good.
100 //
101 // tested answers look right in "issue" 2014-08-08-testing-rect-grid
102 if (n_good_tiles() != n_tiles()) {
103 int newn = 0;
104 for (unsigned i = 0; i < scalar_pt.size(); i++) {
105 if (tile_is_good(i)) {
106 // clang gets confused with the SharedPtr swap if we don't
107 // have std:: here
108 std::swap(scalar_pt[i],scalar_pt[newn]);
109 newn++;
110 }
111 }
112 scalar_pt.resize(newn);
113 }
114
115 // in all cases, carry on with the computation of rho
116 //
117 // first sort
118 sort(scalar_pt.begin(), scalar_pt.end());
119
120 // then compute rho
121 //
122 // watch out: by definition, our sigma is the standard deviation of
123 // the pt density multiplied by the square root of the cell area
124 double p50 = _percentile(scalar_pt, 0.5);
125 _rho = p50 / mean_tile_area();
126 _sigma = (p50-_percentile(scalar_pt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area());
127
128 _has_particles = true;
129}
130
131
132//----------------------------------------------------------------------
133// retrieving fundamental information
134//----------------------------------------------------------------------
135// get rho, the median background density per unit area
136double GridMedianBackgroundEstimator::rho() const {
137 verify_particles_set();
138 return _rho;
139}
140
141
142//----------------------------------------------------------------------
143// get sigma, the background fluctuations per unit area; must be
144// multipled by sqrt(area) to get fluctuations for a region of a
145// given area.
146double GridMedianBackgroundEstimator::sigma() const{
147 verify_particles_set();
148 return _sigma;
149}
150
151//----------------------------------------------------------------------
152// get rho, the background density per unit area, locally at the
153// position of a given jet. Note that this is not const, because a
154// user may then wish to query other aspects of the background that
155// could depend on the position of the jet last used for a rho(jet)
156// determination.
157double GridMedianBackgroundEstimator::rho(const PseudoJet & jet) {
158 //verify_particles_set();
159 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
160 return rescaling*rho();
161}
162
163
164//----------------------------------------------------------------------
165// get sigma, the background fluctuations per unit area, locally at
166// the position of a given jet. As for rho(jet), it is non-const.
167double GridMedianBackgroundEstimator::sigma(const PseudoJet & jet){
168 //verify_particles_set();
169 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
170 return rescaling*sigma();
171}
172
173//----------------------------------------------------------------------
174// returns rho_m (particle-masses contribution to the 4-vector density)
175double GridMedianBackgroundEstimator::rho_m() const {
176 if (! _enable_rho_m){
177 throw Error("GridMediamBackgroundEstimator: rho_m requested but rho_m calculation has been disabled.");
178 }
179 verify_particles_set();
180 return _rho_m;
181}
182
183
184//----------------------------------------------------------------------
185// returns sigma_m (particle-masses contribution to the 4-vector
186// density); must be multipled by sqrt(area) to get fluctuations
187// for a region of a given area.
188double GridMedianBackgroundEstimator::sigma_m() const{
189 if (! _enable_rho_m){
190 throw Error("GridMediamBackgroundEstimator: sigma_m requested but rho_m/sigma_m calculation has been disabled.");
191 }
192 verify_particles_set();
193 return _sigma_m;
194}
195
196//----------------------------------------------------------------------
197// returns rho_m locally at the position of a given jet. As for
198// rho(jet), it is non-const.
199double GridMedianBackgroundEstimator::rho_m(const PseudoJet & jet) {
200 //verify_particles_set();
201 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
202 return rescaling*rho_m();
203}
204
205
206//----------------------------------------------------------------------
207// returns sigma_m locally at the position of a given jet. As for
208// rho(jet), it is non-const.
209double GridMedianBackgroundEstimator::sigma_m(const PseudoJet & jet){
210 //verify_particles_set();
211 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
212 return rescaling*sigma_m();
213}
214
215//----------------------------------------------------------------------
216// verify that particles have been set and throw an error if not
217void GridMedianBackgroundEstimator::verify_particles_set() const {
218 if (!_has_particles) throw Error("GridMedianBackgroundEstimator::rho() or sigma() called without particles having been set");
219}
220
221
222//----------------------------------------------------------------------
223// description
224//----------------------------------------------------------------------
225string GridMedianBackgroundEstimator::description() const {
226 ostringstream desc;
227#ifdef FASTJET_GMBGE_USEFJGRID
228 desc << "GridMedianBackgroundEstimator, with " << RectangularGrid::description();
229#else
230 desc << "GridMedianBackgroundEstimator, with grid extension |y| < " << _ymax
231 << ", and grid cells of size dy x dphi = " << _dy << " x " << _dphi
232 << " (requested size = " << _requested_grid_spacing << ")";
233#endif
234 return desc.str();
235}
236
237
238//----------------------------------------------------------------------
239// configuring the behaviour
240//----------------------------------------------------------------------
241// Set a pointer to a class that calculates the rescaling factor as
242// a function of the jet (position). Note that the rescaling factor
243// is used both in the determination of the "global" rho (the pt/A
244// of each jet is divided by this factor) and when asking for a
245// local rho (the result is multiplied by this factor).
246//
247// The BackgroundRescalingYPolynomial class can be used to get a
248// rescaling that depends just on rapidity.
249//
250// Note that this has to be called BEFORE any attempt to do an
251// actual computation
252void GridMedianBackgroundEstimator::set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) {
253 // The rescaling is taken into account when particles are set. So
254 // you need to call set_particles again if you set the rescaling
255 // class. We thus warn if there are already some available
256 // particles
257 if (_has_particles)
258 _warning_rescaling.warn("GridMedianBackgroundEstimator::set_rescaling_class(): trying to set the rescaling class when there are already particles that have been set is dangerous: the rescaling will not affect the already existing particles resulting in mis-estimation of rho. You need to call set_particles() again before proceeding with any background estimation.");
259
260 BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
261}
262
263
264#ifndef FASTJET_GMBGE_USEFJGRID
265//----------------------------------------------------------------------
266// protected material
267//----------------------------------------------------------------------
268// configure the grid
269void GridMedianBackgroundEstimator::setup_grid() {
270
271 // since we've exchanged the arguments of the grid constructor,
272 // there's a danger of calls with exchanged ymax,spacing arguments --
273 // the following check should catch most such situations.
274 assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing);
275
276 // this grid-definition code is becoming repetitive -- it should
277 // probably be moved somewhere central...
278 double ny_double = (_ymax-_ymin) / _requested_grid_spacing;
279 _ny = int(ny_double+0.5);
280 _dy = (_ymax-_ymin) / _ny;
281
282 _nphi = int (twopi / _requested_grid_spacing + 0.5);
283 _dphi = twopi / _nphi;
284
285 // some sanity checking (could throw a fastjet::Error)
286 assert(_ny >= 1 && _nphi >= 1);
287
288 _ntotal = _nphi * _ny;
289 //_scalar_pt.resize(_ntotal);
290 _tile_area = _dy * _dphi;
291}
292
293
294//----------------------------------------------------------------------
295// retrieve the grid tile index for a given PseudoJet
296int GridMedianBackgroundEstimator::tile_index(const PseudoJet & p) const {
297 // directly taking int does not work for values between -1 and 0
298 // so use floor instead
299 // double iy_double = (p.rap() - _ymin) / _dy;
300 // if (iy_double < 0.0) return -1;
301 // int iy = int(iy_double);
302 // if (iy >= _ny) return -1;
303
304 // writing it as below gives a huge speed gain (factor two!). Even
305 // though answers are identical and the routine here is not the
306 // speed-critical step. It's not at all clear why.
307 int iy = int(floor( (p.rap() - _ymin) / _dy ));
308 if (iy < 0 || iy >= _ny) return -1;
309
310 int iphi = int( p.phi()/_dphi );
311 assert(iphi >= 0 && iphi <= _nphi);
312 if (iphi == _nphi) iphi = 0; // just in case of rounding errors
313
314 int index_res = iy*_nphi + iphi;
315 assert (index_res >= 0 && index_res < _ny*_nphi);
316 return index_res;
317}
318#endif // FASTJET_GMBGE_USEFJGRID
319
320
321
322FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.