1 | //STARTHEADER |
---|
2 | // $Id: GridMedianBackgroundEstimator.cc 999 2013-03-04 11:48:06Z pavel $ |
---|
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 |
---|
61 | double GridMedianBackgroundEstimator::rho() const { |
---|
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. |
---|
71 | double GridMedianBackgroundEstimator::sigma() const{ |
---|
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 | //---------------------------------------------------------------------- |
---|
112 | string GridMedianBackgroundEstimator::description() const { |
---|
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 |
---|