FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequence_Delaunay.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 00030 #include "fastjet/Error.hh" 00031 #include "fastjet/PseudoJet.hh" 00032 #include "fastjet/ClusterSequence.hh" 00033 #include<iostream> 00034 #include<sstream> 00035 #include<cmath> 00036 #include <cstdlib> 00037 #include<cassert> 00038 #include<memory> 00039 // 00040 #ifndef DROP_CGAL // in case we do not have the code for CGAL 00041 #include "fastjet/internal/Dnn4piCylinder.hh" 00042 #include "fastjet/internal/Dnn3piCylinder.hh" 00043 #include "fastjet/internal/Dnn2piCylinder.hh" 00044 #endif // DROP_CGAL 00045 00046 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00047 00048 using namespace std; 00049 00050 00051 //---------------------------------------------------------------------- 00052 /// Run the clustering using a Hierarchical Delaunay triangulation and 00053 /// STL maps to achieve O(N*ln N) behaviour. 00054 /// 00055 /// There may be internally asserted assumptions about absence of 00056 /// points with coincident eta-phi coordinates. 00057 void ClusterSequence::_delaunay_cluster () { 00058 00059 int n = _jets.size(); 00060 00061 vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double> 00062 for (int i = 0; i < n; i++) { 00063 points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi()); 00064 points[i].sanitize(); // make sure things are in the right range 00065 } 00066 00067 // initialise our DNN structure with the set of points 00068 auto_ptr<DynamicNearestNeighbours> DNN; 00069 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL... 00070 bool verbose = false; 00071 bool ignore_nearest_is_mirror = (_Rparam < twopi); 00072 if (_strategy == NlnN4pi) { 00073 DNN.reset(new Dnn4piCylinder(points,verbose)); 00074 } else if (_strategy == NlnN3pi) { 00075 DNN.reset(new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose)); 00076 } else if (_strategy == NlnN) { 00077 DNN.reset(new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose)); 00078 } else 00079 #else 00080 if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) { 00081 ostringstream err; 00082 err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl; 00083 err << " supported because FastJet was compiled without CGAL"<<endl; 00084 throw Error(err.str()); 00085 //assert(false); 00086 } 00087 #endif // DROP_CGAL 00088 { 00089 ostringstream err; 00090 err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl; 00091 assert(false); 00092 throw Error(err.str()); 00093 } 00094 00095 // We will find nearest neighbour for each vertex, and include 00096 // distance in map (NB DistMap is a typedef given in the .h file) 00097 DistMap DijMap; 00098 00099 // fill the map with the minimal (as far as we know) subset of Dij 00100 // distances (i.e. nearest neighbour ones). 00101 for (int ii = 0; ii < n; ii++) { 00102 _add_ktdistance_to_map(ii, DijMap, DNN.get()); 00103 } 00104 00105 // run the clustering (go up to i=n-1, but then will stop half-way down, 00106 // when we reach that point -- it will be the final beam jet and there 00107 // will be no nearest neighbours to find). 00108 for (int i=0;i<n;i++) { 00109 // find nearest vertices 00110 // NB: skip cases where the point is not there anymore! 00111 TwoVertices SmallestDijPair; 00112 int jet_i, jet_j; 00113 double SmallestDij; 00114 bool Valid2; 00115 bool recombine_with_beam; 00116 do { 00117 SmallestDij = DijMap.begin()->first; 00118 SmallestDijPair = DijMap.begin()->second; 00119 jet_i = SmallestDijPair.first; 00120 jet_j = SmallestDijPair.second; 00121 // distance is immediately removed regardless of whether or not 00122 // it is used. 00123 // Some temporary testing code relating to problems with the gcc-3.2 compiler 00124 //cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n"; 00125 //cout << jet_i << " "<< jet_j<<"\n"; 00126 DijMap.erase(DijMap.begin()); 00127 //cout << "got beyond here\n"; 00128 00129 // need to "prime" the validity of jet_j in such a way that 00130 // if it corresponds to the beam then it is automatically valid. 00131 recombine_with_beam = (jet_j == BeamJet); 00132 if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);} 00133 else {Valid2 = true;} 00134 00135 } while ( !DNN->Valid(jet_i) || !Valid2); 00136 00137 00138 // The following part acts just on jet momenta and on the history. 00139 // The action on the nearest-neighbour structures takes place 00140 // later (only if at least 2 jets are around). 00141 if (! recombine_with_beam) { 00142 int nn; // will be index of new jet 00143 _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn); 00144 //OBS // merge the two jets, add new jet, remove old ones 00145 //OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]); 00146 //OBS 00147 //OBS int nn = _jets.size()-1; 00148 //OBS _jets[nn].set_cluster_hist_index(n+i); 00149 //OBS 00150 //OBS // get corresponding indices in history structure 00151 //OBS int hist_i = _jets[jet_i].cluster_hist_index(); 00152 //OBS int hist_j = _jets[jet_j].cluster_hist_index(); 00153 //OBS 00154 //OBS 00155 //OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j), 00156 //OBS _jets.size()-1, SmallestDij); 00157 00158 // add new point to points vector 00159 EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi()); 00160 newpoint.sanitize(); // make sure it is in correct range 00161 points.push_back(newpoint); 00162 } else { 00163 // recombine the jet with the beam 00164 _do_iB_recombination_step(jet_i, SmallestDij); 00165 //OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet, 00166 //OBS Invalid, SmallestDij); 00167 } 00168 00169 // exit the loop because we do not want to look for nearest neighbours 00170 // etc. of zero partons 00171 if (i == n-1) {break;} 00172 00173 vector<int> updated_neighbours; 00174 if (! recombine_with_beam) { 00175 // update DNN 00176 int point3; 00177 DNN->RemoveCombinedAddCombination(jet_i, jet_j, 00178 points[points.size()-1], point3, 00179 updated_neighbours); 00180 // C++ beginners' comment: static_cast to unsigned int is necessary 00181 // to do away with warnings about type mismatch between point3 (int) 00182 // and points.size (unsigned int) 00183 if (static_cast<unsigned int> (point3) != points.size()-1) { 00184 throw Error("INTERNAL ERROR: point3 != points.size()-1");} 00185 } else { 00186 // update DNN 00187 DNN->RemovePoint(jet_i, updated_neighbours); 00188 } 00189 00190 // update map 00191 vector<int>::iterator it = updated_neighbours.begin(); 00192 for (; it != updated_neighbours.end(); ++it) { 00193 int ii = *it; 00194 _add_ktdistance_to_map(ii, DijMap, DNN.get()); 00195 } 00196 00197 } // end clustering loop 00198 00199 } 00200 00201 00202 //---------------------------------------------------------------------- 00203 /// Add the current kt distance for particle ii to the map (DijMap) 00204 /// using information from the DNN object. Work as follows: 00205 /// 00206 /// . if the kt is zero then it's nearest neighbour is taken to be the 00207 /// the beam jet and the distance is zero. 00208 /// 00209 /// . if cylinder distance to nearest neighbour > _Rparam then it is 00210 /// yiB that is smallest and this is added to map. 00211 /// 00212 /// . otherwise if the nearest neighbour jj has a larger kt then add 00213 /// dij to the map. 00214 /// 00215 /// . otherwise do nothing 00216 /// 00217 void ClusterSequence::_add_ktdistance_to_map( 00218 const int & ii, 00219 DistMap & DijMap, 00220 const DynamicNearestNeighbours * DNN) { 00221 00222 double yiB = jet_scale_for_algorithm(_jets[ii]); 00223 if (yiB == 0.0) { 00224 // in this case convention is that we do not worry about distances 00225 // but directly state that nearest neighbour is beam 00226 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1))); 00227 } else { 00228 double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2; 00229 // Logic of following bit is: only add point to map if it has 00230 // smaller kt2 than nearest neighbour j (if it has larger kt, 00231 // then: either it is j's nearest neighbour and then we will 00232 // include dij when we come to j; or it is not j's nearest 00233 // neighbour and j will recombine with someone else). 00234 00235 // If DeltaR2 > 1.0 then in any case it will recombine with beam rather 00236 // than with any neighbours. 00237 // (put general normalisation here at some point) 00238 if (DeltaR2 > 1.0) { 00239 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1))); 00240 } else { 00241 double kt2i = jet_scale_for_algorithm(_jets[ii]); 00242 int jj = DNN->NearestNeighbourIndex(ii); 00243 if (kt2i <= jet_scale_for_algorithm(_jets[jj])) { 00244 double dij = DeltaR2 * kt2i; 00245 DijMap.insert(DijEntry(dij, TwoVertices(ii,jj))); 00246 } 00247 } 00248 } 00249 } 00250 00251 00252 FASTJET_END_NAMESPACE 00253