FastJet  3.3.0
TilingExtent.cc
1 //FJSTARTHEADER
2 // $Id: TilingExtent.cc 4034 2016-03-02 00:20:27Z soyez $
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 TilingExtent::TilingExtent(const vector<PseudoJet> &particles) {
45  _determine_rapidity_extent(particles);
46 }
47 
48 void 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 
156 FASTJET_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...