Fork me on GitHub

source: git/external/fastjet/ClusterSequence_CP2DChan.cc@ 5b822e5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 5b822e5 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 12.2 KB
Line 
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
36FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37
38using namespace std;
39
40// place for things we don't want outside world to run into
41namespace 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
61using 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
67void 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.
200void 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
216void 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
230void 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//----------------------------------------------------------------------
344void 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
353FASTJET_END_NAMESPACE
354
Note: See TracBrowser for help on using the repository browser.