FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ClusterSequence_CP2DChan.cc
1 //STARTHEADER
2 // $Id: ClusterSequence_CP2DChan.cc 2577 2011-09-13 15:11:38Z salam $
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