1 #ifndef __GRID_MEDIAN_BACKGROUND_ESTIMATOR_HH__ 2 #define __GRID_MEDIAN_BACKGROUND_ESTIMATOR_HH__ 3 4 //FJSTARTHEADER 5 // $Id: GridMedianBackgroundEstimator.hh 4442 2020-05-05 07:50:11Z soyez $ 6 // 7 // Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 8 // 9 //---------------------------------------------------------------------- 10 // This file is part of FastJet. 11 // 12 // FastJet is free software; you can redistribute it and/or modify 13 // it under the terms of the GNU General Public License as published by 14 // the Free Software Foundation; either version 2 of the License, or 15 // (at your option) any later version. 16 // 17 // The algorithms that underlie FastJet have required considerable 18 // development. They are described in the original FastJet paper, 19 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use 20 // FastJet as part of work towards a scientific publication, please 21 // quote the version you use and include a citation to the manual and 22 // optionally also to hep-ph/0512210. 23 // 24 // FastJet is distributed in the hope that it will be useful, 25 // but WITHOUT ANY WARRANTY; without even the implied warranty of 26 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27 // GNU General Public License for more details. 28 // 29 // You should have received a copy of the GNU General Public License 30 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 31 //---------------------------------------------------------------------- 32 //FJENDHEADER 33 34 35 #include "fastjet/tools/BackgroundEstimatorBase.hh" 36 37 // if defined then we'll use the RectangularGrid class 38 // 39 // (For FastJet 3.2, maybe remove the symbol and simply clean up the 40 // code below to use exclusively the RectangularGrid) 41 #define FASTJET_GMBGE_USEFJGRID 42 43 #ifdef FASTJET_GMBGE_USEFJGRID 44 #include "fastjet/RectangularGrid.hh" 45 #endif 46 47 48 49 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 50 51 /// @ingroup tools_background 52 /// \class GridMedianBackgroundEstimator 53 /// 54 /// Background Estimator based on the median pt/area of a set of grid 55 /// cells. 56 /// 57 /// Description of the method: 58 /// This background estimator works by projecting the event onto a 59 /// grid in rapidity and azimuth. In each grid cell, the scalar pt 60 /// sum of the particles in the cell is computed. The background 61 /// density is then estimated by the median of (scalar pt sum/cell 62 /// area) for all cells. 63 /// 64 /// Parameters: 65 /// The class takes 2 arguments: the absolute rapidity extent of the 66 /// cells and the size of the grid cells. Note that the size of the cell 67 /// will be adjusted in azimuth to satisfy the 2pi periodicity and 68 /// in rapidity to match the requested rapidity extent. 69 /// 70 /// Rescaling: 71 /// It is possible to use a rescaling profile. In this case, the 72 /// profile needs to be set before setting the particles and it will 73 /// be applied to each particle (i.e. not to each cell). 74 /// Note also that in this case one needs to call rho(jet) instead of 75 /// rho() [Without rescaling, they are identical] 76 /// 77 class GridMedianBackgroundEstimator : public BackgroundEstimatorBase 78 #ifdef FASTJET_GMBGE_USEFJGRID 79 , public RectangularGrid 80 #endif 81 { 82 83 public: 84 /// @name constructors and destructors 85 //\{ 86 #ifdef FASTJET_GMBGE_USEFJGRID 87 //---------------------------------------------------------------- 88 /// \param ymax maximal absolute rapidity extent of the grid 89 /// \param requested_grid_spacing size of the grid cell. The 90 /// "real" cell size could differ due e.g. to the 2pi 91 /// periodicity in azimuthal angle (size, not area) GridMedianBackgroundEstimator(double ymax,double requested_grid_spacing)92 GridMedianBackgroundEstimator(double ymax, double requested_grid_spacing) : 93 RectangularGrid(ymax, requested_grid_spacing), 94 _has_particles(false), _enable_rho_m(true) {} 95 96 //---------------------------------------------------------------- 97 /// Constructor based on a user's fully specified RectangularGrid GridMedianBackgroundEstimator(const RectangularGrid & grid)98 GridMedianBackgroundEstimator(const RectangularGrid & grid) : 99 RectangularGrid(grid), 100 _has_particles(false), _enable_rho_m(true) { 101 if (!RectangularGrid::is_initialised()) 102 throw Error("attempt to construct GridMedianBackgroundEstimator with uninitialised RectangularGrid"); 103 } 104 105 //---------------------------------------------------------------- 106 /// Constructor with the explicit parameters for the underlying 107 /// RectangularGrid 108 /// 109 /// \param rapmin the minimum rapidity extent of the grid 110 /// \param rapmax the maximum rapidity extent of the grid 111 /// \param drap the grid spacing in rapidity 112 /// \param dphi the grid spacing in azimuth 113 /// \param tile_selector optional (geometric) selector to specify 114 /// which tiles are good; a tile is good if 115 /// a massless 4-vector at the center of the tile passes 116 /// the selection GridMedianBackgroundEstimator(double rapmin_in,double rapmax_in,double drap_in,double dphi_in,Selector tile_selector=Selector ())117 GridMedianBackgroundEstimator(double rapmin_in, double rapmax_in, double drap_in, double dphi_in, 118 Selector tile_selector = Selector()) : 119 RectangularGrid(rapmin_in, rapmax_in, drap_in, dphi_in, tile_selector), 120 _has_particles(false), _enable_rho_m(true) {} 121 122 #else // alternative in old framework where we didn't have the rectangular grid 123 GridMedianBackgroundEstimator(double ymax, double requested_grid_spacing) : 124 _ymin(-ymax), _ymax(ymax), 125 _requested_grid_spacing(requested_grid_spacing), 126 _has_particles(false), _enable_rho_m(true) 127 { 128 setup_grid(); 129 } 130 #endif // FASTJET_GMBGE_USEFJGRID 131 132 //\} 133 134 135 /// @name setting a new event 136 //\{ 137 //---------------------------------------------------------------- 138 139 /// tell the background estimator that it has a new event, composed 140 /// of the specified particles. 141 void set_particles(const std::vector<PseudoJet> & particles); 142 143 /// determine whether the automatic calculation of rho_m and sigma_m 144 /// is enabled (by default true) set_compute_rho_m(bool enable)145 void set_compute_rho_m(bool enable){ _enable_rho_m = enable;} 146 147 //\} 148 149 /// @name retrieving fundamental information 150 //\{ 151 //---------------------------------------------------------------- 152 153 /// returns rho, the median background density per unit area 154 double rho() const; 155 156 /// returns sigma, the background fluctuations per unit area; must be 157 /// multipled by sqrt(area) to get fluctuations for a region of a 158 /// given area. 159 double sigma() const; 160 161 /// returns rho, the background density per unit area, locally at the 162 /// position of a given jet. Note that this is not const, because a 163 /// user may then wish to query other aspects of the background that 164 /// could depend on the position of the jet last used for a rho(jet) 165 /// determination. 166 double rho(const PseudoJet & jet); 167 168 /// returns sigma, the background fluctuations per unit area, locally at 169 /// the position of a given jet. As for rho(jet), it is non-const. 170 double sigma(const PseudoJet & jet); 171 172 /// returns true if this background estimator has support for 173 /// determination of sigma has_sigma()174 bool has_sigma() {return true;} 175 176 //----------------------------------------------------------------- 177 /// Returns rho_m, the purely longitudinal, particle-mass-induced 178 /// component of the background density per unit area 179 double rho_m() const; 180 181 /// returns sigma_m, a measure of the fluctuations in the purely 182 /// longitudinal, particle-mass-induced component of the background 183 /// density per unit area; must be multipled by sqrt(area) to get 184 /// fluctuations for a region of a given area. 185 double sigma_m() const; 186 187 /// Returns rho_m locally at the jet position. As for rho(jet), it is non-const. 188 double rho_m(const PseudoJet & jet); 189 190 /// Returns sigma_m locally at the jet position. As for rho(jet), it is non-const. 191 double sigma_m(const PseudoJet & jet); 192 193 /// Returns true if this background estimator has support for 194 /// determination of rho_m. 195 /// 196 /// Note that support for sigma_m is automatic if one has sigma and 197 /// rho_m support. has_rho_m() const198 bool has_rho_m() const {return _enable_rho_m;} 199 200 201 /// returns the area of the grid cells (all identical, but 202 /// referred to as "mean" area for uniformity with JetMedianBGE). mean_area() const203 double mean_area() const {return mean_tile_area();} 204 //\} 205 206 /// @name configuring the behaviour 207 //\{ 208 //---------------------------------------------------------------- 209 210 /// Set a pointer to a class that calculates the rescaling factor as 211 /// a function of the jet (position). Note that the rescaling factor 212 /// is used both in the determination of the "global" rho (the pt/A 213 /// of each jet is divided by this factor) and when asking for a 214 /// local rho (the result is multiplied by this factor). 215 /// 216 /// The BackgroundRescalingYPolynomial class can be used to get a 217 /// rescaling that depends just on rapidity. 218 /// 219 /// Note that this has to be called BEFORE any attempt to do an 220 /// actual computation 221 /// 222 /// The same profile will be used for both pt and mt (this is 223 /// probabaly a good approximation since the particle density 224 /// changes is what dominates the rapidity profile) 225 virtual void set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class); 226 227 //\} 228 229 /// @name description 230 //\{ 231 //---------------------------------------------------------------- 232 233 /// returns a textual description of the background estimator 234 std::string description() const; 235 236 //\} 237 238 239 private: 240 241 #ifndef FASTJET_GMBGE_USEFJGRID 242 243 /// configure the grid 244 void setup_grid(); 245 246 /// retrieve the grid cell index for a given PseudoJet 247 int tile_index(const PseudoJet & p) const; 248 249 // information about the grid 250 double _ymin, _ymax, _dy, _dphi, _requested_grid_spacing, _tile_area; 251 int _ny, _nphi, _ntotal; 252 n_tiles() const253 int n_tiles() const {return _ntotal;} n_good_tiles() const254 int n_good_tiles() const {return n_tiles();} tile_is_good(int) const255 int tile_is_good(int /* itile */) const {return true;} 256 mean_tile_area() const257 double mean_tile_area() const {return _tile_area;} 258 #endif // FASTJET_GMBGE_USEFJGRID 259 260 261 /// verify that particles have been set and throw an error if not 262 void verify_particles_set() const; 263 264 // information abotu the event 265 //std::vector<double> _scalar_pt; 266 double _rho, _sigma, _rho_m, _sigma_m; 267 bool _has_particles; 268 bool _enable_rho_m; 269 270 // various warnings to inform people of potential dangers 271 LimitedWarning _warning_rho_of_jet; 272 LimitedWarning _warning_rescaling; 273 }; 274 275 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh 276 277 #endif // __GRID_MEDIAN_BACKGROUND_ESTIMATOR_HH__ 278