FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
TilingExtent.cc
1 //FJSTARTHEADER
2 // $Id: TilingExtent.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2005-2014, 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"
35 using namespace std;
36 
37 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 TilingExtent::TilingExtent(ClusterSequence & cs) {
41  _determine_rapidity_extent(cs.jets());
42 }
43 
44 void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
45  // have a binning of rapidity that goes from -nrap to nrap
46  // in bins of size 1; the left and right-most bins include
47  // include overflows from smaller/larger rapidities
48  int nrap = 20;
49  int nbins = 2*nrap;
50  vector<double> counts(nbins, 0);
51 
52  // get the minimum and maximum rapidities and at the same time bin
53  // the multiplicities as a function of rapidity to help decide how
54  // far out it's worth going
55  _minrap = numeric_limits<double>::max();
56  _maxrap = -numeric_limits<double>::max();
57  int ibin;
58  for (unsigned i = 0; i < particles.size(); i++) {
59  // ignore particles with infinite rapidity
60  if (particles[i].E() == abs(particles[i].pz())) continue;
61  double rap = particles[i].rap();
62  if (rap < _minrap) _minrap = rap;
63  if (rap > _maxrap) _maxrap = rap;
64  // now bin the rapidity to decide how far to go with the tiling.
65  // Remember the bins go from ibin=0 (rap=-infinity..-19)
66  // to ibin = nbins-1 (rap=19..infinity for nrap=20)
67  ibin = int(rap+nrap);
68  if (ibin < 0) ibin = 0;
69  if (ibin >= nbins) ibin = nbins - 1;
70  counts[ibin]++;
71  }
72 
73  // now figure out the particle count in the busiest bin
74  double max_in_bin = 0;
75  for (ibin = 0; ibin < nbins; ibin++) {
76  if (max_in_bin < counts[ibin]) max_in_bin = counts[ibin];
77  }
78 
79  // and find _minrap, _maxrap such that edge bin never contains more
80  // than some fraction of busiest, and at least a few particles; first do
81  // it from left. NB: the thresholds chosen here are largely
82  // guesstimates as to what might work.
83  //
84  // 2014-07-17: in some tests at high multiplicity (100k) and particles going up to
85  // about 7.3, anti-kt R=0.4, we found that 0.25 gave 20% better run times
86  // than the original value of 0.5.
87  const double allowed_max_fraction = 0.25;
88  // the edge bins should also contain at least min_multiplicity particles
89  const double min_multiplicity = 4;
90  // now calculate how much we can accumulate into an edge bin
91  double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
92  // make sure we don't require more particles in a bin than max_in_bin
93  if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
94 
95  // start scan over rapidity bins from the left, to find out minimum rapidity of tiling
96  double cumul_lo = 0;
97  _cumul2 = 0;
98  for (ibin = 0; ibin < nbins; ibin++) {
99  cumul_lo += counts[ibin];
100  if (cumul_lo >= allowed_max_cumul) {
101  double y = ibin-nrap;
102  if (y > _minrap) _minrap = y;
103  break;
104  }
105  }
106  assert(ibin != nbins); // internal consistency check that you found a bin
107  _cumul2 += cumul_lo*cumul_lo;
108 
109  // ibin_lo is the index of the leftmost bin that should be considered
110  int ibin_lo = ibin;
111 
112  // then do it from right, to find out maximum rapidity of tiling
113  double cumul_hi = 0;
114  for (ibin = nbins-1; ibin >= 0; ibin--) {
115  cumul_hi += counts[ibin];
116  if (cumul_hi >= allowed_max_cumul) {
117  double y = ibin-nrap+1; // +1 here is the rapidity bin width
118  if (y < _maxrap) _maxrap = y;
119  break;
120  }
121  }
122  assert(ibin >= 0); // internal consistency check that you found a bin
123 
124  // ibin_hi is the index of the rightmost bin that should be considered
125  int ibin_hi = ibin;
126 
127  // consistency check
128  assert(ibin_hi >= ibin_lo);
129 
130  // now work out cumul2
131  if (ibin_hi == ibin_lo) {
132  // if there is a single bin (potentially including overflows
133  // from both sides), cumul2 is the square of the total contents
134  // of that bin, which we obtain from cumul_lo and cumul_hi minus
135  // the double counting of part that is contained in both
136  // (putting double
137  _cumul2 = pow(double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
138  } else {
139  // otherwise we have a straightforward sum of squares of bin
140  // contents
141  _cumul2 += cumul_hi*cumul_hi;
142 
143  // now get the rest of the squared bin contents
144  for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
145  _cumul2 += counts[ibin]*counts[ibin];
146  }
147  }
148 
149 }
150 
151 
152 FASTJET_END_NAMESPACE
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
deals with clustering