Fork me on GitHub

source: svn/trunk/external/fastjet/tools/GridMedianBackgroundEstimator.cc@ 1204

Last change on this file since 1204 was 999, checked in by Pavel Demin, 12 years ago

add fastjet/tools

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 8.1 KB
Line 
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"
31using namespace std;
32
33FASTJET_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.
40void 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
61double 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.
71double 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.
86double 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.
96double 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
104void 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//----------------------------------------------------------------------
112string 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
134void 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
150void 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
177int 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
201FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.