Fork me on GitHub

source: git/external/fastjet/ClusterSequence_CP2DChan.cc@ 0e6b14e

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