29 #include "fastjet/ClusterSequence.hh" 30 #include "fastjet/internal/ClosestPair2D.hh" 36 FASTJET_BEGIN_NAMESPACE
46 MirrorInfo(
int a,
int b) : orig(a), mirror(b) {};
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;}
61 using namespace Private;
67 void ClusterSequence::_CP2DChan_limited_cluster (
double Dlim) {
69 unsigned int n = _initial_n;
71 vector<MirrorInfo> coordIDs(2*n);
72 vector<int> jetIDs(2*n);
73 vector<Coord2D> coords(2*n);
81 double Dlim4mirror = min(Dlim,pi);
84 double minrap = numeric_limits<double>::max();
85 double maxrap = -minrap;
88 for (
unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
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)
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);
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;
111 coordIDs[jet_i].mirror = Invalid;
116 coords.resize(coord_index+1);
119 Coord2D left_edge(minrap-1.0, -3.15);
120 Coord2D right_edge(maxrap+1.0, 9.45);
125 ClosestPair2D cp(coords, left_edge, right_edge);
132 vector<Coord2D> new_points(2);
133 vector<unsigned int> cIDs_to_remove(4);
134 vector<unsigned int> new_cIDs(2);
138 unsigned int cID1, cID2;
140 cp.closest_pair(cID1,cID2,distance2);
143 if (distance2 > Dlim*Dlim) {
break;}
149 int jet_i = jetIDs[cID1];
150 int jet_j = jetIDs[cID2];
151 assert (jet_i != jet_j);
153 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
157 if (--n_active == 1) {
break;}
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);
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);
174 cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
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;}
200 void ClusterSequence::_CP2DChan_cluster_2pi2R () {
202 if (_jet_algorithm != cambridge_algorithm)
throw Error(
"CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
206 _CP2DChan_limited_cluster(_Rparam);
209 _do_Cambridge_inclusive_jets();
216 void ClusterSequence::_CP2DChan_cluster_2piMultD () {
219 if (_Rparam >= 0.39) {
220 _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
224 _CP2DChan_cluster_2pi2R ();
230 void ClusterSequence::_CP2DChan_cluster () {
232 if (_jet_algorithm != cambridge_algorithm)
throw Error(
"_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
234 unsigned int n = _jets.size();
236 vector<MirrorInfo> coordIDs(2*n);
237 vector<int> jetIDs(2*n);
238 vector<Coord2D> coords(2*n);
241 double minrap = numeric_limits<double>::max();
242 double maxrap = -minrap;
244 for (
unsigned i = 0; i < n; i++) {
246 if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
247 coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
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);
261 for (
unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
264 coords.resize(coord_index);
267 Coord2D left_edge(minrap-1.0, 0.0);
268 Coord2D right_edge(maxrap+1.0, 2*twopi);
271 ClosestPair2D cp(coords, left_edge, right_edge);
274 vector<Coord2D> new_points(2);
275 vector<unsigned int> cIDs_to_remove(4);
276 vector<unsigned int> new_cIDs(2);
279 unsigned int cID1, cID2;
281 cp.closest_pair(cID1,cID2,distance2);
285 if (distance2 > 1.0) {
break;}
288 int jet_i = jetIDs[cID1];
289 int jet_j = jetIDs[cID2];
290 assert (jet_i != jet_j);
292 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
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);
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]);
307 coordIDs[jet_i].orig = Invalid;
308 coordIDs[jet_j].orig = Invalid;
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;
331 _do_Cambridge_inclusive_jets();
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);
353 FASTJET_END_NAMESPACE