Loading...
Searching...
No Matches
BackgroundEstimatorBase.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/BackgroundEstimatorBase.hh"
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
LimitedWarning BackgroundEstimatorBase::_warnings_empty_area;
#ifdef FASTJET_HAVE_THREAD_SAFETY
BackgroundEstimatorBase::BackgroundEstimatorBase(const BackgroundEstimatorBase &other_bge){
_rescaling_class = other_bge._rescaling_class;
_cached_estimate = other_bge._cached_estimate;
_cache_available = other_bge._cache_available;
_writing_to_cache.store(other_bge._writing_to_cache.load());;
}
#endif
//----------------------------------------------------------------------
// given a quantity in a vector (e.g. pt_over_area) and knowledge
// about the number of empty jets, calculate the median and
// stand_dev_if_gaussian (roughly from the 16th percentile)
//
// If do_fj2_calculation is set to true then this performs FastJet
// 2.X estimation of the standard deviation, which has a spurious
// offset in the limit of a small number of jets.
void BackgroundEstimatorBase::_median_and_stddev(const vector<double> & quantity_vector,
double n_empty_jets,
double & median,
double & stand_dev_if_gaussian,
bool do_fj2_calculation) const {
// this check is redundant (the code below behaves sensibly even
// with a zero size), but serves as a reminder of what happens if
// the quantity vector is zero-sized
if (quantity_vector.size() == 0) {
median = 0;
stand_dev_if_gaussian = 0;
return;
}
vector<double> sorted_quantity_vector = quantity_vector;
sort(sorted_quantity_vector.begin(), sorted_quantity_vector.end());
// empty area can sometimes be negative; with small ranges this can
// become pathological, so warn the user
int n_jets_used = sorted_quantity_vector.size();
if (n_empty_jets < -n_jets_used/4.0)
_warnings_empty_area.warn("BackgroundEstimatorBase::_median_and_stddev(...): the estimated empty area is suspiciously large and negative and may lead to an over-estimation of rho. This may be due to (i) a rare statistical fluctuation or (ii) too small a range used to estimate the background properties.");
// now get the median & error, accounting for empty jets;
// define the fractions of distribution at median, median-1sigma
double posn[2] = {0.5, (1.0-0.6827)/2.0};
double res[2];
for (int i = 0; i < 2; i++) {
res[i] = _percentile(sorted_quantity_vector, posn[i], n_empty_jets,
do_fj2_calculation);
}
median = res[0];
stand_dev_if_gaussian = res[0] - res[1];
}
//----------------------------------------------------------------------
// computes a percentile of a given _sorted_ vector of quantities
// - sorted_quantities the (sorted) vector contains the data sample
// - percentile the percentile (defined between 0 and 1) to compute
// - nempty an additional number of 0's
// (considered at the beginning of
// the quantity vector)
// - do_fj2_calculation carry out the calculation as it
// was done in fj2 (suffers from "edge effects")
double BackgroundEstimatorBase::_percentile(const vector<double> & sorted_quantities,
const double percentile,
const double nempty,
const bool do_fj2_calculation
) const {
assert(percentile >= 0.0 && percentile <= 1.0);
int quantities_size = sorted_quantities.size();
if (quantities_size == 0) return 0;
double total_njets = quantities_size + nempty;
double percentile_pos;
if (do_fj2_calculation) {
percentile_pos = (total_njets-1)*percentile - nempty;
} else {
percentile_pos = (total_njets)*percentile - nempty - 0.5;
}
double result;
if (percentile_pos >= 0 && quantities_size > 1) {
int int_percentile_pos = int(percentile_pos);
// avoid potential overflow issues
if (int_percentile_pos+1 > quantities_size-1){
int_percentile_pos = quantities_size-2;
percentile_pos = quantities_size-1;
}
result =
sorted_quantities[int_percentile_pos] * (int_percentile_pos+1-percentile_pos)
+ sorted_quantities[int_percentile_pos+1] * (percentile_pos - int_percentile_pos);
} else if (percentile_pos > -0.5 && quantities_size >= 1
&& !do_fj2_calculation) {
// in the LHS of this "bin", just keep a constant value (we could have
// interpolated to zero, but this might misbehave in cases where all jets
// are active, because it would go to zero too fast)
result = sorted_quantities[0];
} else {
result = 0.0;
}
return result;
}
void BackgroundEstimatorBase::_lock_if_needed() const{
#ifdef FASTJET_HAVE_THREAD_SAFETY
bool expected;
// the following waits until the cache_writing status is "false" and sets it to "true"
do {
expected = false;
} while (!_writing_to_cache.compare_exchange_strong(expected, true,
memory_order_seq_cst,
memory_order_relaxed));
#endif // FASTJET_HAVE_THREAD_SAFETY
}
void BackgroundEstimatorBase::_unlock_if_needed() const{
#ifdef FASTJET_HAVE_THREAD_SAFETY
// release the "write-in-progress" lock
_writing_to_cache = false;
#endif
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::BackgroundEstimatorBase::_median_and_stddev
void _median_and_stddev(const std::vector< double > &quantity_vector, double n_empty_jets, double &median, double &stand_dev_if_gaussian, bool do_fj2_calculation=false) const
given a quantity in a vector (e.g.
Definition BackgroundEstimatorBase.cc:57
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