FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequence_CP2DChan.cc 2577 2011-09-13 15:11:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 #include "fastjet/ClusterSequence.hh" 00030 #include "fastjet/internal/ClosestPair2D.hh" 00031 #include<limits> 00032 #include<vector> 00033 #include<cmath> 00034 #include<iostream> 00035 00036 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00037 00038 using namespace std; 00039 00040 // place for things we don't want outside world to run into 00041 namespace Private { 00042 /// class for helping us deal with mirror-image particles. 00043 class MirrorInfo{ 00044 public: 00045 int orig, mirror; 00046 MirrorInfo(int a, int b) : orig(a), mirror(b) {}; 00047 MirrorInfo() {}; 00048 }; 00049 00050 /// if there is a need for a mirror when looking for closest pairs 00051 /// up to distance D, then return true and turn the supplied point 00052 /// into its mirror copy 00053 bool make_mirror(Coord2D & point, double Dlim) { 00054 if (point.y < Dlim) {point.y += twopi; return true;} 00055 if (twopi-point.y < Dlim) {point.y -= twopi; return true;} 00056 return false; 00057 } 00058 00059 } 00060 00061 using namespace Private; 00062 00063 00064 //---------------------------------------------------------------------- 00065 /// clusters only up to a distance Dlim -- does not deal with "inclusive" jets 00066 /// -- these are left to some other part of the program 00067 void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) { 00068 00069 unsigned int n = _initial_n; 00070 00071 vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID 00072 vector<int> jetIDs(2*n); // jet ID for a given coord ID 00073 vector<Coord2D> coords(2*n); // our coordinates (and copies) 00074 00075 // particles within a distance Dlim of the phi edges (phi<Dlim || 00076 // phi>2pi-Dli;) will be mirrored. For Dlim>pi, this could lead to 00077 // particles copies outside the fixed range in phi which is 00078 // [-pi:3pi] (see make_mirror above). Since in that case all 00079 // particles get copied anywaym we can just copy particles up to a 00080 // distance "pi" from the edges 00081 double Dlim4mirror = min(Dlim,pi); 00082 00083 // start things off... 00084 double minrap = numeric_limits<double>::max(); 00085 double maxrap = -minrap; 00086 int coord_index = -1; 00087 int n_active = 0; 00088 for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) { 00089 00090 // skip jets that already have children or that have infinite 00091 // rapidity 00092 if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid || 00093 (_jets[jet_i].E() == abs(_jets[jet_i].pz()) && 00094 _jets[jet_i].perp2() == 0.0) 00095 ) {continue;} 00096 00097 n_active++; 00098 00099 coordIDs[jet_i].orig = ++coord_index; 00100 coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi()); 00101 jetIDs[coord_index] = jet_i; 00102 minrap = min(coords[coord_index].x,minrap); 00103 maxrap = max(coords[coord_index].x,maxrap); 00104 00105 Coord2D mirror_point(coords[coord_index]); 00106 if (make_mirror(mirror_point, Dlim4mirror)) { 00107 coordIDs[jet_i].mirror = ++coord_index; 00108 coords[coord_index] = mirror_point; 00109 jetIDs[coord_index] = jet_i; 00110 } else { 00111 coordIDs[jet_i].mirror = Invalid; 00112 } 00113 } 00114 00115 // set them to their actual size... 00116 coords.resize(coord_index+1); 00117 00118 // establish limits (with some leeway on rapidity) 00119 Coord2D left_edge(minrap-1.0, -3.15); // a security margin below -pi 00120 Coord2D right_edge(maxrap+1.0, 9.45); // a security margin above 3*pi 00121 00122 //cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl; 00123 00124 // now create the closest pair search object 00125 ClosestPair2D cp(coords, left_edge, right_edge); 00126 00127 // cross check to see what's going on... 00128 //cerr << "Tree depths before: "; 00129 //cp.print_tree_depths(cerr); 00130 00131 // and start iterating... 00132 vector<Coord2D> new_points(2); 00133 vector<unsigned int> cIDs_to_remove(4); 00134 vector<unsigned int> new_cIDs(2); 00135 00136 do { 00137 // find out which pair we recombine 00138 unsigned int cID1, cID2; 00139 double distance2; 00140 cp.closest_pair(cID1,cID2,distance2); 00141 00142 // if the closest distance > Dlim, we exit and go to "inclusive" stage 00143 if (distance2 > Dlim*Dlim) {break;} 00144 00145 // normalise distance as we like it 00146 distance2 *= _invR2; 00147 00148 // do the recombination... 00149 int jet_i = jetIDs[cID1]; 00150 int jet_j = jetIDs[cID2]; 00151 assert (jet_i != jet_j); // to catch issue of recombining with mirror point 00152 int newjet_k; 00153 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k); 00154 00155 // don't bother with any further action if only one active particle 00156 // is left (also avoid closest-pair error [cannot remove last particle]). 00157 if (--n_active == 1) {break;} 00158 00159 // now prepare operations on CP structure 00160 cIDs_to_remove.resize(0); 00161 cIDs_to_remove.push_back(coordIDs[jet_i].orig); 00162 cIDs_to_remove.push_back(coordIDs[jet_j].orig); 00163 if (coordIDs[jet_i].mirror != Invalid) 00164 cIDs_to_remove.push_back(coordIDs[jet_i].mirror); 00165 if (coordIDs[jet_j].mirror != Invalid) 00166 cIDs_to_remove.push_back(coordIDs[jet_j].mirror); 00167 00168 Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()); 00169 new_points.resize(0); 00170 new_points.push_back(new_point); 00171 if (make_mirror(new_point, Dlim4mirror)) new_points.push_back(new_point); //< same warning as before concerning the mirroring 00172 00173 // carry out actions on search tree 00174 cp.replace_many(cIDs_to_remove, new_points, new_cIDs); 00175 00176 // now fill in info for new points... 00177 coordIDs[newjet_k].orig = new_cIDs[0]; 00178 jetIDs[new_cIDs[0]] = newjet_k; 00179 if (new_cIDs.size() == 2) { 00180 coordIDs[newjet_k].mirror = new_cIDs[1]; 00181 jetIDs[new_cIDs[1]] = newjet_k; 00182 } else {coordIDs[newjet_k].mirror = Invalid;} 00183 00184 //// if we've reached one "active" jet we should exit... 00185 //n_active--; 00186 //if (n_active == 1) {break;} 00187 00188 } while(true); 00189 00190 // cross check to see what's going on... 00191 //cerr << "Max tree depths after: "; 00192 //cp.print_tree_depths(cerr); 00193 00194 } 00195 00196 00197 //---------------------------------------------------------------------- 00198 /// a variant of the closest pair clustering which uses a region of 00199 /// size 2pi+2R in phi. 00200 void ClusterSequence::_CP2DChan_cluster_2pi2R () { 00201 00202 if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm"); 00203 00204 // run the clustering with mirror copies kept such that only things 00205 // within _Rparam of a border are mirrored 00206 _CP2DChan_limited_cluster(_Rparam); 00207 00208 // 00209 _do_Cambridge_inclusive_jets(); 00210 } 00211 00212 00213 //---------------------------------------------------------------------- 00214 /// a variant of the closest pair clustering which uses a region of 00215 /// size 2pi + 2*0.3 and then carries on with 2pi+2*R 00216 void ClusterSequence::_CP2DChan_cluster_2piMultD () { 00217 00218 // do a first run of clustering up to a small distance parameter, 00219 if (_Rparam >= 0.39) { 00220 _CP2DChan_limited_cluster(min(_Rparam/2,0.3)); 00221 } 00222 00223 // and then the final round of clustering 00224 _CP2DChan_cluster_2pi2R (); 00225 } 00226 00227 00228 //---------------------------------------------------------------------- 00229 /// a 4pi variant of the closest pair clustering 00230 void ClusterSequence::_CP2DChan_cluster () { 00231 00232 if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm"); 00233 00234 unsigned int n = _jets.size(); 00235 00236 vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices 00237 vector<int> jetIDs(2*n); // link from mirror to original indices 00238 vector<Coord2D> coords(2*n); // our coordinates (and copies) 00239 00240 // start things off... 00241 double minrap = numeric_limits<double>::max(); 00242 double maxrap = -minrap; 00243 int coord_index = 0; 00244 for (unsigned i = 0; i < n; i++) { 00245 // separate out points with infinite rapidity 00246 if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) { 00247 coordIDs[i] = MirrorInfo(BeamJet,BeamJet); 00248 } else { 00249 coordIDs[i].orig = coord_index; 00250 coordIDs[i].mirror = coord_index+1; 00251 coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()); 00252 coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi); 00253 jetIDs[coord_index] = i; 00254 jetIDs[coord_index+1] = i; 00255 minrap = min(coords[coord_index].x,minrap); 00256 maxrap = max(coords[coord_index].x,maxrap); 00257 coord_index += 2; 00258 } 00259 } 00260 // label remaining "mirror info" as saying that there are no jets 00261 for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;} 00262 00263 // set them to their actual size... 00264 coords.resize(coord_index); 00265 00266 // establish limits (with some leeway on rapidity) 00267 Coord2D left_edge(minrap-1.0, 0.0); 00268 Coord2D right_edge(maxrap+1.0, 2*twopi); 00269 00270 // now create the closest pair search object 00271 ClosestPair2D cp(coords, left_edge, right_edge); 00272 00273 // and start iterating... 00274 vector<Coord2D> new_points(2); 00275 vector<unsigned int> cIDs_to_remove(4); 00276 vector<unsigned int> new_cIDs(2); 00277 do { 00278 // find out which pair we recombine 00279 unsigned int cID1, cID2; 00280 double distance2; 00281 cp.closest_pair(cID1,cID2,distance2); 00282 distance2 *= _invR2; 00283 00284 // if the closest distance > 1, we exit and go to "inclusive" stage 00285 if (distance2 > 1.0) {break;} 00286 00287 // do the recombination... 00288 int jet_i = jetIDs[cID1]; 00289 int jet_j = jetIDs[cID2]; 00290 assert (jet_i != jet_j); // to catch issue of recombining with mirror point 00291 int newjet_k; 00292 _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k); 00293 00294 // now prepare operations on CP structure 00295 cIDs_to_remove[0] = coordIDs[jet_i].orig; 00296 cIDs_to_remove[1] = coordIDs[jet_i].mirror; 00297 cIDs_to_remove[2] = coordIDs[jet_j].orig; 00298 cIDs_to_remove[3] = coordIDs[jet_j].mirror; 00299 new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()); 00300 new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi); 00301 // carry out the CP operation 00302 //cp.replace_many(cIDs_to_remove, new_points, new_cIDs); 00303 // remarkable the following is faster... 00304 new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]); 00305 new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]); 00306 // signal that the following jets no longer exist 00307 coordIDs[jet_i].orig = Invalid; 00308 coordIDs[jet_j].orig = Invalid; 00309 // and do bookkeeping for new jet 00310 coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]); 00311 jetIDs[new_cIDs[0]] = newjet_k; 00312 jetIDs[new_cIDs[1]] = newjet_k; 00313 00314 // if we've reached one jet we should exit... 00315 n--; 00316 if (n == 1) {break;} 00317 00318 } while(true); 00319 00320 00321 // now do "beam" recombinations 00322 //for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) { 00323 // // if we have a predeclared beam jet or a valid set of mirror 00324 // // coordinates, recombine then recombine this jet with the beam 00325 // if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) { 00326 // // diB = 1 _always_ (for the cambridge alg.) 00327 // _do_iB_recombination_step(jet_i, 1.0); 00328 // } 00329 //} 00330 00331 _do_Cambridge_inclusive_jets(); 00332 00333 //n = _history.size(); 00334 //for (unsigned int hist_i = 0; hist_i < n; hist_i++) { 00335 // if (_history[hist_i].child == Invalid) { 00336 // _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0); 00337 // } 00338 //} 00339 00340 } 00341 00342 00343 //---------------------------------------------------------------------- 00344 void ClusterSequence::_do_Cambridge_inclusive_jets () { 00345 unsigned int n = _history.size(); 00346 for (unsigned int hist_i = 0; hist_i < n; hist_i++) { 00347 if (_history[hist_i].child == Invalid) { 00348 _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0); 00349 } 00350 } 00351 } 00352 00353 FASTJET_END_NAMESPACE 00354