FastJet 3.0.2
NNH.hh
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__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends