FastJet  3.3.3
RectangularGrid.cc
1 //FJSTARTHEADER
2 // $Id: RectangularGrid.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, 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 
32 
33 #include "fastjet/RectangularGrid.hh"
34 #include <sstream>
35 using namespace std;
36 
37 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38 
39  /// dummy ctor (will give an unusable grid)
40 RectangularGrid::RectangularGrid()
41  : _ymax(-1.0), _ymin(1.0), _requested_drap(-1.0), _requested_dphi(-1.0) {
42  // put in nonsensical values for the other variables too, to keep coverity happy
43  _ntotal = -1;
44  _ngood = -1;
45  _dy = 0.0;
46  _dphi = 0.0;
47  _cell_area = 0.0;
48  _inverse_dy = 0;
49  _inverse_dphi = 0;
50  _ny = 0;
51  _nphi = 0;
52 }
53 
54 
55 int RectangularGrid::tile_index(const PseudoJet & p) const {
56  // the code below has seem some degree of optimization: don't change
57  // it without testing the speed again
58 
59  // new version as of 2014-08-04
60  double y_minus_ymin = p.rap() - _ymin;
61  if (y_minus_ymin < 0) return -1;
62  int iy = int(y_minus_ymin * _inverse_dy); // guaranteed positive, so int is safe
63  if (iy >= _ny) return -1;
64 
65  // old version: gives a SoftKiller that's about 10% slower on Gavin's retina mac.
66  // (though having it hard coded inside SoftKiller returns that advantage)
67  // BUT: some comments said that this was a factor of two faster than
68  // something similar to the version above. What is going on?
69  // int iy = int(floor( (p.rap() - _ymin) * _inverse_dy ));
70  // if (iy < 0 || iy >= _ny) return -1;
71 
72  int iphi = int( p.phi() * _inverse_dphi );
73  if (iphi == _nphi) iphi = 0; // just in case of rounding errors
74 
75  return iy*_nphi + iphi;
76 }
77 
78 
79 //----------------------------------------------------------------------
80 void RectangularGrid::_setup_grid() {
81  // initial sanity checks
82  assert(_ymax > _ymin);
83  assert(_requested_drap > 0);
84  assert(_requested_dphi > 0);
85 
86  double ny_double = (_ymax-_ymin) / _requested_drap;
87  _ny = max(int(ny_double+0.5),1);
88  _dy = (_ymax-_ymin) / _ny;
89  _inverse_dy = _ny/(_ymax-_ymin);
90 
91  _nphi = int (twopi / _requested_dphi + 0.5);
92  _dphi = twopi / _nphi;
93  _inverse_dphi = _nphi/twopi;
94 
95  // some sanity checking (could throw a fastjet::Error)
96  assert(_ny >= 1 && _nphi >= 1);
97 
98  _ntotal = _nphi * _ny;
99  //_max_pt.resize(_ntotal);
100  _cell_area = _dy * _dphi;
101 
102  // if we have a selector, establish which tiles are good;
103  // apply the selector to a 4-vector at the tile's centre
104  if (_tile_selector.worker()) {
105  _is_good.resize(n_tiles());
106  _ngood = 0;
107  for (int i = 0; i < n_tiles(); i++) {
108  int iphi = i % _nphi;
109  int irap = i / _nphi;
110  double phi = (iphi + 0.5)*_dphi;
111  double rap = (irap + 0.5)*_dy + _ymin;
112  _is_good[i] = _tile_selector.pass(PtYPhiM(1.0, rap, phi));
113  if (_is_good[i]) _ngood++;
114  }
115  } else {
116  _ngood = n_tiles();
117  }
118 }
119 
120 //----------------------------------------------------------------------
122  if (! is_initialised())
123  return "Uninitialised rectangular grid";
124 
125  ostringstream oss;
126  oss << "rectangular grid with rapidity extent " << _ymin << " < rap < " << _ymax
127  << ", tile size drap x dphi = " << _dy << " x " << _dphi;
128 
129  if (_tile_selector.worker()) {
130  oss << ", good tiles are those that pass selector " << _tile_selector.description();
131  }
132  return oss.str();
133 }
134 
135 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition: PseudoJet.cc:365
virtual int n_tiles() const override
returns the total number of tiles in the tiling; valid tile indices run from 0 ...
virtual std::string description() const override
returns a textual description of the grid
virtual bool is_initialised() const override
returns true if the grid is in a suitably initialised state
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker&#39;s shared pointer
Definition: Selector.hh:277
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:178
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
virtual int tile_index(const PseudoJet &p) const override
returns the index of the tile in which p is located, or -1 if p is outside the tiling region ...
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:237
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67