Loading...
Searching...
No Matches
TilingExtent.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 <iomanip>
#include <limits>
#include <cmath>
#include "fastjet/internal/TilingExtent.hh"
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
TilingExtent::TilingExtent(ClusterSequence & cs) {
_determine_rapidity_extent(cs.jets());
}
TilingExtent::TilingExtent(const vector<PseudoJet> &particles) {
_determine_rapidity_extent(particles);
}
void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
// have a binning of rapidity that goes from -nrap to nrap
// in bins of size 1; the left and right-most bins include
// include overflows from smaller/larger rapidities
int nrap = 20;
int nbins = 2*nrap;
vector<double> counts(nbins, 0);
// get the minimum and maximum rapidities and at the same time bin
// the multiplicities as a function of rapidity to help decide how
// far out it's worth going
_minrap = numeric_limits<double>::max();
_maxrap = -numeric_limits<double>::max();
int ibin;
for (unsigned i = 0; i < particles.size(); i++) {
// ignore particles with infinite rapidity
if (particles[i].E() == abs(particles[i].pz())) continue;
double rap = particles[i].rap();
if (rap < _minrap) _minrap = rap;
if (rap > _maxrap) _maxrap = rap;
// now bin the rapidity to decide how far to go with the tiling.
// Remember the bins go from ibin=0 (rap=-infinity..-19)
// to ibin = nbins-1 (rap=19..infinity for nrap=20)
ibin = int(rap+nrap);
if (ibin < 0) ibin = 0;
if (ibin >= nbins) ibin = nbins - 1;
counts[ibin]++;
}
// now figure out the particle count in the busiest bin
double max_in_bin = 0;
for (ibin = 0; ibin < nbins; ibin++) {
if (max_in_bin < counts[ibin]) max_in_bin = counts[ibin];
}
// and find _minrap, _maxrap such that edge bin never contains more
// than some fraction of busiest, and at least a few particles; first do
// it from left. NB: the thresholds chosen here are largely
// guesstimates as to what might work.
//
// 2014-07-17: in some tests at high multiplicity (100k) and particles going up to
// about 7.3, anti-kt R=0.4, we found that 0.25 gave 20% better run times
// than the original value of 0.5.
const double allowed_max_fraction = 0.25;
// the edge bins should also contain at least min_multiplicity particles
const double min_multiplicity = 4;
// now calculate how much we can accumulate into an edge bin
double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
// make sure we don't require more particles in a bin than max_in_bin
if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
// start scan over rapidity bins from the left, to find out minimum rapidity of tiling
double cumul_lo = 0;
_cumul2 = 0;
for (ibin = 0; ibin < nbins; ibin++) {
cumul_lo += counts[ibin];
if (cumul_lo >= allowed_max_cumul) {
double y = ibin-nrap;
if (y > _minrap) _minrap = y;
break;
}
}
assert(ibin != nbins); // internal consistency check that you found a bin
_cumul2 += cumul_lo*cumul_lo;
// ibin_lo is the index of the leftmost bin that should be considered
int ibin_lo = ibin;
// then do it from right, to find out maximum rapidity of tiling
double cumul_hi = 0;
for (ibin = nbins-1; ibin >= 0; ibin--) {
cumul_hi += counts[ibin];
if (cumul_hi >= allowed_max_cumul) {
double y = ibin-nrap+1; // +1 here is the rapidity bin width
if (y < _maxrap) _maxrap = y;
break;
}
}
assert(ibin >= 0); // internal consistency check that you found a bin
// ibin_hi is the index of the rightmost bin that should be considered
int ibin_hi = ibin;
// consistency check
assert(ibin_hi >= ibin_lo);
// now work out cumul2
if (ibin_hi == ibin_lo) {
// if there is a single bin (potentially including overflows
// from both sides), cumul2 is the square of the total contents
// of that bin, which we obtain from cumul_lo and cumul_hi minus
// the double counting of part that is contained in both
// (putting double
_cumul2 = pow(double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
} else {
// otherwise we have a straightforward sum of squares of bin
// contents
_cumul2 += cumul_hi*cumul_hi;
// now get the rest of the squared bin contents
for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
_cumul2 += counts[ibin]*counts[ibin];
}
}
}
FASTJET_END_NAMESPACE
fastjet::ClusterSequence
deals with clustering
Definition ClusterSequence.hh:63
fastjet::ClusterSequence::jets
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Definition ClusterSequence.hh:1003
fastjet::TilingExtent::TilingExtent
TilingExtent(ClusterSequence &cs)
constructor that takes a ClusterSequence in a state where the initial particles have been set up,...
Definition TilingExtent.cc:40