FastJet  3.3.0
GridMedianBackgroundEstimator.cc
1 //FJSTARTHEADER
2 // $Id: GridMedianBackgroundEstimator.cc 3555 2014-08-11 09:56:35Z salam $
3 //
4 // Copyright (c) 2005-2014, 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"
33 using namespace std;
34 
35 FASTJET_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.
43 void 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
136 double 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.
146 double 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.
157 double 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.
167 double 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)
175 double 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.
188 double 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.
199 double 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.
209 double 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
217 void 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 //----------------------------------------------------------------------
225 string 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
252 void 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
269 void 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
296 int 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 
322 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
double sigma(const PseudoJet &jet)
returns sigma, the background fluctuations per unit area, locally at the position of a given jet...
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67