FastJet 3.4.3
Loading...
Searching...
No Matches
TilingExtent.cc
1//FJSTARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2024, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5//
6//----------------------------------------------------------------------
7// This file is part of FastJet.
8//
9// FastJet is free software; you can redistribute it and/or modify
10// it under the terms of the GNU General Public License as published by
11// the Free Software Foundation; either version 2 of the License, or
12// (at your option) any later version.
13//
14// The algorithms that underlie FastJet have required considerable
15// development. They are described in the original FastJet paper,
16// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17// FastJet as part of work towards a scientific publication, please
18// quote the version you use and include a citation to the manual and
19// optionally also to hep-ph/0512210.
20//
21// FastJet is distributed in the hope that it will be useful,
22// but WITHOUT ANY WARRANTY; without even the implied warranty of
23// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24// GNU General Public License for more details.
25//
26// You should have received a copy of the GNU General Public License
27// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28//----------------------------------------------------------------------
29//FJENDHEADER
30
31#include <iomanip>
32#include <limits>
33#include <cmath>
34#include "fastjet/internal/TilingExtent.hh"
35using namespace std;
36
37
38FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39
40TilingExtent::TilingExtent(ClusterSequence & cs) {
41 _determine_rapidity_extent(cs.jets());
42}
43
44TilingExtent::TilingExtent(const vector<PseudoJet> &particles) {
45 _determine_rapidity_extent(particles);
46}
47
48void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
49 // have a binning of rapidity that goes from -nrap to nrap
50 // in bins of size 1; the left and right-most bins include
51 // include overflows from smaller/larger rapidities
52 int nrap = 20;
53 int nbins = 2*nrap;
54 vector<double> counts(nbins, 0);
55
56 // get the minimum and maximum rapidities and at the same time bin
57 // the multiplicities as a function of rapidity to help decide how
58 // far out it's worth going
59 _minrap = numeric_limits<double>::max();
60 _maxrap = -numeric_limits<double>::max();
61 int ibin;
62 for (unsigned i = 0; i < particles.size(); i++) {
63 // ignore particles with infinite rapidity
64 if (particles[i].E() == abs(particles[i].pz())) continue;
65 double rap = particles[i].rap();
66 if (rap < _minrap) _minrap = rap;
67 if (rap > _maxrap) _maxrap = rap;
68 // now bin the rapidity to decide how far to go with the tiling.
69 // Remember the bins go from ibin=0 (rap=-infinity..-19)
70 // to ibin = nbins-1 (rap=19..infinity for nrap=20)
71 ibin = int(rap+nrap);
72 if (ibin < 0) ibin = 0;
73 if (ibin >= nbins) ibin = nbins - 1;
74 counts[ibin]++;
75 }
76
77 // now figure out the particle count in the busiest bin
78 double max_in_bin = 0;
79 for (ibin = 0; ibin < nbins; ibin++) {
80 if (max_in_bin < counts[ibin]) max_in_bin = counts[ibin];
81 }
82
83 // and find _minrap, _maxrap such that edge bin never contains more
84 // than some fraction of busiest, and at least a few particles; first do
85 // it from left. NB: the thresholds chosen here are largely
86 // guesstimates as to what might work.
87 //
88 // 2014-07-17: in some tests at high multiplicity (100k) and particles going up to
89 // about 7.3, anti-kt R=0.4, we found that 0.25 gave 20% better run times
90 // than the original value of 0.5.
91 const double allowed_max_fraction = 0.25;
92 // the edge bins should also contain at least min_multiplicity particles
93 const double min_multiplicity = 4;
94 // now calculate how much we can accumulate into an edge bin
95 double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
96 // make sure we don't require more particles in a bin than max_in_bin
97 if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
98
99 // start scan over rapidity bins from the left, to find out minimum rapidity of tiling
100 double cumul_lo = 0;
101 _cumul2 = 0;
102 for (ibin = 0; ibin < nbins; ibin++) {
103 cumul_lo += counts[ibin];
104 if (cumul_lo >= allowed_max_cumul) {
105 double y = ibin-nrap;
106 if (y > _minrap) _minrap = y;
107 break;
108 }
109 }
110 assert(ibin != nbins); // internal consistency check that you found a bin
111 _cumul2 += cumul_lo*cumul_lo;
112
113 // ibin_lo is the index of the leftmost bin that should be considered
114 int ibin_lo = ibin;
115
116 // then do it from right, to find out maximum rapidity of tiling
117 double cumul_hi = 0;
118 for (ibin = nbins-1; ibin >= 0; ibin--) {
119 cumul_hi += counts[ibin];
120 if (cumul_hi >= allowed_max_cumul) {
121 double y = ibin-nrap+1; // +1 here is the rapidity bin width
122 if (y < _maxrap) _maxrap = y;
123 break;
124 }
125 }
126 assert(ibin >= 0); // internal consistency check that you found a bin
127
128 // ibin_hi is the index of the rightmost bin that should be considered
129 int ibin_hi = ibin;
130
131 // consistency check
132 assert(ibin_hi >= ibin_lo);
133
134 // now work out cumul2
135 if (ibin_hi == ibin_lo) {
136 // if there is a single bin (potentially including overflows
137 // from both sides), cumul2 is the square of the total contents
138 // of that bin, which we obtain from cumul_lo and cumul_hi minus
139 // the double counting of part that is contained in both
140 // (putting double
141 _cumul2 = pow(double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
142 } else {
143 // otherwise we have a straightforward sum of squares of bin
144 // contents
145 _cumul2 += cumul_hi*cumul_hi;
146
147 // now get the rest of the squared bin contents
148 for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
149 _cumul2 += counts[ibin]*counts[ibin];
150 }
151 }
152
153}
154
155
156FASTJET_END_NAMESPACE
deals with clustering
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...