FastJet 3.0.2
GridMedianBackgroundEstimator.cc
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends