FastJet  3.3.0
GridJetPlugin.cc
1 //FJSTARTHEADER
2 // $Id: GridJetPlugin.cc 2268 2011-06-20 15:12:26Z 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 // fastjet stuff
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/GridJetPlugin.hh"
34 
35 // other stuff
36 #include <vector>
37 #include <sstream>
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 using namespace std;
42 
43 //----------------------------------------------------------------------
44 GridJetPlugin::GridJetPlugin (double ymax,
45  double requested_grid_spacing,
46  const JetDefinition & post_jet_def) :
47 #ifdef FASTJET_GRIDJET_USEFJGRID
48  RectangularGrid(ymax, requested_grid_spacing), _post_jet_def(post_jet_def) {
49 }
50 #else
51  _ymin(-ymax), _ymax(ymax),
52  _requested_grid_spacing(requested_grid_spacing) ,
53  _post_jet_def(post_jet_def)
54 {
55  setup_grid();
56 }
57 #endif
58 
59 #ifdef FASTJET_GRIDJET_USEFJGRID
61  const JetDefinition & post_jet_def) :
62  RectangularGrid(grid), _post_jet_def(post_jet_def) {
64  throw Error("attempt to construct GridJetPlugin with uninitialised RectangularGrid");
65 }
66 #endif // FASTJET_GRIDJET_USEFJGRID
67 
68 #ifndef FASTJET_GRIDJET_USEFJGRID
69 void GridJetPlugin::setup_grid() {
70  // since we've exchanged the arguments of the constructor,
71  // there's a danger of calls with exchanged ymax,spacing arguments --
72  // the following check should catch most such situations.
73  assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing);
74 
75  double ny_double = (_ymax-_ymin) / _requested_grid_spacing;
76  _ny = int(ny_double+0.49999);
77  _dy = (_ymax-_ymin) / _ny;
78 
79  _nphi = int (twopi / _requested_grid_spacing + 0.5);
80  _dphi = twopi / _nphi;
81 
82  // some sanity checking (could throw a fastjet::Error)
83  assert(_ny >= 1 && _nphi >= 1);
84 
85  _ntotal = _nphi * _ny;
86 }
87 
88 //----------------------------------------------------------------------
89 int GridJetPlugin::tile_index(const PseudoJet & p) const {
90  // directly taking int does not work for values between -1 and 0
91  // so use floor instead
92  // double iy_double = (p.rap() - _ymin) / _dy;
93  // if (iy_double < 0.0) return -1;
94  // int iy = int(iy_double);
95  // if (iy >= _ny) return -1;
96 
97  // writing it as below gives a huge speed gain (factor two!). Even
98  // though answers are identical and the routine here is not the
99  // speed-critical step. It's not at all clear why.
100  int iy = int(floor( (p.rap() - _ymin) / _dy ));
101  if (iy < 0 || iy >= _ny) return -1;
102 
103  int iphi = int( p.phi()/_dphi );
104  assert(iphi >= 0 && iphi <= _nphi);
105  if (iphi == _nphi) iphi = 0; // just in case of rounding errors
106 
107  int igrid_res = iy*_nphi + iphi;
108  assert (igrid_res >= 0 && igrid_res < _ny*_nphi);
109  return igrid_res;
110 }
111 #endif // not FASTJET_GRIDJET_USEFJGRID
112 
113 
114 //----------------------------------------------------------------------
115 string GridJetPlugin::description () const {
116  ostringstream desc;
117  desc << "GridJetPlugin plugin with ";
118 #ifndef FASTJET_GRIDJET_USEFJGRID
119  desc << "ymax = " << _ymax << ", dy = " << _dy << ", dphi = " << _dphi << " (requested grid spacing was " << _requested_grid_spacing << ")";
120 #else
122 #endif
123  if (_post_jet_def.jet_algorithm() != undefined_jet_algorithm) {
124  desc << ", followed by " << _post_jet_def.description();
125  }
126  return desc.str();
127 }
128 
129 
130 //----------------------------------------------------------------------
131 double GridJetPlugin::R() const {return sqrt(drap()*dphi()/pi);}
132 
133 
134 //----------------------------------------------------------------------
136 
137  // we will create a grid;
138  // * -1 will indicate there is no jet here currently
139  // * a number >= 0 will mean that particle indicated by the index
140  // is currently the jet on the grid
141  vector<int> grid(n_tiles(), -1);
142 
143  int nparticles = cs.jets().size();
144  double dij_or_diB = 1.0;
145 
146  int ngrid_active = 0;
147 
148  // combine particles with whatever is in the grid
149  for (int i = 0; i < nparticles; i++) {
150  int igrd = tile_index(cs.jets()[i]);
151  //cout << i << " " << cs.jets()[i].rap() << " " << cs.jets()[i].phi()
152  // << " " << igrd << " " << grid.size() << " " << _ntotal << endl;
153  if (igrd < 0) continue;
154  assert(igrd <= n_tiles());
155  if (grid[igrd] == -1) {
156  grid[igrd] = i; // jet index of initial particle i is i
157  ngrid_active++;
158  } else {
159  int k;
160  cs.plugin_record_ij_recombination(grid[igrd], i, dij_or_diB, k);
161  grid[igrd] = k; // grid takes jet index of new particle
162  //cout << " res: " << cs.jets()[k].rap() << " " << cs.jets()[k].phi() << endl;
163  }
164  }
165 
166  if (_post_jet_def.jet_algorithm() == undefined_jet_algorithm) {
167  // make the final jets via iB recombinations
168  for (unsigned igrd = 0; igrd < grid.size(); igrd++) {
169  if (grid[igrd] != -1 && tile_is_good(igrd))
170  cs.plugin_record_iB_recombination(grid[igrd], dij_or_diB);
171  }
172  } else {
173  // otherwise post-cluster the grid elements with a normal jet algorithm
174  vector<PseudoJet> inputs;
175  vector<int> cs_indices;
176  inputs.reserve(ngrid_active);
177  cs_indices.reserve(2*ngrid_active);
178  for (unsigned igrd = 0; igrd < grid.size(); igrd++) {
179  if (grid[igrd] != -1) {
180  inputs.push_back(cs.jets()[grid[igrd]]);
181  cs_indices.push_back(grid[igrd]);
182  }
183  }
184  ClusterSequence post_cs(inputs, _post_jet_def);
185  const vector<ClusterSequence::history_element> & post_history = post_cs.history();
186  const vector<PseudoJet> & post_jets = post_cs.jets();
187  for (unsigned ihist = ngrid_active; ihist < post_history.size(); ihist++) {
188  const ClusterSequence::history_element & hist = post_history[ihist];
189  int post_ij1 = post_history[hist.parent1].jetp_index;
190  int ij1 = cs_indices[post_ij1];
191  if (hist.parent2 >= 0) {
192  int post_ij2 = post_history[hist.parent2].jetp_index;
193  int ij2 = cs_indices[post_ij2];
194  int k;
195  cs.plugin_record_ij_recombination(ij1, ij2, hist.dij, post_jets[hist.jetp_index], k);
196  assert(int(cs_indices.size()) == hist.jetp_index);
197  cs_indices.push_back(k);
198  } else {
199  cs.plugin_record_iB_recombination(ij1, hist.dij);
200  }
201  }
202 
203  }
204 }
205 
206 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
GridJetPlugin(double ymax, double requested_grid_spacing, const JetDefinition &post_jet_def=JetDefinition())
Basic constructor for the GridJetPlugin Plugin class.
JetAlgorithm jet_algorithm() const
return information about the definition...
deals with clustering
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined ...
virtual int n_tiles() const override
returns the total number of tiles in the tiling; valid tile indices run from 0 ...
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
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
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
double dij
index in the _jets vector where we will find the
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam...
Class that holds a generic rectangular tiling.
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
virtual double R() const
This returns the sqrt(dphi*dy/pi) – i.e.
virtual bool tile_is_good(int itile) const override
returns whether a given tile is good
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
virtual std::string description() const
return a textual description of the jet-definition implemented in this plugin
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
a single element in the clustering history
virtual void run_clustering(ClusterSequence &) const
given a ClusterSequence that has been filled up with initial particles, the following function should...
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 ...
double drap() const
returns the spacing of the grid in rapidity
std::string description() const
return a textual description of the current jet definition
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
class that is intended to hold a full definition of the jet clusterer
double dphi() const
returns the spacing of the grid in azimuth