Loading...
Searching...
No Matches
GridMedianBackgroundEstimator.cc
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
//----------------------------------------------------------------------
// setting a new event
//----------------------------------------------------------------------
// tell the background estimator that it has a new event, composed
// of the specified particles.
void GridMedianBackgroundEstimator::set_particles(const vector<PseudoJet> & particles) {
vector<double> scalar_pt(n_tiles(), 0.0);
assert(all_tiles_equal_area());
//assert(n_good_tiles() == n_tiles()); // not needed now that we have an implementation
_cached_estimate.reset();
_cached_estimate.set_has_sigma(true);
_cached_estimate.set_mean_area(mean_tile_area());
// check if we need to compute only rho or both rho and rho_m
if (_enable_rho_m){
// both rho and rho_m
//
// this requires a few other variables
vector<double> scalar_dt(n_tiles(), 0.0);
double pt, dt;
for (unsigned i = 0; i < particles.size(); i++) {
int j = tile_index(particles[i]);
if (j >= 0){
pt = particles[i].pt();
dt = particles[i].mt() - pt;
if (_rescaling_class == 0){
scalar_pt[j] += pt;
scalar_dt[j] += dt;
} else {
double r = (*_rescaling_class)(particles[i]);
scalar_pt[j] += pt/r;
scalar_dt[j] += dt/r;
}
}
}
// sort things for _percentile
sort(scalar_dt.begin(), scalar_dt.end());
// compute rho_m and sigma_m (see comment below for the
// normaliosation of sigma)
double p50 = _percentile(scalar_dt, 0.5);
_cached_estimate.set_has_rho_m(true);
_cached_estimate.set_rho_m(p50 / mean_tile_area());
_cached_estimate.set_sigma_m((p50-_percentile(scalar_dt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area()));
} else {
// only rho
//fill(_scalar_pt.begin(), _scalar_pt.end(), 0.0);
for (unsigned i = 0; i < particles.size(); i++) {
int j = tile_index(particles[i]);
if (j >= 0){
if (_rescaling_class == 0){
scalar_pt[j] += particles[i].pt();
} else {
scalar_pt[j] += particles[i].pt()/(*_rescaling_class)(particles[i]);
}
}
}
}
// if there are some "bad" tiles, then we need to exclude them from
// the calculation of the median. We'll do this by condensing the
// scalar_pt vector down to just the values for the tiles that are
// good.
//
// tested answers look right in "issue" 2014-08-08-testing-rect-grid
if (n_good_tiles() != n_tiles()) {
int newn = 0;
for (unsigned i = 0; i < scalar_pt.size(); i++) {
if (tile_is_good(i)) {
// clang gets confused with the SharedPtr swap if we don't
// have std:: here
std::swap(scalar_pt[i],scalar_pt[newn]);
newn++;
}
}
scalar_pt.resize(newn);
}
// in all cases, carry on with the computation of rho
//
// first sort
sort(scalar_pt.begin(), scalar_pt.end());
// then compute rho
//
// watch out: by definition, our sigma is the standard deviation of
// the pt density multiplied by the square root of the cell area
double p50 = _percentile(scalar_pt, 0.5);
_cached_estimate.set_rho(p50 / mean_tile_area());
_cached_estimate.set_sigma((p50-_percentile(scalar_pt, (1.0-0.6827)/2.0))/sqrt(mean_tile_area()));
_cache_available = true;
}
//----------------------------------------------------------------------
// retrieving fundamental information
//----------------------------------------------------------------------
// get the full set of background properties
BackgroundEstimate GridMedianBackgroundEstimator::estimate() const{
verify_particles_set();
return _cached_estimate;
}
// get the full set of background properties for a given reference jet
BackgroundEstimate GridMedianBackgroundEstimator::estimate(const PseudoJet &jet) const{
verify_particles_set();
if (_rescaling_class == 0)
return _cached_estimate;
BackgroundEstimate local_estimate = _cached_estimate;
local_estimate.apply_rescaling_factor((*_rescaling_class)(jet));
return local_estimate;
}
// get rho, the median background density per unit area
double GridMedianBackgroundEstimator::rho() const {
verify_particles_set();
return _cached_estimate.rho();
}
//----------------------------------------------------------------------
// get sigma, the background fluctuations per unit area; must be
// multipled by sqrt(area) to get fluctuations for a region of a
// given area.
double GridMedianBackgroundEstimator::sigma() const{
verify_particles_set();
return _cached_estimate.sigma();
}
//----------------------------------------------------------------------
// get rho, the background density per unit area, locally at the
// position of a given jet. Note that this is not const, because a
// user may then wish to query other aspects of the background that
// could depend on the position of the jet last used for a rho(jet)
// determination.
double GridMedianBackgroundEstimator::rho(const PseudoJet & jet) {
//verify_particles_set();
double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
return rescaling*rho();
}
//----------------------------------------------------------------------
// get sigma, the background fluctuations per unit area, locally at
// the position of a given jet. As for rho(jet), it is non-const.
double GridMedianBackgroundEstimator::sigma(const PseudoJet & jet){
//verify_particles_set();
double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
return rescaling*sigma();
}
//----------------------------------------------------------------------
// returns rho_m (particle-masses contribution to the 4-vector density)
double GridMedianBackgroundEstimator::rho_m() const {
if (! _enable_rho_m){
throw Error("GridMediamBackgroundEstimator: rho_m requested but rho_m calculation has been disabled.");
}
verify_particles_set();
return _cached_estimate.rho_m();
}
//----------------------------------------------------------------------
// returns sigma_m (particle-masses contribution to the 4-vector
// density); must be multipled by sqrt(area) to get fluctuations
// for a region of a given area.
double GridMedianBackgroundEstimator::sigma_m() const{
if (! _enable_rho_m){
throw Error("GridMediamBackgroundEstimator: sigma_m requested but rho_m/sigma_m calculation has been disabled.");
}
verify_particles_set();
return _cached_estimate.sigma_m();
}
//----------------------------------------------------------------------
// returns rho_m locally at the position of a given jet. As for
// rho(jet), it is non-const.
double GridMedianBackgroundEstimator::rho_m(const PseudoJet & jet) {
//verify_particles_set();
double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
return rescaling*rho_m();
}
//----------------------------------------------------------------------
// returns sigma_m locally at the position of a given jet. As for
// rho(jet), it is non-const.
double GridMedianBackgroundEstimator::sigma_m(const PseudoJet & jet){
//verify_particles_set();
double rescaling = (_rescaling_class == 0) ? 1.0 : (*_rescaling_class)(jet);
return rescaling*sigma_m();
}
//----------------------------------------------------------------------
// verify that particles have been set and throw an error if not
void GridMedianBackgroundEstimator::verify_particles_set() const {
if (!_cache_available) throw Error("GridMedianBackgroundEstimator::rho() or sigma() called without particles having been set");
}
//----------------------------------------------------------------------
// description
//----------------------------------------------------------------------
string GridMedianBackgroundEstimator::description() const {
ostringstream desc;
desc << "GridMedianBackgroundEstimator, with " << RectangularGrid::description();
return desc.str();
}
//----------------------------------------------------------------------
// configuring the behaviour
//----------------------------------------------------------------------
// Set a pointer to a class that calculates the rescaling factor as
// a function of the jet (position). Note that the rescaling factor
// is used both in the determination of the "global" rho (the pt/A
// of each jet is divided by this factor) and when asking for a
// local rho (the result is multiplied by this factor).
//
// The BackgroundRescalingYPolynomial class can be used to get a
// rescaling that depends just on rapidity.
//
// Note that this has to be called BEFORE any attempt to do an
// actual computation
void GridMedianBackgroundEstimator::set_rescaling_class(const FunctionOfPseudoJet<double> * rescaling_class_in) {
// The rescaling is taken into account when particles are set. So
// you need to call set_particles again if you set the rescaling
// class. We thus warn if there are already some available
// particles
if (_cache_available)
_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.");
BackgroundEstimatorBase::set_rescaling_class(rescaling_class_in);
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::BackgroundEstimate
/// a class that holds the result of the calculation
Definition BackgroundEstimatorBase.hh:53
fastjet::BackgroundEstimate::apply_rescaling_factor
void apply_rescaling_factor(double rescaling_factor)
apply a rescaling factor (to rho, rho_m, sigma, sigma_m)
Definition BackgroundEstimatorBase.hh:144
fastjet::BackgroundEstimatorBase::set_rescaling_class
virtual void set_rescaling_class(const FunctionOfPseudoJet< double > *rescaling_class_in)
Set a pointer to a class that calculates the rescaling factor as a function of the jet (position).
Definition BackgroundEstimatorBase.hh:311
fastjet::BackgroundEstimatorBase::_percentile
double _percentile(const std::vector< double > &sorted_quantity_vector, const double percentile, const double nempty=0.0, const bool do_fj2_calculation=false) const
computes a percentile of a given sorted vector
Definition BackgroundEstimatorBase.cc:104
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::FunctionOfPseudoJet
base class providing interface for a generic function of a PseudoJet
Definition FunctionOfPseudoJet.hh:48
fastjet::GridMedianBackgroundEstimator::sigma_m
double sigma_m() const override
returns sigma_m, a measure of the fluctuations in the purely longitudinal, particle-mass-induced comp...
Definition GridMedianBackgroundEstimator.cc:210
fastjet::GridMedianBackgroundEstimator::estimate
BackgroundEstimate estimate() const override
get the full set of background properties
Definition GridMedianBackgroundEstimator.cc:140
fastjet::GridMedianBackgroundEstimator::sigma
double sigma() const override
returns sigma, the background fluctuations per unit area; must be multipled by sqrt(area) to get fluc...
Definition GridMedianBackgroundEstimator.cc:168
fastjet::GridMedianBackgroundEstimator::set_rescaling_class
virtual void set_rescaling_class(const FunctionOfPseudoJet< double > *rescaling_class) override
Set a pointer to a class that calculates the rescaling factor as a function of the jet (position).
Definition GridMedianBackgroundEstimator.cc:268
fastjet::GridMedianBackgroundEstimator::set_particles
void set_particles(const std::vector< PseudoJet > &particles) override
tell the background estimator that it has a new event, composed of the specified particles.
Definition GridMedianBackgroundEstimator.cc:43
fastjet::GridMedianBackgroundEstimator::rho_m
double rho_m() const override
Returns rho_m, the purely longitudinal, particle-mass-induced component of the background density per...
Definition GridMedianBackgroundEstimator.cc:197
fastjet::GridMedianBackgroundEstimator::rho
double rho() const override
returns rho, the median background density per unit area
Definition GridMedianBackgroundEstimator.cc:158
fastjet::GridMedianBackgroundEstimator::description
std::string description() const override
returns a textual description of the background estimator
Definition GridMedianBackgroundEstimator.cc:247
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::RectangularGrid::mean_tile_area
virtual double mean_tile_area() const override
returns the mean area of tiles.
Definition RectangularGrid.hh:140
fastjet::RectangularGrid::tile_index
virtual int tile_index(const PseudoJet &p) const override
returns the index of the tile in which p is located, or -1 if p is outside the tiling region
Definition RectangularGrid.cc:55
fastjet::RectangularGrid::tile_is_good
virtual bool tile_is_good(int itile) const override
returns whether a given tile is good
Definition RectangularGrid.hh:130
fastjet::RectangularGrid::n_tiles
virtual int n_tiles() const override
returns the total number of tiles in the tiling; valid tile indices run from 0 ......
Definition RectangularGrid.hh:120
fastjet::RectangularGrid::description
virtual std::string description() const override
returns a textual description of the grid
Definition RectangularGrid.cc:121
fastjet::RectangularGrid::n_good_tiles
virtual int n_good_tiles() const override
returns the number of tiles that are "good"; i.e.
Definition RectangularGrid.hh:122
fastjet::TilingBase::all_tiles_equal_area
virtual bool all_tiles_equal_area() const
returns true if all tiles have the same area
Definition RectangularGrid.hh:65