FastJet  3.2.1
ClusterSequence_CP2DChan.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence_CP2DChan.cc 4045 2016-03-03 10:01:55Z 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 #include "fastjet/ClusterSequence.hh"
32 #include "fastjet/internal/ClosestPair2D.hh"
33 #include<limits>
34 #include<vector>
35 #include<cmath>
36 #include<iostream>
37 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 using namespace std;
41 
42 // place for things we don't want outside world to run into
43 namespace 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 
63 using 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
69 void 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.
202 void 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
218 void 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
232 void 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 //----------------------------------------------------------------------
346 void 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 
355 FASTJET_END_NAMESPACE
356