30 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
33 FASTJET_BEGIN_NAMESPACE
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]);
45 if (_rescaling_class == 0)
46 _scalar_pt[j] += particles[i].perp();
48 _scalar_pt[j] += particles[i].perp()/(*_rescaling_class)(particles[i]);
51 sort(_scalar_pt.begin(), _scalar_pt.end());
53 _has_particles =
true;
62 verify_particles_set();
63 return _percentile(_scalar_pt, 0.5) / _cell_area;
72 verify_particles_set();
75 return (_percentile(_scalar_pt, 0.5) -
76 _percentile(_scalar_pt, (1.0-0.6827)/2.0)
87 verify_particles_set();
88 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
89 return rescaling*rho();
97 verify_particles_set();
98 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
99 return rescaling*sigma();
104 void GridMedianBackgroundEstimator::verify_particles_set()
const {
105 if (!_has_particles)
throw Error(
"GridMedianBackgroundEstimator::rho() or sigma() called without particles having been set");
114 desc <<
"GridMedianBackgroundEstimator, with grid extension |y| < " << _ymax
115 <<
" and requested grid spacing = " << _requested_grid_spacing;
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.");
142 BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
150 void GridMedianBackgroundEstimator::setup_grid() {
155 assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing);
159 double ny_double = (_ymax-_ymin) / _requested_grid_spacing;
160 _ny = int(ny_double+0.5);
161 _dy = (_ymax-_ymin) / _ny;
163 _nphi = int (twopi / _requested_grid_spacing + 0.5);
164 _dphi = twopi / _nphi;
167 assert(_ny >= 1 && _nphi >= 1);
169 _ntotal = _nphi * _ny;
170 _scalar_pt.resize(_ntotal);
171 _cell_area = _dy * _dphi;
177 int GridMedianBackgroundEstimator::igrid(
const PseudoJet & p)
const {
188 int iy = int(floor( (p.rap() - _ymin) / _dy ));
189 if (iy < 0 || iy >= _ny)
return -1;
191 int iphi = int( p.phi()/_dphi );
192 assert(iphi >= 0 && iphi <= _nphi);
193 if (iphi == _nphi) iphi = 0;
195 int igrid_res = iy*_nphi + iphi;
196 assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
201 FASTJET_END_NAMESPACE