Fork me on GitHub

source: git/external/fastjet/plugins/GridJet/GridJetPlugin.cc@ 94083e7

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 94083e7 was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 7.4 KB
Line 
1//FJSTARTHEADER
2// $Id: GridJetPlugin.cc 2268 2011-06-20 15:12:26Z salam $
3//
4// Copyright (c) 2005-2018, 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
39FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40
41using namespace std;
42
43//----------------------------------------------------------------------
44GridJetPlugin::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
60GridJetPlugin::GridJetPlugin (const RectangularGrid & grid,
61 const JetDefinition & post_jet_def) :
62 RectangularGrid(grid), _post_jet_def(post_jet_def) {
63 if (!RectangularGrid::is_initialised())
64 throw Error("attempt to construct GridJetPlugin with uninitialised RectangularGrid");
65}
66#endif // FASTJET_GRIDJET_USEFJGRID
67
68#ifndef FASTJET_GRIDJET_USEFJGRID
69void 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//----------------------------------------------------------------------
89int 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//----------------------------------------------------------------------
115string 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
121 desc << RectangularGrid::description();
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//----------------------------------------------------------------------
131double GridJetPlugin::R() const {return sqrt(drap()*dphi()/pi);}
132
133
134//----------------------------------------------------------------------
135void GridJetPlugin::run_clustering(ClusterSequence & cs) const {
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
206FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.