FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
GridMedianBackgroundEstimator.cc
1 //STARTHEADER
2 // $Id: GridMedianBackgroundEstimator.cc 2689 2011-11-14 14:51:06Z soyez $
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 
30 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
31 using namespace std;
32 
33 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
34 
35 //----------------------------------------------------------------------
36 // setting a new event
37 //----------------------------------------------------------------------
38 // tell the background estimator that it has a new event, composed
39 // of the specified particles.
40 void GridMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
41  fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0);
42  for (unsigned i = 0; i < particles.size(); i++) {
43  int j = igrid(particles[i]);
44  if (j >= 0){
45  if (_rescaling_class == 0)
46  _scalar_pt[j] += particles[i].perp();
47  else
48  _scalar_pt[j] += particles[i].perp()/(*_rescaling_class)(particles[i]);
49  }
50  }
51  sort(_scalar_pt.begin(), _scalar_pt.end());
52 
53  _has_particles = true;
54 }
55 
56 
57 //----------------------------------------------------------------------
58 // retrieving fundamental information
59 //----------------------------------------------------------------------
60 // get rho, the median background density per unit area
62  verify_particles_set();
63  return _percentile(_scalar_pt, 0.5) / _cell_area;
64 }
65 
66 
67 //----------------------------------------------------------------------
68 // get sigma, the background fluctuations per unit area; must be
69 // multipled by sqrt(area) to get fluctuations for a region of a
70 // given area.
72  verify_particles_set();
73  // watch out: by definition, our sigma is the standard deviation of
74  // the pt density multiplied by the square root of the cell area
75  return (_percentile(_scalar_pt, 0.5) -
76  _percentile(_scalar_pt, (1.0-0.6827)/2.0)
77  )/sqrt(_cell_area);
78 }
79 
80 //----------------------------------------------------------------------
81 // get rho, the background density per unit area, locally at the
82 // position of a given jet. Note that this is not const, because a
83 // user may then wish to query other aspects of the background that
84 // could depend on the position of the jet last used for a rho(jet)
85 // determination.
86 double GridMedianBackgroundEstimator::rho(const PseudoJet & jet) {
87  verify_particles_set();
88  double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
89  return rescaling*rho();
90 }
91 
92 
93 //----------------------------------------------------------------------
94 // get sigma, the background fluctuations per unit area, locally at
95 // the position of a given jet. As for rho(jet), it is non-const.
96 double GridMedianBackgroundEstimator::sigma(const PseudoJet & jet){
97  verify_particles_set();
98  double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
99  return rescaling*sigma();
100 }
101 
102 //----------------------------------------------------------------------
103 // verify that particles have been set and throw an error if not
104 void GridMedianBackgroundEstimator::verify_particles_set() const {
105  if (!_has_particles) throw Error("GridMedianBackgroundEstimator::rho() or sigma() called without particles having been set");
106 }
107 
108 
109 //----------------------------------------------------------------------
110 // description
111 //----------------------------------------------------------------------
113  ostringstream desc;
114  desc << "GridMedianBackgroundEstimator, with grid extension |y| < " << _ymax
115  << " and requested grid spacing = " << _requested_grid_spacing;
116  return desc.str();
117 }
118 
119 
120 //----------------------------------------------------------------------
121 // configuring the behaviour
122 //----------------------------------------------------------------------
123 // Set a pointer to a class that calculates the rescaling factor as
124 // a function of the jet (position). Note that the rescaling factor
125 // is used both in the determination of the "global" rho (the pt/A
126 // of each jet is divided by this factor) and when asking for a
127 // local rho (the result is multiplied by this factor).
128 //
129 // The BackgroundRescalingYPolynomial class can be used to get a
130 // rescaling that depends just on rapidity.
131 //
132 // Note that this has to be called BEFORE any attempt to do an
133 // actual computation
134 void GridMedianBackgroundEstimator::set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) {
135  // The rescaling is taken into account when particles are set. So
136  // you need to call set_particles again if you set the rescaling
137  // class. We thus warn if there are already some available
138  // particles
139  if (_has_particles)
140  _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.");
141 
142  BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
143 }
144 
145 
146 //----------------------------------------------------------------------
147 // protected material
148 //----------------------------------------------------------------------
149 // configure the grid
150 void GridMedianBackgroundEstimator::setup_grid() {
151 
152  // since we've exchanged the arguments of the grid constructor,
153  // there's a danger of calls with exchanged ymax,spacing arguments --
154  // the following check should catch most such situations.
155  assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing);
156 
157  // this grid-definition code is becoming repetitive -- it should
158  // probably be moved somewhere central...
159  double ny_double = (_ymax-_ymin) / _requested_grid_spacing;
160  _ny = int(ny_double+0.5);
161  _dy = (_ymax-_ymin) / _ny;
162 
163  _nphi = int (twopi / _requested_grid_spacing + 0.5);
164  _dphi = twopi / _nphi;
165 
166  // some sanity checking (could throw a fastjet::Error)
167  assert(_ny >= 1 && _nphi >= 1);
168 
169  _ntotal = _nphi * _ny;
170  _scalar_pt.resize(_ntotal);
171  _cell_area = _dy * _dphi;
172 }
173 
174 
175 //----------------------------------------------------------------------
176 // retrieve the grid cell index for a given PseudoJet
177 int GridMedianBackgroundEstimator::igrid(const PseudoJet & p) const {
178  // directly taking int does not work for values between -1 and 0
179  // so use floor instead
180  // double iy_double = (p.rap() - _ymin) / _dy;
181  // if (iy_double < 0.0) return -1;
182  // int iy = int(iy_double);
183  // if (iy >= _ny) return -1;
184 
185  // writing it as below gives a huge speed gain (factor two!). Even
186  // though answers are identical and the routine here is not the
187  // speed-critical step. It's not at all clear why.
188  int iy = int(floor( (p.rap() - _ymin) / _dy ));
189  if (iy < 0 || iy >= _ny) return -1;
190 
191  int iphi = int( p.phi()/_dphi );
192  assert(iphi >= 0 && iphi <= _nphi);
193  if (iphi == _nphi) iphi = 0; // just in case of rounding errors
194 
195  int igrid_res = iy*_nphi + iphi;
196  assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
197  return igrid_res;
198 }
199 
200 
201 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh