FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: GridMedianBackgroundEstimator.cc 2689 2011-11-14 14:51:06Z soyez $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 00030 #include "fastjet/tools/GridMedianBackgroundEstimator.hh" 00031 using namespace std; 00032 00033 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00034 00035 //---------------------------------------------------------------------- 00036 // setting a new event 00037 //---------------------------------------------------------------------- 00038 // tell the background estimator that it has a new event, composed 00039 // of the specified particles. 00040 void GridMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) { 00041 fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0); 00042 for (unsigned i = 0; i < particles.size(); i++) { 00043 int j = igrid(particles[i]); 00044 if (j >= 0){ 00045 if (_rescaling_class == 0) 00046 _scalar_pt[j] += particles[i].perp(); 00047 else 00048 _scalar_pt[j] += particles[i].perp()/(*_rescaling_class)(particles[i]); 00049 } 00050 } 00051 sort(_scalar_pt.begin(), _scalar_pt.end()); 00052 00053 _has_particles = true; 00054 } 00055 00056 00057 //---------------------------------------------------------------------- 00058 // retrieving fundamental information 00059 //---------------------------------------------------------------------- 00060 // get rho, the median background density per unit area 00061 double GridMedianBackgroundEstimator::rho() const { 00062 verify_particles_set(); 00063 return _percentile(_scalar_pt, 0.5) / _cell_area; 00064 } 00065 00066 00067 //---------------------------------------------------------------------- 00068 // get sigma, the background fluctuations per unit area; must be 00069 // multipled by sqrt(area) to get fluctuations for a region of a 00070 // given area. 00071 double GridMedianBackgroundEstimator::sigma() const{ 00072 verify_particles_set(); 00073 // watch out: by definition, our sigma is the standard deviation of 00074 // the pt density multiplied by the square root of the cell area 00075 return (_percentile(_scalar_pt, 0.5) - 00076 _percentile(_scalar_pt, (1.0-0.6827)/2.0) 00077 )/sqrt(_cell_area); 00078 } 00079 00080 //---------------------------------------------------------------------- 00081 // get rho, the background density per unit area, locally at the 00082 // position of a given jet. Note that this is not const, because a 00083 // user may then wish to query other aspects of the background that 00084 // could depend on the position of the jet last used for a rho(jet) 00085 // determination. 00086 double GridMedianBackgroundEstimator::rho(const PseudoJet & jet) { 00087 verify_particles_set(); 00088 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet); 00089 return rescaling*rho(); 00090 } 00091 00092 00093 //---------------------------------------------------------------------- 00094 // get sigma, the background fluctuations per unit area, locally at 00095 // the position of a given jet. As for rho(jet), it is non-const. 00096 double GridMedianBackgroundEstimator::sigma(const PseudoJet & jet){ 00097 verify_particles_set(); 00098 double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet); 00099 return rescaling*sigma(); 00100 } 00101 00102 //---------------------------------------------------------------------- 00103 // verify that particles have been set and throw an error if not 00104 void GridMedianBackgroundEstimator::verify_particles_set() const { 00105 if (!_has_particles) throw Error("GridMedianBackgroundEstimator::rho() or sigma() called without particles having been set"); 00106 } 00107 00108 00109 //---------------------------------------------------------------------- 00110 // description 00111 //---------------------------------------------------------------------- 00112 string GridMedianBackgroundEstimator::description() const { 00113 ostringstream desc; 00114 desc << "GridMedianBackgroundEstimator, with grid extension |y| < " << _ymax 00115 << " and requested grid spacing = " << _requested_grid_spacing; 00116 return desc.str(); 00117 } 00118 00119 00120 //---------------------------------------------------------------------- 00121 // configuring the behaviour 00122 //---------------------------------------------------------------------- 00123 // Set a pointer to a class that calculates the rescaling factor as 00124 // a function of the jet (position). Note that the rescaling factor 00125 // is used both in the determination of the "global" rho (the pt/A 00126 // of each jet is divided by this factor) and when asking for a 00127 // local rho (the result is multiplied by this factor). 00128 // 00129 // The BackgroundRescalingYPolynomial class can be used to get a 00130 // rescaling that depends just on rapidity. 00131 // 00132 // Note that this has to be called BEFORE any attempt to do an 00133 // actual computation 00134 void GridMedianBackgroundEstimator::set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) { 00135 // The rescaling is taken into account when particles are set. So 00136 // you need to call set_particles again if you set the rescaling 00137 // class. We thus warn if there are already some available 00138 // particles 00139 if (_has_particles) 00140 _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."); 00141 00142 BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in); 00143 } 00144 00145 00146 //---------------------------------------------------------------------- 00147 // protected material 00148 //---------------------------------------------------------------------- 00149 // configure the grid 00150 void GridMedianBackgroundEstimator::setup_grid() { 00151 00152 // since we've exchanged the arguments of the grid constructor, 00153 // there's a danger of calls with exchanged ymax,spacing arguments -- 00154 // the following check should catch most such situations. 00155 assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing); 00156 00157 // this grid-definition code is becoming repetitive -- it should 00158 // probably be moved somewhere central... 00159 double ny_double = (_ymax-_ymin) / _requested_grid_spacing; 00160 _ny = int(ny_double+0.5); 00161 _dy = (_ymax-_ymin) / _ny; 00162 00163 _nphi = int (twopi / _requested_grid_spacing + 0.5); 00164 _dphi = twopi / _nphi; 00165 00166 // some sanity checking (could throw a fastjet::Error) 00167 assert(_ny >= 1 && _nphi >= 1); 00168 00169 _ntotal = _nphi * _ny; 00170 _scalar_pt.resize(_ntotal); 00171 _cell_area = _dy * _dphi; 00172 } 00173 00174 00175 //---------------------------------------------------------------------- 00176 // retrieve the grid cell index for a given PseudoJet 00177 int GridMedianBackgroundEstimator::igrid(const PseudoJet & p) const { 00178 // directly taking int does not work for values between -1 and 0 00179 // so use floor instead 00180 // double iy_double = (p.rap() - _ymin) / _dy; 00181 // if (iy_double < 0.0) return -1; 00182 // int iy = int(iy_double); 00183 // if (iy >= _ny) return -1; 00184 00185 // writing it as below gives a huge speed gain (factor two!). Even 00186 // though answers are identical and the routine here is not the 00187 // speed-critical step. It's not at all clear why. 00188 int iy = int(floor( (p.rap() - _ymin) / _dy )); 00189 if (iy < 0 || iy >= _ny) return -1; 00190 00191 int iphi = int( p.phi()/_dphi ); 00192 assert(iphi >= 0 && iphi <= _nphi); 00193 if (iphi == _nphi) iphi = 0; // just in case of rounding errors 00194 00195 int igrid_res = iy*_nphi + iphi; 00196 assert (igrid_res >= 0 && igrid_res < _ny*_nphi); 00197 return igrid_res; 00198 } 00199 00200 00201 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh