00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/ClusterSequence.hh"
00032 #include "fastjet/internal/ClosestPair2D.hh"
00033 #include<limits>
00034 #include<vector>
00035 #include<cmath>
00036 #include<iostream>
00037
00038 FASTJET_BEGIN_NAMESPACE
00039
00040 using namespace std;
00041
00042
00043 namespace Private {
00045 class MirrorInfo{
00046 public:
00047 int orig, mirror;
00048 MirrorInfo(int a, int b) : orig(a), mirror(b) {};
00049 MirrorInfo() {};
00050 };
00051
00055 bool make_mirror(Coord2D & point, double Dlim) {
00056 if (point.y < Dlim) {point.y += twopi; return true;}
00057 if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
00058 return false;
00059 }
00060
00061 }
00062
00063 using namespace Private;
00064
00065
00066
00069 void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) {
00070
00071 unsigned int n = _initial_n;
00072
00073 vector<MirrorInfo> coordIDs(2*n);
00074 vector<int> jetIDs(2*n);
00075 vector<Coord2D> coords(2*n);
00076
00077
00078 double minrap = numeric_limits<double>::max();
00079 double maxrap = -minrap;
00080 int coord_index = -1;
00081 int n_active = 0;
00082 for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
00083
00084
00085
00086 if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
00087 (_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
00088 _jets[jet_i].perp2() == 0.0)
00089 ) {continue;}
00090
00091 n_active++;
00092
00093 coordIDs[jet_i].orig = ++coord_index;
00094 coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
00095 jetIDs[coord_index] = jet_i;
00096 minrap = min(coords[coord_index].x,minrap);
00097 maxrap = max(coords[coord_index].x,maxrap);
00098
00099 Coord2D mirror_point(coords[coord_index]);
00100 if (make_mirror(mirror_point, Dlim)) {
00101 coordIDs[jet_i].mirror = ++coord_index;
00102 coords[coord_index] = mirror_point;
00103 jetIDs[coord_index] = jet_i;
00104 } else {
00105 coordIDs[jet_i].mirror = Invalid;
00106 }
00107 }
00108
00109
00110 coords.resize(coord_index+1);
00111
00112
00113 Coord2D left_edge(minrap-1.0, -pi);
00114 Coord2D right_edge(maxrap+1.0, 3*pi);
00115
00116
00117
00118
00119 ClosestPair2D cp(coords, left_edge, right_edge);
00120
00121
00122
00123
00124
00125
00126 vector<Coord2D> new_points(2);
00127 vector<unsigned int> cIDs_to_remove(4);
00128 vector<unsigned int> new_cIDs(2);
00129
00130 do {
00131
00132 unsigned int cID1, cID2;
00133 double distance2;
00134 cp.closest_pair(cID1,cID2,distance2);
00135
00136
00137 if (distance2 > Dlim*Dlim) {break;}
00138
00139
00140 distance2 *= _invR2;
00141
00142
00143 int jet_i = jetIDs[cID1];
00144 int jet_j = jetIDs[cID2];
00145 int newjet_k;
00146 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00147
00148
00149
00150 if (--n_active == 1) {break;}
00151
00152
00153 cIDs_to_remove.resize(0);
00154 cIDs_to_remove.push_back(coordIDs[jet_i].orig);
00155 cIDs_to_remove.push_back(coordIDs[jet_j].orig);
00156 if (coordIDs[jet_i].mirror != Invalid)
00157 cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
00158 if (coordIDs[jet_j].mirror != Invalid)
00159 cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
00160
00161 Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00162 new_points.resize(0);
00163 new_points.push_back(new_point);
00164 if (make_mirror(new_point, Dlim)) new_points.push_back(new_point);
00165
00166
00167 cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
00168
00169
00170 coordIDs[newjet_k].orig = new_cIDs[0];
00171 jetIDs[new_cIDs[0]] = newjet_k;
00172 if (new_cIDs.size() == 2) {
00173 coordIDs[newjet_k].mirror = new_cIDs[1];
00174 jetIDs[new_cIDs[1]] = newjet_k;
00175 } else {coordIDs[newjet_k].mirror = Invalid;}
00176
00178
00179
00180
00181 } while(true);
00182
00183
00184
00185
00186
00187 }
00188
00189
00190
00193 void ClusterSequence::_CP2DChan_cluster_2pi2R () {
00194
00195 if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
00196
00197
00198
00199 _CP2DChan_limited_cluster(_Rparam);
00200
00201
00202 _do_Cambridge_inclusive_jets();
00203 }
00204
00205
00206
00209 void ClusterSequence::_CP2DChan_cluster_2piMultD () {
00210
00211
00212 if (_Rparam >= 0.39) {
00213 _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
00214 }
00215
00216
00217 _CP2DChan_cluster_2pi2R ();
00218 }
00219
00220
00221
00223 void ClusterSequence::_CP2DChan_cluster () {
00224
00225 if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
00226
00227 unsigned int n = _jets.size();
00228
00229 vector<MirrorInfo> coordIDs(2*n);
00230 vector<int> jetIDs(2*n);
00231 vector<Coord2D> coords(2*n);
00232
00233
00234 double minrap = numeric_limits<double>::max();
00235 double maxrap = -minrap;
00236 int coord_index = 0;
00237 for (unsigned i = 0; i < n; i++) {
00238
00239 if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
00240 coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
00241 } else {
00242 coordIDs[i].orig = coord_index;
00243 coordIDs[i].mirror = coord_index+1;
00244 coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
00245 coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
00246 jetIDs[coord_index] = i;
00247 jetIDs[coord_index+1] = i;
00248 minrap = min(coords[coord_index].x,minrap);
00249 maxrap = max(coords[coord_index].x,maxrap);
00250 coord_index += 2;
00251 }
00252 }
00253
00254 for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
00255
00256
00257 coords.resize(coord_index);
00258
00259
00260 Coord2D left_edge(minrap-1.0, 0.0);
00261 Coord2D right_edge(maxrap+1.0, 2*twopi);
00262
00263
00264 ClosestPair2D cp(coords, left_edge, right_edge);
00265
00266
00267 vector<Coord2D> new_points(2);
00268 vector<unsigned int> cIDs_to_remove(4);
00269 vector<unsigned int> new_cIDs(2);
00270 do {
00271
00272 unsigned int cID1, cID2;
00273 double distance2;
00274 cp.closest_pair(cID1,cID2,distance2);
00275 distance2 *= _invR2;
00276
00277
00278 if (distance2 > 1.0) {break;}
00279
00280
00281 int jet_i = jetIDs[cID1];
00282 int jet_j = jetIDs[cID2];
00283 int newjet_k;
00284 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00285
00286
00287 cIDs_to_remove[0] = coordIDs[jet_i].orig;
00288 cIDs_to_remove[1] = coordIDs[jet_i].mirror;
00289 cIDs_to_remove[2] = coordIDs[jet_j].orig;
00290 cIDs_to_remove[3] = coordIDs[jet_j].mirror;
00291 new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00292 new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
00293
00294
00295
00296 new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
00297 new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
00298
00299 coordIDs[jet_i].orig = Invalid;
00300 coordIDs[jet_j].orig = Invalid;
00301
00302 coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
00303 jetIDs[new_cIDs[0]] = newjet_k;
00304 jetIDs[new_cIDs[1]] = newjet_k;
00305
00306
00307 n--;
00308 if (n == 1) {break;}
00309
00310 } while(true);
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323 _do_Cambridge_inclusive_jets();
00324
00325
00326
00327
00328
00329
00330
00331
00332 }
00333
00334
00335
00336 void ClusterSequence::_do_Cambridge_inclusive_jets () {
00337 unsigned int n = _history.size();
00338 for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
00339 if (_history[hist_i].child == Invalid) {
00340 _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
00341 }
00342 }
00343 }
00344
00345 FASTJET_END_NAMESPACE
00346