1 | //STARTHEADER
|
---|
2 | // $Id$
|
---|
3 | //
|
---|
4 | // Copyright (c) 2005-2011, 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 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 | #include "fastjet/ClusterSequence.hh"
|
---|
30 | #include "fastjet/internal/ClosestPair2D.hh"
|
---|
31 | #include<limits>
|
---|
32 | #include<vector>
|
---|
33 | #include<cmath>
|
---|
34 | #include<iostream>
|
---|
35 |
|
---|
36 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
---|
37 |
|
---|
38 | using namespace std;
|
---|
39 |
|
---|
40 | // place for things we don't want outside world to run into
|
---|
41 | namespace Private {
|
---|
42 | /// class for helping us deal with mirror-image particles.
|
---|
43 | class MirrorInfo{
|
---|
44 | public:
|
---|
45 | int orig, mirror;
|
---|
46 | MirrorInfo(int a, int b) : orig(a), mirror(b) {};
|
---|
47 | MirrorInfo() {};
|
---|
48 | };
|
---|
49 |
|
---|
50 | /// if there is a need for a mirror when looking for closest pairs
|
---|
51 | /// up to distance D, then return true and turn the supplied point
|
---|
52 | /// into its mirror copy
|
---|
53 | bool make_mirror(Coord2D & point, double Dlim) {
|
---|
54 | if (point.y < Dlim) {point.y += twopi; return true;}
|
---|
55 | if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
|
---|
56 | return false;
|
---|
57 | }
|
---|
58 |
|
---|
59 | }
|
---|
60 |
|
---|
61 | using namespace Private;
|
---|
62 |
|
---|
63 |
|
---|
64 | //----------------------------------------------------------------------
|
---|
65 | /// clusters only up to a distance Dlim -- does not deal with "inclusive" jets
|
---|
66 | /// -- these are left to some other part of the program
|
---|
67 | void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) {
|
---|
68 |
|
---|
69 | unsigned int n = _initial_n;
|
---|
70 |
|
---|
71 | vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID
|
---|
72 | vector<int> jetIDs(2*n); // jet ID for a given coord ID
|
---|
73 | vector<Coord2D> coords(2*n); // our coordinates (and copies)
|
---|
74 |
|
---|
75 | // particles within a distance Dlim of the phi edges (phi<Dlim ||
|
---|
76 | // phi>2pi-Dli;) will be mirrored. For Dlim>pi, this could lead to
|
---|
77 | // particles copies outside the fixed range in phi which is
|
---|
78 | // [-pi:3pi] (see make_mirror above). Since in that case all
|
---|
79 | // particles get copied anywaym we can just copy particles up to a
|
---|
80 | // distance "pi" from the edges
|
---|
81 | double Dlim4mirror = min(Dlim,pi);
|
---|
82 |
|
---|
83 | // start things off...
|
---|
84 | double minrap = numeric_limits<double>::max();
|
---|
85 | double maxrap = -minrap;
|
---|
86 | int coord_index = -1;
|
---|
87 | int n_active = 0;
|
---|
88 | for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
|
---|
89 |
|
---|
90 | // skip jets that already have children or that have infinite
|
---|
91 | // rapidity
|
---|
92 | if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
|
---|
93 | (_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
|
---|
94 | _jets[jet_i].perp2() == 0.0)
|
---|
95 | ) {continue;}
|
---|
96 |
|
---|
97 | n_active++;
|
---|
98 |
|
---|
99 | coordIDs[jet_i].orig = ++coord_index;
|
---|
100 | coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
|
---|
101 | jetIDs[coord_index] = jet_i;
|
---|
102 | minrap = min(coords[coord_index].x,minrap);
|
---|
103 | maxrap = max(coords[coord_index].x,maxrap);
|
---|
104 |
|
---|
105 | Coord2D mirror_point(coords[coord_index]);
|
---|
106 | if (make_mirror(mirror_point, Dlim4mirror)) {
|
---|
107 | coordIDs[jet_i].mirror = ++coord_index;
|
---|
108 | coords[coord_index] = mirror_point;
|
---|
109 | jetIDs[coord_index] = jet_i;
|
---|
110 | } else {
|
---|
111 | coordIDs[jet_i].mirror = Invalid;
|
---|
112 | }
|
---|
113 | }
|
---|
114 |
|
---|
115 | // set them to their actual size...
|
---|
116 | coords.resize(coord_index+1);
|
---|
117 |
|
---|
118 | // establish limits (with some leeway on rapidity)
|
---|
119 | Coord2D left_edge(minrap-1.0, -3.15); // a security margin below -pi
|
---|
120 | Coord2D right_edge(maxrap+1.0, 9.45); // a security margin above 3*pi
|
---|
121 |
|
---|
122 | //cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl;
|
---|
123 |
|
---|
124 | // now create the closest pair search object
|
---|
125 | ClosestPair2D cp(coords, left_edge, right_edge);
|
---|
126 |
|
---|
127 | // cross check to see what's going on...
|
---|
128 | //cerr << "Tree depths before: ";
|
---|
129 | //cp.print_tree_depths(cerr);
|
---|
130 |
|
---|
131 | // and start iterating...
|
---|
132 | vector<Coord2D> new_points(2);
|
---|
133 | vector<unsigned int> cIDs_to_remove(4);
|
---|
134 | vector<unsigned int> new_cIDs(2);
|
---|
135 |
|
---|
136 | do {
|
---|
137 | // find out which pair we recombine
|
---|
138 | unsigned int cID1, cID2;
|
---|
139 | double distance2;
|
---|
140 | cp.closest_pair(cID1,cID2,distance2);
|
---|
141 |
|
---|
142 | // if the closest distance > Dlim, we exit and go to "inclusive" stage
|
---|
143 | if (distance2 > Dlim*Dlim) {break;}
|
---|
144 |
|
---|
145 | // normalise distance as we like it
|
---|
146 | distance2 *= _invR2;
|
---|
147 |
|
---|
148 | // do the recombination...
|
---|
149 | int jet_i = jetIDs[cID1];
|
---|
150 | int jet_j = jetIDs[cID2];
|
---|
151 | assert (jet_i != jet_j); // to catch issue of recombining with mirror point
|
---|
152 | int newjet_k;
|
---|
153 | _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
|
---|
154 |
|
---|
155 | // don't bother with any further action if only one active particle
|
---|
156 | // is left (also avoid closest-pair error [cannot remove last particle]).
|
---|
157 | if (--n_active == 1) {break;}
|
---|
158 |
|
---|
159 | // now prepare operations on CP structure
|
---|
160 | cIDs_to_remove.resize(0);
|
---|
161 | cIDs_to_remove.push_back(coordIDs[jet_i].orig);
|
---|
162 | cIDs_to_remove.push_back(coordIDs[jet_j].orig);
|
---|
163 | if (coordIDs[jet_i].mirror != Invalid)
|
---|
164 | cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
|
---|
165 | if (coordIDs[jet_j].mirror != Invalid)
|
---|
166 | cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
|
---|
167 |
|
---|
168 | Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
|
---|
169 | new_points.resize(0);
|
---|
170 | new_points.push_back(new_point);
|
---|
171 | if (make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point); //< same warning as before concerning the mirroring
|
---|
172 |
|
---|
173 | // carry out actions on search tree
|
---|
174 | cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
|
---|
175 |
|
---|
176 | // now fill in info for new points...
|
---|
177 | coordIDs[newjet_k].orig = new_cIDs[0];
|
---|
178 | jetIDs[new_cIDs[0]] = newjet_k;
|
---|
179 | if (new_cIDs.size() == 2) {
|
---|
180 | coordIDs[newjet_k].mirror = new_cIDs[1];
|
---|
181 | jetIDs[new_cIDs[1]] = newjet_k;
|
---|
182 | } else {coordIDs[newjet_k].mirror = Invalid;}
|
---|
183 |
|
---|
184 | //// if we've reached one "active" jet we should exit...
|
---|
185 | //n_active--;
|
---|
186 | //if (n_active == 1) {break;}
|
---|
187 |
|
---|
188 | } while(true);
|
---|
189 |
|
---|
190 | // cross check to see what's going on...
|
---|
191 | //cerr << "Max tree depths after: ";
|
---|
192 | //cp.print_tree_depths(cerr);
|
---|
193 |
|
---|
194 | }
|
---|
195 |
|
---|
196 |
|
---|
197 | //----------------------------------------------------------------------
|
---|
198 | /// a variant of the closest pair clustering which uses a region of
|
---|
199 | /// size 2pi+2R in phi.
|
---|
200 | void ClusterSequence::_CP2DChan_cluster_2pi2R () {
|
---|
201 |
|
---|
202 | if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
|
---|
203 |
|
---|
204 | // run the clustering with mirror copies kept such that only things
|
---|
205 | // within _Rparam of a border are mirrored
|
---|
206 | _CP2DChan_limited_cluster(_Rparam);
|
---|
207 |
|
---|
208 | //
|
---|
209 | _do_Cambridge_inclusive_jets();
|
---|
210 | }
|
---|
211 |
|
---|
212 |
|
---|
213 | //----------------------------------------------------------------------
|
---|
214 | /// a variant of the closest pair clustering which uses a region of
|
---|
215 | /// size 2pi + 2*0.3 and then carries on with 2pi+2*R
|
---|
216 | void ClusterSequence::_CP2DChan_cluster_2piMultD () {
|
---|
217 |
|
---|
218 | // do a first run of clustering up to a small distance parameter,
|
---|
219 | if (_Rparam >= 0.39) {
|
---|
220 | _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
|
---|
221 | }
|
---|
222 |
|
---|
223 | // and then the final round of clustering
|
---|
224 | _CP2DChan_cluster_2pi2R ();
|
---|
225 | }
|
---|
226 |
|
---|
227 |
|
---|
228 | //----------------------------------------------------------------------
|
---|
229 | /// a 4pi variant of the closest pair clustering
|
---|
230 | void ClusterSequence::_CP2DChan_cluster () {
|
---|
231 |
|
---|
232 | if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
|
---|
233 |
|
---|
234 | unsigned int n = _jets.size();
|
---|
235 |
|
---|
236 | vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices
|
---|
237 | vector<int> jetIDs(2*n); // link from mirror to original indices
|
---|
238 | vector<Coord2D> coords(2*n); // our coordinates (and copies)
|
---|
239 |
|
---|
240 | // start things off...
|
---|
241 | double minrap = numeric_limits<double>::max();
|
---|
242 | double maxrap = -minrap;
|
---|
243 | int coord_index = 0;
|
---|
244 | for (unsigned i = 0; i < n; i++) {
|
---|
245 | // separate out points with infinite rapidity
|
---|
246 | if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
|
---|
247 | coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
|
---|
248 | } else {
|
---|
249 | coordIDs[i].orig = coord_index;
|
---|
250 | coordIDs[i].mirror = coord_index+1;
|
---|
251 | coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
|
---|
252 | coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
|
---|
253 | jetIDs[coord_index] = i;
|
---|
254 | jetIDs[coord_index+1] = i;
|
---|
255 | minrap = min(coords[coord_index].x,minrap);
|
---|
256 | maxrap = max(coords[coord_index].x,maxrap);
|
---|
257 | coord_index += 2;
|
---|
258 | }
|
---|
259 | }
|
---|
260 | // label remaining "mirror info" as saying that there are no jets
|
---|
261 | for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
|
---|
262 |
|
---|
263 | // set them to their actual size...
|
---|
264 | coords.resize(coord_index);
|
---|
265 |
|
---|
266 | // establish limits (with some leeway on rapidity)
|
---|
267 | Coord2D left_edge(minrap-1.0, 0.0);
|
---|
268 | Coord2D right_edge(maxrap+1.0, 2*twopi);
|
---|
269 |
|
---|
270 | // now create the closest pair search object
|
---|
271 | ClosestPair2D cp(coords, left_edge, right_edge);
|
---|
272 |
|
---|
273 | // and start iterating...
|
---|
274 | vector<Coord2D> new_points(2);
|
---|
275 | vector<unsigned int> cIDs_to_remove(4);
|
---|
276 | vector<unsigned int> new_cIDs(2);
|
---|
277 | do {
|
---|
278 | // find out which pair we recombine
|
---|
279 | unsigned int cID1, cID2;
|
---|
280 | double distance2;
|
---|
281 | cp.closest_pair(cID1,cID2,distance2);
|
---|
282 | distance2 *= _invR2;
|
---|
283 |
|
---|
284 | // if the closest distance > 1, we exit and go to "inclusive" stage
|
---|
285 | if (distance2 > 1.0) {break;}
|
---|
286 |
|
---|
287 | // do the recombination...
|
---|
288 | int jet_i = jetIDs[cID1];
|
---|
289 | int jet_j = jetIDs[cID2];
|
---|
290 | assert (jet_i != jet_j); // to catch issue of recombining with mirror point
|
---|
291 | int newjet_k;
|
---|
292 | _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
|
---|
293 |
|
---|
294 | // now prepare operations on CP structure
|
---|
295 | cIDs_to_remove[0] = coordIDs[jet_i].orig;
|
---|
296 | cIDs_to_remove[1] = coordIDs[jet_i].mirror;
|
---|
297 | cIDs_to_remove[2] = coordIDs[jet_j].orig;
|
---|
298 | cIDs_to_remove[3] = coordIDs[jet_j].mirror;
|
---|
299 | new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
|
---|
300 | new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
|
---|
301 | // carry out the CP operation
|
---|
302 | //cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
|
---|
303 | // remarkable the following is faster...
|
---|
304 | new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
|
---|
305 | new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
|
---|
306 | // signal that the following jets no longer exist
|
---|
307 | coordIDs[jet_i].orig = Invalid;
|
---|
308 | coordIDs[jet_j].orig = Invalid;
|
---|
309 | // and do bookkeeping for new jet
|
---|
310 | coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
|
---|
311 | jetIDs[new_cIDs[0]] = newjet_k;
|
---|
312 | jetIDs[new_cIDs[1]] = newjet_k;
|
---|
313 |
|
---|
314 | // if we've reached one jet we should exit...
|
---|
315 | n--;
|
---|
316 | if (n == 1) {break;}
|
---|
317 |
|
---|
318 | } while(true);
|
---|
319 |
|
---|
320 |
|
---|
321 | // now do "beam" recombinations
|
---|
322 | //for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) {
|
---|
323 | // // if we have a predeclared beam jet or a valid set of mirror
|
---|
324 | // // coordinates, recombine then recombine this jet with the beam
|
---|
325 | // if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) {
|
---|
326 | // // diB = 1 _always_ (for the cambridge alg.)
|
---|
327 | // _do_iB_recombination_step(jet_i, 1.0);
|
---|
328 | // }
|
---|
329 | //}
|
---|
330 |
|
---|
331 | _do_Cambridge_inclusive_jets();
|
---|
332 |
|
---|
333 | //n = _history.size();
|
---|
334 | //for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
|
---|
335 | // if (_history[hist_i].child == Invalid) {
|
---|
336 | // _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
|
---|
337 | // }
|
---|
338 | //}
|
---|
339 |
|
---|
340 | }
|
---|
341 |
|
---|
342 |
|
---|
343 | //----------------------------------------------------------------------
|
---|
344 | void ClusterSequence::_do_Cambridge_inclusive_jets () {
|
---|
345 | unsigned int n = _history.size();
|
---|
346 | for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
|
---|
347 | if (_history[hist_i].child == Invalid) {
|
---|
348 | _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
|
---|
349 | }
|
---|
350 | }
|
---|
351 | }
|
---|
352 |
|
---|
353 | FASTJET_END_NAMESPACE
|
---|
354 |
|
---|