| 1 | //FJSTARTHEADER
|
|---|
| 2 | // $Id: ClusterSequence_CP2DChan.cc 3433 2014-07-23 08:17:03Z salam $
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
|
|---|
| 16 | // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
|
|---|
| 17 | // FastJet as part of work towards a scientific publication, please
|
|---|
| 18 | // quote the version you use and include a citation to the manual and
|
|---|
| 19 | // optionally also to hep-ph/0512210.
|
|---|
| 20 | //
|
|---|
| 21 | // FastJet is distributed in the hope that it will be useful,
|
|---|
| 22 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 23 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 24 | // GNU General Public License for more details.
|
|---|
| 25 | //
|
|---|
| 26 | // You should have received a copy of the GNU General Public License
|
|---|
| 27 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 28 | //----------------------------------------------------------------------
|
|---|
| 29 | //FJENDHEADER
|
|---|
| 30 |
|
|---|
| 31 | #include "fastjet/ClusterSequence.hh"
|
|---|
| 32 | #include "fastjet/internal/ClosestPair2D.hh"
|
|---|
| 33 | #include<limits>
|
|---|
| 34 | #include<vector>
|
|---|
| 35 | #include<cmath>
|
|---|
| 36 | #include<iostream>
|
|---|
| 37 |
|
|---|
| 38 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 39 |
|
|---|
| 40 | using namespace std;
|
|---|
| 41 |
|
|---|
| 42 | // place for things we don't want outside world to run into
|
|---|
| 43 | namespace Private {
|
|---|
| 44 | /// class for helping us deal with mirror-image particles.
|
|---|
| 45 | class MirrorInfo{
|
|---|
| 46 | public:
|
|---|
| 47 | int orig, mirror;
|
|---|
| 48 | MirrorInfo(int a, int b) : orig(a), mirror(b) {};
|
|---|
| 49 | MirrorInfo() {};
|
|---|
| 50 | };
|
|---|
| 51 |
|
|---|
| 52 | /// if there is a need for a mirror when looking for closest pairs
|
|---|
| 53 | /// up to distance D, then return true and turn the supplied point
|
|---|
| 54 | /// into its mirror copy
|
|---|
| 55 | bool make_mirror(Coord2D & point, double Dlim) {
|
|---|
| 56 | if (point.y < Dlim) {point.y += twopi; return true;}
|
|---|
| 57 | if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
|
|---|
| 58 | return false;
|
|---|
| 59 | }
|
|---|
| 60 |
|
|---|
| 61 | }
|
|---|
| 62 |
|
|---|
| 63 | using namespace Private;
|
|---|
| 64 |
|
|---|
| 65 |
|
|---|
| 66 | //----------------------------------------------------------------------
|
|---|
| 67 | /// clusters only up to a distance Dlim -- does not deal with "inclusive" jets
|
|---|
| 68 | /// -- these are left to some other part of the program
|
|---|
| 69 | void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) {
|
|---|
| 70 |
|
|---|
| 71 | unsigned int n = _initial_n;
|
|---|
| 72 |
|
|---|
| 73 | vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID
|
|---|
| 74 | vector<int> jetIDs(2*n); // jet ID for a given coord ID
|
|---|
| 75 | vector<Coord2D> coords(2*n); // our coordinates (and copies)
|
|---|
| 76 |
|
|---|
| 77 | // particles within a distance Dlim of the phi edges (phi<Dlim ||
|
|---|
| 78 | // phi>2pi-Dli;) will be mirrored. For Dlim>pi, this could lead to
|
|---|
| 79 | // particles copies outside the fixed range in phi which is
|
|---|
| 80 | // [-pi:3pi] (see make_mirror above). Since in that case all
|
|---|
| 81 | // particles get copied anywaym we can just copy particles up to a
|
|---|
| 82 | // distance "pi" from the edges
|
|---|
| 83 | double Dlim4mirror = min(Dlim,pi);
|
|---|
| 84 |
|
|---|
| 85 | // start things off...
|
|---|
| 86 | double minrap = numeric_limits<double>::max();
|
|---|
| 87 | double maxrap = -minrap;
|
|---|
| 88 | int coord_index = -1;
|
|---|
| 89 | int n_active = 0;
|
|---|
| 90 | for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
|
|---|
| 91 |
|
|---|
| 92 | // skip jets that already have children or that have infinite
|
|---|
| 93 | // rapidity
|
|---|
| 94 | if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
|
|---|
| 95 | (_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
|
|---|
| 96 | _jets[jet_i].perp2() == 0.0)
|
|---|
| 97 | ) {continue;}
|
|---|
| 98 |
|
|---|
| 99 | n_active++;
|
|---|
| 100 |
|
|---|
| 101 | coordIDs[jet_i].orig = ++coord_index;
|
|---|
| 102 | coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
|
|---|
| 103 | jetIDs[coord_index] = jet_i;
|
|---|
| 104 | minrap = min(coords[coord_index].x,minrap);
|
|---|
| 105 | maxrap = max(coords[coord_index].x,maxrap);
|
|---|
| 106 |
|
|---|
| 107 | Coord2D mirror_point(coords[coord_index]);
|
|---|
| 108 | if (make_mirror(mirror_point, Dlim4mirror)) {
|
|---|
| 109 | coordIDs[jet_i].mirror = ++coord_index;
|
|---|
| 110 | coords[coord_index] = mirror_point;
|
|---|
| 111 | jetIDs[coord_index] = jet_i;
|
|---|
| 112 | } else {
|
|---|
| 113 | coordIDs[jet_i].mirror = Invalid;
|
|---|
| 114 | }
|
|---|
| 115 | }
|
|---|
| 116 |
|
|---|
| 117 | // set them to their actual size...
|
|---|
| 118 | coords.resize(coord_index+1);
|
|---|
| 119 |
|
|---|
| 120 | // establish limits (with some leeway on rapidity)
|
|---|
| 121 | Coord2D left_edge(minrap-1.0, -3.15); // a security margin below -pi
|
|---|
| 122 | Coord2D right_edge(maxrap+1.0, 9.45); // a security margin above 3*pi
|
|---|
| 123 |
|
|---|
| 124 | //cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl;
|
|---|
| 125 |
|
|---|
| 126 | // now create the closest pair search object
|
|---|
| 127 | ClosestPair2D cp(coords, left_edge, right_edge);
|
|---|
| 128 |
|
|---|
| 129 | // cross check to see what's going on...
|
|---|
| 130 | //cerr << "Tree depths before: ";
|
|---|
| 131 | //cp.print_tree_depths(cerr);
|
|---|
| 132 |
|
|---|
| 133 | // and start iterating...
|
|---|
| 134 | vector<Coord2D> new_points(2);
|
|---|
| 135 | vector<unsigned int> cIDs_to_remove(4);
|
|---|
| 136 | vector<unsigned int> new_cIDs(2);
|
|---|
| 137 |
|
|---|
| 138 | do {
|
|---|
| 139 | // find out which pair we recombine
|
|---|
| 140 | unsigned int cID1, cID2;
|
|---|
| 141 | double distance2;
|
|---|
| 142 | cp.closest_pair(cID1,cID2,distance2);
|
|---|
| 143 |
|
|---|
| 144 | // if the closest distance > Dlim, we exit and go to "inclusive" stage
|
|---|
| 145 | if (distance2 > Dlim*Dlim) {break;}
|
|---|
| 146 |
|
|---|
| 147 | // normalise distance as we like it
|
|---|
| 148 | distance2 *= _invR2;
|
|---|
| 149 |
|
|---|
| 150 | // do the recombination...
|
|---|
| 151 | int jet_i = jetIDs[cID1];
|
|---|
| 152 | int jet_j = jetIDs[cID2];
|
|---|
| 153 | assert (jet_i != jet_j); // to catch issue of recombining with mirror point
|
|---|
| 154 | int newjet_k;
|
|---|
| 155 | _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
|
|---|
| 156 |
|
|---|
| 157 | // don't bother with any further action if only one active particle
|
|---|
| 158 | // is left (also avoid closest-pair error [cannot remove last particle]).
|
|---|
| 159 | if (--n_active == 1) {break;}
|
|---|
| 160 |
|
|---|
| 161 | // now prepare operations on CP structure
|
|---|
| 162 | cIDs_to_remove.resize(0);
|
|---|
| 163 | cIDs_to_remove.push_back(coordIDs[jet_i].orig);
|
|---|
| 164 | cIDs_to_remove.push_back(coordIDs[jet_j].orig);
|
|---|
| 165 | if (coordIDs[jet_i].mirror != Invalid)
|
|---|
| 166 | cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
|
|---|
| 167 | if (coordIDs[jet_j].mirror != Invalid)
|
|---|
| 168 | cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
|
|---|
| 169 |
|
|---|
| 170 | Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
|
|---|
| 171 | new_points.resize(0);
|
|---|
| 172 | new_points.push_back(new_point);
|
|---|
| 173 | if (make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point); //< same warning as before concerning the mirroring
|
|---|
| 174 |
|
|---|
| 175 | // carry out actions on search tree
|
|---|
| 176 | cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
|
|---|
| 177 |
|
|---|
| 178 | // now fill in info for new points...
|
|---|
| 179 | coordIDs[newjet_k].orig = new_cIDs[0];
|
|---|
| 180 | jetIDs[new_cIDs[0]] = newjet_k;
|
|---|
| 181 | if (new_cIDs.size() == 2) {
|
|---|
| 182 | coordIDs[newjet_k].mirror = new_cIDs[1];
|
|---|
| 183 | jetIDs[new_cIDs[1]] = newjet_k;
|
|---|
| 184 | } else {coordIDs[newjet_k].mirror = Invalid;}
|
|---|
| 185 |
|
|---|
| 186 | //// if we've reached one "active" jet we should exit...
|
|---|
| 187 | //n_active--;
|
|---|
| 188 | //if (n_active == 1) {break;}
|
|---|
| 189 |
|
|---|
| 190 | } while(true);
|
|---|
| 191 |
|
|---|
| 192 | // cross check to see what's going on...
|
|---|
| 193 | //cerr << "Max tree depths after: ";
|
|---|
| 194 | //cp.print_tree_depths(cerr);
|
|---|
| 195 |
|
|---|
| 196 | }
|
|---|
| 197 |
|
|---|
| 198 |
|
|---|
| 199 | //----------------------------------------------------------------------
|
|---|
| 200 | /// a variant of the closest pair clustering which uses a region of
|
|---|
| 201 | /// size 2pi+2R in phi.
|
|---|
| 202 | void ClusterSequence::_CP2DChan_cluster_2pi2R () {
|
|---|
| 203 |
|
|---|
| 204 | if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
|
|---|
| 205 |
|
|---|
| 206 | // run the clustering with mirror copies kept such that only things
|
|---|
| 207 | // within _Rparam of a border are mirrored
|
|---|
| 208 | _CP2DChan_limited_cluster(_Rparam);
|
|---|
| 209 |
|
|---|
| 210 | //
|
|---|
| 211 | _do_Cambridge_inclusive_jets();
|
|---|
| 212 | }
|
|---|
| 213 |
|
|---|
| 214 |
|
|---|
| 215 | //----------------------------------------------------------------------
|
|---|
| 216 | /// a variant of the closest pair clustering which uses a region of
|
|---|
| 217 | /// size 2pi + 2*0.3 and then carries on with 2pi+2*R
|
|---|
| 218 | void ClusterSequence::_CP2DChan_cluster_2piMultD () {
|
|---|
| 219 |
|
|---|
| 220 | // do a first run of clustering up to a small distance parameter,
|
|---|
| 221 | if (_Rparam >= 0.39) {
|
|---|
| 222 | _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | // and then the final round of clustering
|
|---|
| 226 | _CP2DChan_cluster_2pi2R ();
|
|---|
| 227 | }
|
|---|
| 228 |
|
|---|
| 229 |
|
|---|
| 230 | //----------------------------------------------------------------------
|
|---|
| 231 | /// a 4pi variant of the closest pair clustering
|
|---|
| 232 | void ClusterSequence::_CP2DChan_cluster () {
|
|---|
| 233 |
|
|---|
| 234 | if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
|
|---|
| 235 |
|
|---|
| 236 | unsigned int n = _jets.size();
|
|---|
| 237 |
|
|---|
| 238 | vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices
|
|---|
| 239 | vector<int> jetIDs(2*n); // link from mirror to original indices
|
|---|
| 240 | vector<Coord2D> coords(2*n); // our coordinates (and copies)
|
|---|
| 241 |
|
|---|
| 242 | // start things off...
|
|---|
| 243 | double minrap = numeric_limits<double>::max();
|
|---|
| 244 | double maxrap = -minrap;
|
|---|
| 245 | int coord_index = 0;
|
|---|
| 246 | for (unsigned i = 0; i < n; i++) {
|
|---|
| 247 | // separate out points with infinite rapidity
|
|---|
| 248 | if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
|
|---|
| 249 | coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
|
|---|
| 250 | } else {
|
|---|
| 251 | coordIDs[i].orig = coord_index;
|
|---|
| 252 | coordIDs[i].mirror = coord_index+1;
|
|---|
| 253 | coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
|
|---|
| 254 | coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
|
|---|
| 255 | jetIDs[coord_index] = i;
|
|---|
| 256 | jetIDs[coord_index+1] = i;
|
|---|
| 257 | minrap = min(coords[coord_index].x,minrap);
|
|---|
| 258 | maxrap = max(coords[coord_index].x,maxrap);
|
|---|
| 259 | coord_index += 2;
|
|---|
| 260 | }
|
|---|
| 261 | }
|
|---|
| 262 | // label remaining "mirror info" as saying that there are no jets
|
|---|
| 263 | for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
|
|---|
| 264 |
|
|---|
| 265 | // set them to their actual size...
|
|---|
| 266 | coords.resize(coord_index);
|
|---|
| 267 |
|
|---|
| 268 | // establish limits (with some leeway on rapidity)
|
|---|
| 269 | Coord2D left_edge(minrap-1.0, 0.0);
|
|---|
| 270 | Coord2D right_edge(maxrap+1.0, 2*twopi);
|
|---|
| 271 |
|
|---|
| 272 | // now create the closest pair search object
|
|---|
| 273 | ClosestPair2D cp(coords, left_edge, right_edge);
|
|---|
| 274 |
|
|---|
| 275 | // and start iterating...
|
|---|
| 276 | vector<Coord2D> new_points(2);
|
|---|
| 277 | vector<unsigned int> cIDs_to_remove(4);
|
|---|
| 278 | vector<unsigned int> new_cIDs(2);
|
|---|
| 279 | do {
|
|---|
| 280 | // find out which pair we recombine
|
|---|
| 281 | unsigned int cID1, cID2;
|
|---|
| 282 | double distance2;
|
|---|
| 283 | cp.closest_pair(cID1,cID2,distance2);
|
|---|
| 284 | distance2 *= _invR2;
|
|---|
| 285 |
|
|---|
| 286 | // if the closest distance > 1, we exit and go to "inclusive" stage
|
|---|
| 287 | if (distance2 > 1.0) {break;}
|
|---|
| 288 |
|
|---|
| 289 | // do the recombination...
|
|---|
| 290 | int jet_i = jetIDs[cID1];
|
|---|
| 291 | int jet_j = jetIDs[cID2];
|
|---|
| 292 | assert (jet_i != jet_j); // to catch issue of recombining with mirror point
|
|---|
| 293 | int newjet_k;
|
|---|
| 294 | _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
|
|---|
| 295 |
|
|---|
| 296 | // now prepare operations on CP structure
|
|---|
| 297 | cIDs_to_remove[0] = coordIDs[jet_i].orig;
|
|---|
| 298 | cIDs_to_remove[1] = coordIDs[jet_i].mirror;
|
|---|
| 299 | cIDs_to_remove[2] = coordIDs[jet_j].orig;
|
|---|
| 300 | cIDs_to_remove[3] = coordIDs[jet_j].mirror;
|
|---|
| 301 | new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
|
|---|
| 302 | new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
|
|---|
| 303 | // carry out the CP operation
|
|---|
| 304 | //cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
|
|---|
| 305 | // remarkable the following is faster...
|
|---|
| 306 | new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
|
|---|
| 307 | new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
|
|---|
| 308 | // signal that the following jets no longer exist
|
|---|
| 309 | coordIDs[jet_i].orig = Invalid;
|
|---|
| 310 | coordIDs[jet_j].orig = Invalid;
|
|---|
| 311 | // and do bookkeeping for new jet
|
|---|
| 312 | coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
|
|---|
| 313 | jetIDs[new_cIDs[0]] = newjet_k;
|
|---|
| 314 | jetIDs[new_cIDs[1]] = newjet_k;
|
|---|
| 315 |
|
|---|
| 316 | // if we've reached one jet we should exit...
|
|---|
| 317 | n--;
|
|---|
| 318 | if (n == 1) {break;}
|
|---|
| 319 |
|
|---|
| 320 | } while(true);
|
|---|
| 321 |
|
|---|
| 322 |
|
|---|
| 323 | // now do "beam" recombinations
|
|---|
| 324 | //for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) {
|
|---|
| 325 | // // if we have a predeclared beam jet or a valid set of mirror
|
|---|
| 326 | // // coordinates, recombine then recombine this jet with the beam
|
|---|
| 327 | // if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) {
|
|---|
| 328 | // // diB = 1 _always_ (for the cambridge alg.)
|
|---|
| 329 | // _do_iB_recombination_step(jet_i, 1.0);
|
|---|
| 330 | // }
|
|---|
| 331 | //}
|
|---|
| 332 |
|
|---|
| 333 | _do_Cambridge_inclusive_jets();
|
|---|
| 334 |
|
|---|
| 335 | //n = _history.size();
|
|---|
| 336 | //for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
|
|---|
| 337 | // if (_history[hist_i].child == Invalid) {
|
|---|
| 338 | // _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
|
|---|
| 339 | // }
|
|---|
| 340 | //}
|
|---|
| 341 |
|
|---|
| 342 | }
|
|---|
| 343 |
|
|---|
| 344 |
|
|---|
| 345 | //----------------------------------------------------------------------
|
|---|
| 346 | void ClusterSequence::_do_Cambridge_inclusive_jets () {
|
|---|
| 347 | unsigned int n = _history.size();
|
|---|
| 348 | for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
|
|---|
| 349 | if (_history[hist_i].child == Invalid) {
|
|---|
| 350 | _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
|
|---|
| 351 | }
|
|---|
| 352 | }
|
|---|
| 353 | }
|
|---|
| 354 |
|
|---|
| 355 | FASTJET_END_NAMESPACE
|
|---|
| 356 |
|
|---|
