1 | //STARTHEADER |
---|
2 | // $Id: GridJetPlugin.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
3 | // |
---|
4 | // Copyright (c) 2011, Matteo Cacciari, Gavin 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 and are described in hep-ph/0512210. If you use |
---|
16 | // FastJet as part of work towards a scientific publication, please |
---|
17 | // include a citation to the FastJet paper. |
---|
18 | // |
---|
19 | // FastJet is distributed in the hope that it will be useful, |
---|
20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
22 | // GNU General Public License for more details. |
---|
23 | // |
---|
24 | // You should have received a copy of the GNU General Public License |
---|
25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
26 | //---------------------------------------------------------------------- |
---|
27 | //ENDHEADER |
---|
28 | |
---|
29 | // fastjet stuff |
---|
30 | #include "fastjet/ClusterSequence.hh" |
---|
31 | #include "fastjet/GridJetPlugin.hh" |
---|
32 | |
---|
33 | // other stuff |
---|
34 | #include <vector> |
---|
35 | #include <sstream> |
---|
36 | |
---|
37 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
38 | |
---|
39 | using namespace std; |
---|
40 | |
---|
41 | //---------------------------------------------------------------------- |
---|
42 | GridJetPlugin::GridJetPlugin (double ymax, |
---|
43 | double requested_grid_spacing, |
---|
44 | const JetDefinition & post_jet_def) : |
---|
45 | _ymin(-ymax), _ymax(ymax), |
---|
46 | _requested_grid_spacing(requested_grid_spacing) , |
---|
47 | _post_jet_def(post_jet_def) |
---|
48 | { |
---|
49 | setup_grid(); |
---|
50 | } |
---|
51 | |
---|
52 | void GridJetPlugin::setup_grid() { |
---|
53 | // since we've exchanged the arguments of the constructor, |
---|
54 | // there's a danger of calls with exchanged ymax,spacing arguments -- |
---|
55 | // the following check should catch most such situations. |
---|
56 | assert(_ymax>0 && _ymax - _ymin >= _requested_grid_spacing); |
---|
57 | |
---|
58 | double ny_double = (_ymax-_ymin) / _requested_grid_spacing; |
---|
59 | _ny = int(ny_double+0.49999); |
---|
60 | _dy = (_ymax-_ymin) / _ny; |
---|
61 | |
---|
62 | _nphi = int (twopi / _requested_grid_spacing + 0.5); |
---|
63 | _dphi = twopi / _nphi; |
---|
64 | |
---|
65 | // some sanity checking (could throw a fastjet::Error) |
---|
66 | assert(_ny >= 1 && _nphi >= 1); |
---|
67 | |
---|
68 | _ntotal = _nphi * _ny; |
---|
69 | } |
---|
70 | |
---|
71 | |
---|
72 | //---------------------------------------------------------------------- |
---|
73 | string GridJetPlugin::description () const { |
---|
74 | ostringstream desc; |
---|
75 | desc << "GridJetPlugin plugin with ymax = " << _ymax << ", dy = " << _dy << ", dphi = " << _dphi << " (requested grid spacing was " << _requested_grid_spacing << ")"; |
---|
76 | if (_post_jet_def.jet_algorithm() != undefined_jet_algorithm) { |
---|
77 | desc << ", followed by " << _post_jet_def.description(); |
---|
78 | } |
---|
79 | return desc.str(); |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | //---------------------------------------------------------------------- |
---|
84 | double GridJetPlugin::R() const {return sqrt(_dy*_dphi/pi);} |
---|
85 | |
---|
86 | |
---|
87 | //---------------------------------------------------------------------- |
---|
88 | int GridJetPlugin::igrid(const PseudoJet & p) const { |
---|
89 | // directly taking int does not work for values between -1 and 0 |
---|
90 | // so use floor instead |
---|
91 | // double iy_double = (p.rap() - _ymin) / _dy; |
---|
92 | // if (iy_double < 0.0) return -1; |
---|
93 | // int iy = int(iy_double); |
---|
94 | // if (iy >= _ny) return -1; |
---|
95 | |
---|
96 | // writing it as below gives a huge speed gain (factor two!). Even |
---|
97 | // though answers are identical and the routine here is not the |
---|
98 | // speed-critical step. It's not at all clear why. |
---|
99 | int iy = int(floor( (p.rap() - _ymin) / _dy )); |
---|
100 | if (iy < 0 || iy >= _ny) return -1; |
---|
101 | |
---|
102 | int iphi = int( p.phi()/_dphi ); |
---|
103 | assert(iphi >= 0 && iphi <= _nphi); |
---|
104 | if (iphi == _nphi) iphi = 0; // just in case of rounding errors |
---|
105 | |
---|
106 | int igrid_res = iy*_nphi + iphi; |
---|
107 | assert (igrid_res >= 0 && igrid_res < _ny*_nphi); |
---|
108 | return igrid_res; |
---|
109 | } |
---|
110 | |
---|
111 | |
---|
112 | //---------------------------------------------------------------------- |
---|
113 | void GridJetPlugin::run_clustering(ClusterSequence & cs) const { |
---|
114 | |
---|
115 | // we will create a grid; |
---|
116 | // * -1 will indicate there is no jet here currently |
---|
117 | // * a number >= 0 will mean that particle indicated by the index |
---|
118 | // is currently the jet on the grid |
---|
119 | vector<int> grid(_ntotal, -1); |
---|
120 | |
---|
121 | int nparticles = cs.jets().size(); |
---|
122 | double dij_or_diB = 1.0; |
---|
123 | |
---|
124 | int ngrid_active = 0; |
---|
125 | |
---|
126 | // combine particles with whatever is in the grid |
---|
127 | for (int i = 0; i < nparticles; i++) { |
---|
128 | int igrd = igrid(cs.jets()[i]); |
---|
129 | //cout << i << " " << cs.jets()[i].rap() << " " << cs.jets()[i].phi() |
---|
130 | // << " " << igrd << " " << grid.size() << " " << _ntotal << endl; |
---|
131 | if (igrd < 0) continue; |
---|
132 | assert(igrd <= _ntotal); |
---|
133 | if (grid[igrd] == -1) { |
---|
134 | grid[igrd] = i; // jet index of initial particle i is i |
---|
135 | ngrid_active++; |
---|
136 | } else { |
---|
137 | int k; |
---|
138 | cs.plugin_record_ij_recombination(grid[igrd], i, dij_or_diB, k); |
---|
139 | grid[igrd] = k; // grid takes jet index of new particle |
---|
140 | //cout << " res: " << cs.jets()[k].rap() << " " << cs.jets()[k].phi() << endl; |
---|
141 | } |
---|
142 | } |
---|
143 | |
---|
144 | if (_post_jet_def.jet_algorithm() == undefined_jet_algorithm) { |
---|
145 | // make the final jets via iB recombinations |
---|
146 | for (unsigned igrd = 0; igrd < grid.size(); igrd++) { |
---|
147 | if (grid[igrd] != -1) cs.plugin_record_iB_recombination(grid[igrd], |
---|
148 | dij_or_diB); |
---|
149 | } |
---|
150 | } else { |
---|
151 | // otherwise post-cluster the grid elements with a normal jet algorithm |
---|
152 | vector<PseudoJet> inputs; |
---|
153 | vector<int> cs_indices; |
---|
154 | inputs.reserve(ngrid_active); |
---|
155 | cs_indices.reserve(2*ngrid_active); |
---|
156 | for (unsigned igrd = 0; igrd < grid.size(); igrd++) { |
---|
157 | if (grid[igrd] != -1) { |
---|
158 | inputs.push_back(cs.jets()[grid[igrd]]); |
---|
159 | cs_indices.push_back(grid[igrd]); |
---|
160 | } |
---|
161 | } |
---|
162 | ClusterSequence post_cs(inputs, _post_jet_def); |
---|
163 | const vector<ClusterSequence::history_element> & post_history = post_cs.history(); |
---|
164 | const vector<PseudoJet> & post_jets = post_cs.jets(); |
---|
165 | for (unsigned ihist = ngrid_active; ihist < post_history.size(); ihist++) { |
---|
166 | const ClusterSequence::history_element & hist = post_history[ihist]; |
---|
167 | int post_ij1 = post_history[hist.parent1].jetp_index; |
---|
168 | int ij1 = cs_indices[post_ij1]; |
---|
169 | if (hist.parent2 >= 0) { |
---|
170 | int post_ij2 = post_history[hist.parent2].jetp_index; |
---|
171 | int ij2 = cs_indices[post_ij2]; |
---|
172 | int k; |
---|
173 | cs.plugin_record_ij_recombination(ij1, ij2, hist.dij, post_jets[hist.jetp_index], k); |
---|
174 | assert(int(cs_indices.size()) == hist.jetp_index); |
---|
175 | cs_indices.push_back(k); |
---|
176 | } else { |
---|
177 | cs.plugin_record_iB_recombination(ij1, hist.dij); |
---|
178 | } |
---|
179 | } |
---|
180 | |
---|
181 | } |
---|
182 | } |
---|
183 | |
---|
184 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh |
---|