[1] | 1 | //STARTHEADER |
---|
| 2 | // $Id: ClusterSequence_CP2DChan.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
| 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 | |
---|