FastJet 3.0.2
|
00001 #ifndef __NNH_HH__ 00002 #define __NNH_HH__ 00003 00004 //STARTHEADER 00005 // $Id: NNH.hh 2696 2011-11-15 09:42:59Z soyez $ 00006 // 00007 // Copyright (c) 2009, Matteo Cacciari, Gavin Salam and Gregory Soyez 00008 // 00009 //---------------------------------------------------------------------- 00010 // This file is part of FastJet. 00011 // 00012 // FastJet is free software; you can redistribute it and/or modify 00013 // it under the terms of the GNU General Public License as published by 00014 // the Free Software Foundation; either version 2 of the License, or 00015 // (at your option) any later version. 00016 // 00017 // The algorithms that underlie FastJet have required considerable 00018 // development and are described in hep-ph/0512210. If you use 00019 // FastJet as part of work towards a scientific publication, please 00020 // include a citation to the FastJet paper. 00021 // 00022 // FastJet is distributed in the hope that it will be useful, 00023 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00024 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00025 // GNU General Public License for more details. 00026 // 00027 // You should have received a copy of the GNU General Public License 00028 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00029 //---------------------------------------------------------------------- 00030 //ENDHEADER 00031 00032 #include<fastjet/ClusterSequence.hh> 00033 00034 00035 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00036 00037 /// @ingroup advanced_usage 00038 /// \class _NoInfo 00039 /// dummy class, used as a default template argument 00040 class _NoInfo {}; 00041 00042 /// @ingroup advanced_usage 00043 /// \class NNHInfo 00044 /// template that will help initialise a BJ with a PseudoJet and extra information 00045 template<class I> class NNHInfo { 00046 public: 00047 NNHInfo() : _info(NULL) {} 00048 NNHInfo(I * info) : _info(info) {} 00049 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);} 00050 private: 00051 I * _info; 00052 }; 00053 00054 /// @ingroup advanced_usage 00055 /// Specialisation of NNHInfo for cases where there is no extra info 00056 template<> class NNHInfo<_NoInfo> { 00057 public: 00058 NNHInfo() {} 00059 NNHInfo(_NoInfo * ) {} 00060 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index);} 00061 }; 00062 00063 00064 //---------------------------------------------------------------------- 00065 /// @ingroup advanced_usage 00066 /// \class NNH 00067 /// Help solve closest pair problems with generic interparticle and 00068 /// beam distance. 00069 /// 00070 /// Class to help solve closest pair problems with generic interparticle 00071 /// distances and a beam distance, using Anderberg's Nearest Neighbour 00072 /// Heuristic. 00073 /// 00074 /// It is templated with a BJ (brief jet) class --- BJ should 00075 /// basically cache the minimal amount of information that is needed 00076 /// to efficiently calculate interparticle distances and particle-beam 00077 /// distances. 00078 /// 00079 /// This class can be used with or without an extra "Information" template, 00080 /// i.e. NNB<BJ> or NNH<BJ,I> 00081 /// 00082 /// For the NNH<BJ> version of the class to function, BJ must provide 00083 /// three member functions 00084 /// 00085 /// - void BJ::init(const PseudoJet & jet); // initialise with a PseudoJet 00086 /// - double BJ::distance(const BJ * other_bj_jet); // distance between this and other_bj_jet 00087 /// - double BJ::beam_distance() ; // distance to the beam 00088 /// 00089 /// For the NNH<BJ,I> version to function, the BJ::init(...) member 00090 /// must accept an extra argument 00091 /// 00092 /// - void BJ::init(const PseudoJet & jet, I * info); // initialise with a PseudoJet + info 00093 /// 00094 /// where info might be a pointer to a class that contains, e.g., information 00095 /// about R, or other parameters of the jet algorithm 00096 /// 00097 /// For an example of how the NNH<BJ> class is used, see the Jade (and 00098 /// EECambridge) plugins 00099 /// 00100 /// NB: the NNH algorithm is expected N^2, but has a worst case of 00101 /// N^3. Many QCD problems tend to place one closer to the N^3 end of 00102 /// the spectrum than one would like. There is scope for further 00103 /// progress (cf Eppstein, Cardinal & Eppstein), nevertheless the 00104 /// current class is already significantly faster than standard N^3 00105 /// implementations. 00106 /// 00107 /// 00108 /// Implementation note: this class derives from NNHInfo, which deals 00109 /// with storing any global information that is needed during the clustering 00110 00111 template<class BJ, class I = _NoInfo> class NNH : public NNHInfo<I> { 00112 public: 00113 00114 /// constructor with an initial set of jets (which will be assigned indices 00115 /// 0 ... jets.size()-1 00116 NNH(const std::vector<PseudoJet> & jets) {start(jets);} 00117 NNH(const std::vector<PseudoJet> & jets, I * info) : NNHInfo<I>(info) {start(jets);} 00118 00119 void start(const std::vector<PseudoJet> & jets); 00120 00121 /// return the dij_min and indices iA, iB, for the corresponding jets. 00122 /// If iB < 0 then iA recombines with the beam 00123 double dij_min(int & iA, int & iB); 00124 00125 /// remove the jet pointed to by index iA 00126 void remove_jet(int iA); 00127 00128 /// merge the jets pointed to by indices A and B and replace them with 00129 /// jet, assigning it an index jet_index. 00130 void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index); 00131 00132 /// a destructor 00133 ~NNH() { 00134 delete[] briefjets; 00135 } 00136 00137 private: 00138 class NNBJ; // forward declaration 00139 00140 /// establish the nearest neighbour for jet, and cross check constistency 00141 /// of distances for the other jets that are encountered. Assumes 00142 /// jet not contained within begin...end 00143 void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end); 00144 00145 /// establish the nearest neighbour for jet; don't cross check other jets' 00146 /// distances and allow jet to be contained within begin...end 00147 void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end); 00148 00149 /// contains the briefjets 00150 NNBJ * briefjets; 00151 00152 /// semaphores for the current extent of our structure 00153 NNBJ * head, * tail; 00154 00155 /// currently active number of jets 00156 int n; 00157 00158 /// where_is[i] contains a pointer to the jet with index i 00159 std::vector<NNBJ *> where_is; 00160 00161 /// a class that wraps around the BJ, supplementing it with extra information 00162 /// such as pointers to neighbours, etc. 00163 class NNBJ : public BJ { 00164 public: 00165 void init(const PseudoJet & jet, int index_in) { 00166 BJ::init(jet); 00167 other_init(index_in); 00168 } 00169 void init(const PseudoJet & jet, int index_in, I * info) { 00170 BJ::init(jet, info); 00171 other_init(index_in); 00172 } 00173 void other_init(int index_in) { 00174 _index = index_in; 00175 NN_dist = BJ::beam_distance(); 00176 NN = NULL; 00177 } 00178 int index() const {return _index;} 00179 00180 double NN_dist; 00181 NNBJ * NN; 00182 00183 private: 00184 int _index; 00185 }; 00186 00187 }; 00188 00189 00190 00191 //---------------------------------------------------------------------- 00192 template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) { 00193 n = jets.size(); 00194 briefjets = new NNBJ[n]; 00195 where_is.resize(2*n); 00196 00197 NNBJ * jetA = briefjets; 00198 00199 // initialise the basic jet info 00200 for (int i = 0; i< n; i++) { 00201 //jetA->init(jets[i], i); 00202 init_jet(jetA, jets[i], i); 00203 where_is[i] = jetA; 00204 jetA++; // move on to next entry of briefjets 00205 } 00206 tail = jetA; // a semaphore for the end of briefjets 00207 head = briefjets; // a nicer way of naming start 00208 00209 // now initialise the NN distances: jetA will run from 1..n-1; and 00210 // jetB from 0..jetA-1 00211 for (jetA = head + 1; jetA != tail; jetA++) { 00212 // set NN info for jetA based on jets running from head..jetA-1, 00213 // checking in the process whether jetA itself is an undiscovered 00214 // NN of one of those jets. 00215 set_NN_crosscheck(jetA, head, jetA); 00216 } 00217 //std::cout << "OOOO " << briefjets[1].NN_dist << " " << briefjets[1].NN - briefjets << std::endl; 00218 } 00219 00220 00221 //---------------------------------------------------------------------- 00222 template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) { 00223 // find the minimum of the diJ on this round 00224 double diJ_min = briefjets[0].NN_dist; 00225 int diJ_min_jet = 0; 00226 for (int i = 1; i < n; i++) { 00227 if (briefjets[i].NN_dist < diJ_min) { 00228 diJ_min_jet = i; 00229 diJ_min = briefjets[i].NN_dist; 00230 } 00231 } 00232 00233 // return information to user about recombination 00234 NNBJ * jetA = & briefjets[diJ_min_jet]; 00235 //std::cout << jetA - briefjets << std::endl; 00236 iA = jetA->index(); 00237 iB = jetA->NN ? jetA->NN->index() : -1; 00238 return diJ_min; 00239 } 00240 00241 00242 //---------------------------------------------------------------------- 00243 // remove jetA from the list 00244 template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) { 00245 NNBJ * jetA = where_is[iA]; 00246 // now update our nearest neighbour info and diJ table 00247 // first reduce size of table 00248 tail--; n--; 00249 // Copy last jet contents and diJ info into position of jetA 00250 *jetA = *tail; 00251 // update the info on where the given index is stored 00252 where_is[jetA->index()] = jetA; 00253 00254 for (NNBJ * jetI = head; jetI != tail; jetI++) { 00255 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN 00256 if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail); 00257 00258 // if jetI's NN is the new tail then relabel it so that it becomes jetA 00259 if (jetI->NN == tail) {jetI->NN = jetA;} 00260 } 00261 } 00262 00263 00264 //---------------------------------------------------------------------- 00265 template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB, 00266 const PseudoJet & jet, int index) { 00267 00268 NNBJ * jetA = where_is[iA]; 00269 NNBJ * jetB = where_is[iB]; 00270 00271 // If necessary relabel A & B to ensure jetB < jetA, that way if 00272 // the larger of them == newtail then that ends up being jetA and 00273 // the new jet that is added as jetB is inserted in a position that 00274 // has a future! 00275 if (jetA < jetB) std::swap(jetA,jetB); 00276 00277 // initialise jetB based on the new jet 00278 //jetB->init(jet, index); 00279 init_jet(jetB, jet, index); 00280 // and record its position (making sure we have the space) 00281 if (index >= int(where_is.size())) where_is.resize(2*index); 00282 where_is[jetB->index()] = jetB; 00283 00284 // now update our nearest neighbour info 00285 // first reduce size of table 00286 tail--; n--; 00287 // Copy last jet contents into position of jetA 00288 *jetA = *tail; 00289 // update the info on where the tail's index is stored 00290 where_is[jetA->index()] = jetA; 00291 00292 for (NNBJ * jetI = head; jetI != tail; jetI++) { 00293 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN 00294 if (jetI->NN == jetA || jetI->NN == jetB) { 00295 set_NN_nocross(jetI, head, tail); 00296 } 00297 00298 // check whether new jetB is closer than jetI's current NN and 00299 // if need be update things 00300 double dist = jetI->distance(jetB); 00301 if (dist < jetI->NN_dist) { 00302 if (jetI != jetB) { 00303 jetI->NN_dist = dist; 00304 jetI->NN = jetB; 00305 } 00306 } 00307 if (dist < jetB->NN_dist) { 00308 if (jetI != jetB) { 00309 jetB->NN_dist = dist; 00310 jetB->NN = jetI; 00311 } 00312 } 00313 00314 // if jetI's NN is the new tail then relabel it so that it becomes jetA 00315 if (jetI->NN == tail) {jetI->NN = jetA;} 00316 } 00317 } 00318 00319 00320 //---------------------------------------------------------------------- 00321 // this function assumes that jet is not contained within begin...end 00322 template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet, 00323 NNBJ * begin, NNBJ * end) { 00324 double NN_dist = jet->beam_distance(); 00325 NNBJ * NN = NULL; 00326 for (NNBJ * jetB = begin; jetB != end; jetB++) { 00327 double dist = jet->distance(jetB); 00328 if (dist < NN_dist) { 00329 NN_dist = dist; 00330 NN = jetB; 00331 } 00332 if (dist < jetB->NN_dist) { 00333 jetB->NN_dist = dist; 00334 jetB->NN = jet; 00335 } 00336 } 00337 jet->NN = NN; 00338 jet->NN_dist = NN_dist; 00339 } 00340 00341 00342 //---------------------------------------------------------------------- 00343 // set the NN for jet without checking whether in the process you might 00344 // have discovered a new nearest neighbour for another jet 00345 template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross( 00346 NNBJ * jet, NNBJ * begin, NNBJ * end) { 00347 double NN_dist = jet->beam_distance(); 00348 NNBJ * NN = NULL; 00349 // if (head < jet) { 00350 // for (NNBJ * jetB = head; jetB != jet; jetB++) { 00351 if (begin < jet) { 00352 for (NNBJ * jetB = begin; jetB != jet; jetB++) { 00353 double dist = jet->distance(jetB); 00354 if (dist < NN_dist) { 00355 NN_dist = dist; 00356 NN = jetB; 00357 } 00358 } 00359 } 00360 // if (tail > jet) { 00361 // for (NNBJ * jetB = jet+1; jetB != tail; jetB++) { 00362 if (end > jet) { 00363 for (NNBJ * jetB = jet+1; jetB != end; jetB++) { 00364 double dist = jet->distance (jetB); 00365 if (dist < NN_dist) { 00366 NN_dist = dist; 00367 NN = jetB; 00368 } 00369 } 00370 } 00371 jet->NN = NN; 00372 jet->NN_dist = NN_dist; 00373 } 00374 00375 00376 00377 00378 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh 00379 00380 00381 #endif // __NNH_HH__