FastJet 3.0.2
ClusterSequence_TiledN2.cc
00001 //STARTHEADER
00002 // $Id: ClusterSequence_TiledN2.cc 2696 2011-11-15 09:42:59Z soyez $
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 // The plain N^2 part of the ClusterSequence class -- separated out
00031 // from the rest of the class implementation so as to speed up
00032 // compilation of this particular part while it is under test.
00033 
00034 #include<iostream>
00035 #include<vector>
00036 #include<cmath>
00037 #include<algorithm>
00038 #include "fastjet/PseudoJet.hh"
00039 #include "fastjet/ClusterSequence.hh"
00040 #include "fastjet/internal/MinHeap.hh"
00041 
00042 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00043 
00044 using namespace std;
00045 
00046 
00047 //----------------------------------------------------------------------
00048 void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) {
00049   Tile * tile = & _tiles[jet->tile_index];
00050 
00051   if (jet->previous == NULL) {
00052     // we are at head of the tile, so reset it.
00053     // If this was the only jet on the tile then tile->head will now be NULL
00054     tile->head = jet->next;
00055   } else {
00056     // adjust link from previous jet in this tile
00057     jet->previous->next = jet->next;
00058   }
00059   if (jet->next != NULL) {
00060     // adjust backwards-link from next jet in this tile
00061     jet->next->previous = jet->previous;
00062   }
00063 }
00064 
00065 //----------------------------------------------------------------------
00066 /// Set up the tiles:
00067 ///  - decide the range in eta
00068 ///  - allocate the tiles
00069 ///  - set up the cross-referencing info between tiles
00070 ///
00071 /// The neighbourhood of a tile is set up as follows
00072 ///
00073 ///           LRR
00074 ///           LXR
00075 ///           LLR
00076 ///
00077 /// such that tiles is an array containing XLLLLRRRR with pointers
00078 ///                                         |   \ RH_tiles
00079 ///                                         \ surrounding_tiles
00080 ///
00081 /// with appropriate precautions when close to the edge of the tiled
00082 /// region.
00083 ///
00084 void ClusterSequence::_initialise_tiles() {
00085 
00086   // first decide tile sizes (with a lower bound to avoid huge memory use with
00087   // very small R)
00088   double default_size = max(0.1,_Rparam);
00089   _tile_size_eta = default_size;
00090   // it makes no sense to go below 3 tiles in phi -- 3 tiles is
00091   // sufficient to make sure all pair-wise combinations up to pi in
00092   // phi are possible
00093   _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
00094   _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
00095 
00096   // always include zero rapidity in the tiling region
00097   _tiles_eta_min = 0.0;
00098   _tiles_eta_max = 0.0;
00099   // but go no further than following
00100   const double maxrap = 7.0;
00101 
00102   // and find out how much further one should go
00103   for(unsigned int i = 0; i < _jets.size(); i++) {
00104     double eta = _jets[i].rap();
00105     // first check if eta is in range -- to avoid taking into account
00106     // very spurious rapidities due to particles with near-zero kt.
00107     if (abs(eta) < maxrap) {
00108       if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
00109       if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
00110     }
00111   }
00112 
00113   // now adjust the values
00114   _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
00115   _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
00116   _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
00117   _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
00118 
00119   // allocate the tiles
00120   _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
00121 
00122   // now set up the cross-referencing between tiles
00123   for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
00124     for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
00125       Tile * tile = & _tiles[_tile_index(ieta,iphi)];
00126       // no jets in this tile yet
00127       tile->head = NULL; // first element of tiles points to itself
00128       tile->begin_tiles[0] =  tile;
00129       Tile ** pptile = & (tile->begin_tiles[0]);
00130       pptile++;
00131       //
00132       // set up L's in column to the left of X
00133       tile->surrounding_tiles = pptile;
00134       if (ieta > _tiles_ieta_min) {
00135         // with the itile subroutine, we can safely run tiles from
00136         // idphi=-1 to idphi=+1, because it takes care of
00137         // negative and positive boundaries
00138         for (int idphi = -1; idphi <=+1; idphi++) {
00139           *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
00140           pptile++;
00141         }       
00142       }
00143       // now set up last L (below X)
00144       *pptile = & _tiles[_tile_index(ieta,iphi-1)];
00145       pptile++;
00146       // set up first R (above X)
00147       tile->RH_tiles = pptile;
00148       *pptile = & _tiles[_tile_index(ieta,iphi+1)];
00149       pptile++;
00150       // set up remaining R's, to the right of X
00151       if (ieta < _tiles_ieta_max) {
00152         for (int idphi = -1; idphi <= +1; idphi++) {
00153           *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
00154           pptile++;
00155         }       
00156       }
00157       // now put semaphore for end tile
00158       tile->end_tiles = pptile;
00159       // finally make sure tiles are untagged
00160       tile->tagged = false;
00161     }
00162   }
00163 
00164 }
00165 
00166 
00167 //----------------------------------------------------------------------
00168 /// return the tile index corresponding to the given eta,phi point
00169 int ClusterSequence::_tile_index(const double & eta, const double & phi) const {
00170   int ieta, iphi;
00171   if      (eta <= _tiles_eta_min) {ieta = 0;}
00172   else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
00173   else {
00174     //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
00175     ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
00176     // following needed in case of rare but nasty rounding errors
00177     if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
00178       ieta = _tiles_ieta_max-_tiles_ieta_min;} 
00179   }
00180   // allow for some extent of being beyond range in calculation of phi
00181   // as well
00182   //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
00183   // with just int and no floor, things run faster but beware
00184   iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
00185   return (iphi + ieta * _n_tiles_phi);
00186 }
00187 
00188 
00189 //----------------------------------------------------------------------
00190 // overloaded version which additionally sets up information regarding the
00191 // tiling
00192 inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
00193                                               const int _jets_index) {
00194   // first call the generic setup
00195   _bj_set_jetinfo<>(jet, _jets_index);
00196 
00197   // Then do the setup specific to the tiled case.
00198 
00199   // Find out which tile it belonds to
00200   jet->tile_index = _tile_index(jet->eta, jet->phi);
00201 
00202   // Insert it into the tile's linked list of jets
00203   Tile * tile = &_tiles[jet->tile_index];
00204   jet->previous   = NULL;
00205   jet->next       = tile->head;
00206   if (jet->next != NULL) {jet->next->previous = jet;}
00207   tile->head      = jet;
00208 }
00209 
00210 
00211 //----------------------------------------------------------------------
00212 /// output the contents of the tiles
00213 void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
00214   for (vector<Tile>::const_iterator tile = _tiles.begin(); 
00215        tile < _tiles.end(); tile++) {
00216     cout << "Tile " << tile - _tiles.begin()<<" = ";
00217     vector<int> list;
00218     for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00219       list.push_back(jetI-briefjets);
00220       //cout <<" "<<jetI-briefjets;
00221     }
00222     sort(list.begin(),list.end());
00223     for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
00224     cout <<"\n";
00225   }
00226 }
00227 
00228 
00229 //----------------------------------------------------------------------
00230 /// Add to the vector tile_union the tiles that are in the neighbourhood
00231 /// of the specified tile_index, including itself -- start adding
00232 /// from position n_near_tiles-1, and increase n_near_tiles as
00233 /// you go along (could have done it more C++ like with vector with reserved
00234 /// space, but fear is that it would have been slower, e.g. checking
00235 /// for end of vector at each stage to decide whether to resize it)
00236 void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index, 
00237                vector<int> & tile_union, int & n_near_tiles) const {
00238   for (Tile * const * near_tile = _tiles[tile_index].begin_tiles; 
00239        near_tile != _tiles[tile_index].end_tiles; near_tile++){
00240     // get the tile number
00241     tile_union[n_near_tiles] = *near_tile - & _tiles[0];
00242     n_near_tiles++;
00243   }
00244 }
00245 
00246 
00247 //----------------------------------------------------------------------
00248 /// Like _add_neighbours_to_tile_union, but only adds neighbours if 
00249 /// their "tagged" status is false; when a neighbour is added its
00250 /// tagged status is set to true.
00251 inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
00252                const int tile_index, 
00253                vector<int> & tile_union, int & n_near_tiles)  {
00254   for (Tile ** near_tile = _tiles[tile_index].begin_tiles; 
00255        near_tile != _tiles[tile_index].end_tiles; near_tile++){
00256     if (! (*near_tile)->tagged) {
00257       (*near_tile)->tagged = true;
00258       // get the tile number
00259       tile_union[n_near_tiles] = *near_tile - & _tiles[0];
00260       n_near_tiles++;
00261     }
00262   }
00263 }
00264 
00265 
00266 //----------------------------------------------------------------------
00267 /// run a tiled clustering
00268 void ClusterSequence::_tiled_N2_cluster() {
00269 
00270   _initialise_tiles();
00271 
00272   int n = _jets.size();
00273   TiledJet * briefjets = new TiledJet[n];
00274   TiledJet * jetA = briefjets, * jetB;
00275   TiledJet oldB;
00276   
00277 
00278   // will be used quite deep inside loops, but declare it here so that
00279   // memory (de)allocation gets done only once
00280   vector<int> tile_union(3*n_tile_neighbours);
00281   
00282   // initialise the basic jet info 
00283   for (int i = 0; i< n; i++) {
00284     _tj_set_jetinfo(jetA, i);
00285     //cout << i<<": "<<jetA->tile_index<<"\n";
00286     jetA++; // move on to next entry of briefjets
00287   }
00288   TiledJet * tail = jetA; // a semaphore for the end of briefjets
00289   TiledJet * head = briefjets; // a nicer way of naming start
00290 
00291   // set up the initial nearest neighbour information
00292   vector<Tile>::const_iterator tile;
00293   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00294     // first do it on this tile
00295     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00296       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00297         double dist = _bj_dist(jetA,jetB);
00298         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00299         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00300       }
00301     }
00302     // then do it for RH tiles
00303     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00304       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00305         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00306           double dist = _bj_dist(jetA,jetB);
00307           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00308           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00309         }
00310       }
00311     }
00312   }
00313   
00314   // now create the diJ (where J is i's NN) table -- remember that 
00315   // we differ from standard normalisation here by a factor of R2
00316   double * diJ = new double[n];
00317   jetA = head;
00318   for (int i = 0; i < n; i++) {
00319     diJ[i] = _bj_diJ(jetA);
00320     jetA++; // have jetA follow i
00321   }
00322 
00323   // now run the recombination loop
00324   int history_location = n-1;
00325   while (tail != head) {
00326 
00327     // find the minimum of the diJ on this round
00328     double diJ_min = diJ[0];
00329     int diJ_min_jet = 0;
00330     for (int i = 1; i < n; i++) {
00331       if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
00332     }
00333 
00334     // do the recombination between A and B
00335     history_location++;
00336     jetA = & briefjets[diJ_min_jet];
00337     jetB = jetA->NN;
00338     // put the normalisation back in
00339     diJ_min *= _invR2; 
00340 
00341     //if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";}
00342 
00343     //cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n";
00344 
00345     if (jetB != NULL) {
00346       // jet-jet recombination
00347       // If necessary relabel A & B to ensure jetB < jetA, that way if
00348       // the larger of them == newtail then that ends up being jetA and 
00349       // the new jet that is added as jetB is inserted in a position that
00350       // has a future!
00351       if (jetA < jetB) {std::swap(jetA,jetB);}
00352 
00353       int nn; // new jet index
00354       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00355 
00356       //OBS// get the two history indices
00357       //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00358       //OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index();
00359       //OBS// create the recombined jet
00360       //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00361       //OBSint nn = _jets.size() - 1;
00362       //OBS_jets[nn].set_cluster_hist_index(history_location);
00363       //OBS// update history
00364       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
00365       //OBS_add_step_to_history(history_location, 
00366       //OBS                min(hist_a,hist_b),max(hist_a,hist_b),
00367       //OBS                        nn, diJ_min);
00368 
00369       // what was jetB will now become the new jet
00370       _bj_remove_from_tiles(jetA);
00371       oldB = * jetB;  // take a copy because we will need it...
00372       _bj_remove_from_tiles(jetB);
00373       _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
00374     } else {
00375       // jet-beam recombination
00376       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00377             
00378       //OBS// get the hist_index
00379       //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00380       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
00381       //OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); 
00382       _bj_remove_from_tiles(jetA);
00383     }
00384 
00385     // first establish the set of tiles over which we are going to 
00386     // have to run searches for updated and new nearest-neighbours
00387     int n_near_tiles = 0;
00388     _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
00389     if (jetB != NULL) {
00390       bool sort_it = false;
00391       if (jetB->tile_index != jetA->tile_index) {
00392         sort_it = true;
00393         _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
00394       }
00395       if (oldB.tile_index != jetA->tile_index && 
00396           oldB.tile_index != jetB->tile_index) {
00397         sort_it = true;
00398         _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
00399       }
00400 
00401       if (sort_it) {
00402         // sort the tiles before then compressing the list
00403         sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
00404         // and now condense the list
00405         int nnn = 1;
00406         for (int i = 1; i < n_near_tiles; i++) {
00407           if (tile_union[i] != tile_union[nnn-1]) {
00408             tile_union[nnn] = tile_union[i]; 
00409             nnn++;
00410           }
00411         }
00412         n_near_tiles = nnn;
00413       }
00414     }
00415 
00416     // now update our nearest neighbour info and diJ table
00417     // first reduce size of table
00418     tail--; n--;
00419     if (jetA == tail) {
00420       // there is nothing to be done
00421     } else {
00422       // Copy last jet contents and diJ info into position of jetA
00423       *jetA = *tail;
00424       diJ[jetA - head] = diJ[tail-head];
00425       // IN the tiling fix pointers to tail and turn them into
00426       // pointers to jetA (from predecessors, successors and the tile
00427       // head if need be)
00428       if (jetA->previous == NULL) {
00429         _tiles[jetA->tile_index].head = jetA;
00430       } else {
00431         jetA->previous->next = jetA;
00432       }
00433       if (jetA->next != NULL) {jetA->next->previous = jetA;}
00434     }
00435 
00436     // Initialise jetB's NN distance as well as updating it for 
00437     // other particles.
00438     for (int itile = 0; itile < n_near_tiles; itile++) {
00439       Tile * tile_ptr = &_tiles[tile_union[itile]];
00440       for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
00441         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00442         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00443           jetI->NN_dist = _R2;
00444           jetI->NN      = NULL;
00445           // now go over tiles that are neighbours of I (include own tile)
00446           for (Tile ** near_tile  = tile_ptr->begin_tiles; 
00447                        near_tile != tile_ptr->end_tiles; near_tile++) {
00448             // and then over the contents of that tile
00449             for (TiledJet * jetJ  = (*near_tile)->head; 
00450                             jetJ != NULL; jetJ = jetJ->next) {
00451               double dist = _bj_dist(jetI,jetJ);
00452               if (dist < jetI->NN_dist && jetJ != jetI) {
00453                 jetI->NN_dist = dist; jetI->NN = jetJ;
00454               }
00455             }
00456           }
00457           diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
00458         }
00459         // check whether new jetB is closer than jetI's current NN and
00460         // if need to update things
00461         if (jetB != NULL) {
00462           double dist = _bj_dist(jetI,jetB);
00463           if (dist < jetI->NN_dist) {
00464             if (jetI != jetB) {
00465               jetI->NN_dist = dist;
00466               jetI->NN = jetB;
00467               diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
00468             }
00469           }
00470           if (dist < jetB->NN_dist) {
00471             if (jetI != jetB) {
00472               jetB->NN_dist = dist;
00473               jetB->NN      = jetI;}
00474           }
00475         }
00476       }
00477     }
00478 
00479 
00480     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00481     //cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
00482 
00483     // remember to update pointers to tail
00484     for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles; 
00485                  near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
00486       // and then the contents of that tile
00487       for (TiledJet * jetJ = (*near_tile)->head; 
00488                      jetJ != NULL; jetJ = jetJ->next) {
00489         if (jetJ->NN == tail) {jetJ->NN = jetA;}
00490       }
00491     }
00492 
00493     //for (int i = 0; i < n; i++) {
00494     //  if (briefjets[i].NN-briefjets >= n && briefjets[i].NN != NULL) {cout <<"YOU MUST BE CRAZY for n ="<<n<<", i = "<<i<<", NN = "<<briefjets[i].NN-briefjets<<"\n";}
00495     //}
00496 
00497 
00498     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00499     //cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
00500 
00501   }
00502 
00503   // final cleaning up;
00504   delete[] diJ;
00505   delete[] briefjets;
00506 }
00507 
00508 
00509 //----------------------------------------------------------------------
00510 /// run a tiled clustering
00511 void ClusterSequence::_faster_tiled_N2_cluster() {
00512 
00513   _initialise_tiles();
00514 
00515   int n = _jets.size();
00516   TiledJet * briefjets = new TiledJet[n];
00517   TiledJet * jetA = briefjets, * jetB;
00518   TiledJet oldB;
00519   
00520 
00521   // will be used quite deep inside loops, but declare it here so that
00522   // memory (de)allocation gets done only once
00523   vector<int> tile_union(3*n_tile_neighbours);
00524   
00525   // initialise the basic jet info 
00526   for (int i = 0; i< n; i++) {
00527     _tj_set_jetinfo(jetA, i);
00528     //cout << i<<": "<<jetA->tile_index<<"\n";
00529     jetA++; // move on to next entry of briefjets
00530   }
00531   TiledJet * head = briefjets; // a nicer way of naming start
00532 
00533   // set up the initial nearest neighbour information
00534   vector<Tile>::const_iterator tile;
00535   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00536     // first do it on this tile
00537     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00538       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00539         double dist = _bj_dist(jetA,jetB);
00540         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00541         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00542       }
00543     }
00544     // then do it for RH tiles
00545     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00546       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00547         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00548           double dist = _bj_dist(jetA,jetB);
00549           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00550           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00551         }
00552       }
00553     }
00554     // no need to do it for LH tiles, since they are implicitly done
00555     // when we set NN for both jetA and jetB on the RH tiles.
00556   }
00557 
00558   
00559   // now create the diJ (where J is i's NN) table -- remember that 
00560   // we differ from standard normalisation here by a factor of R2
00561   // (corrected for at the end). 
00562   struct diJ_plus_link {
00563     double     diJ; // the distance
00564     TiledJet * jet; // the jet (i) for which we've found this distance
00565                     // (whose NN will the J).
00566   };
00567   diJ_plus_link * diJ = new diJ_plus_link[n];
00568   jetA = head;
00569   for (int i = 0; i < n; i++) {
00570     diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
00571     diJ[i].jet = jetA;  // our compact diJ table will not be in      
00572     jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
00573                         // so set up bi-directional correspondence here.
00574     jetA++; // have jetA follow i 
00575   }
00576 
00577   // now run the recombination loop
00578   int history_location = n-1;
00579   while (n > 0) {
00580 
00581     // find the minimum of the diJ on this round
00582     diJ_plus_link * best, *stop; // pointers a bit faster than indices
00583     // could use best to keep track of diJ min, but it turns out to be
00584     // marginally faster to have a separate variable (avoids n
00585     // dereferences at the expense of n/2 assignments).
00586     double diJ_min = diJ[0].diJ; // initialise the best one here.
00587     best = diJ;                  // and here
00588     stop = diJ+n;
00589     for (diJ_plus_link * here = diJ+1; here != stop; here++) {
00590       if (here->diJ < diJ_min) {best = here; diJ_min  = here->diJ;}
00591     }
00592 
00593     // do the recombination between A and B
00594     history_location++;
00595     jetA = best->jet;
00596     jetB = jetA->NN;
00597     // put the normalisation back in
00598     diJ_min *= _invR2; 
00599 
00600     if (jetB != NULL) {
00601       // jet-jet recombination
00602       // If necessary relabel A & B to ensure jetB < jetA, that way if
00603       // the larger of them == newtail then that ends up being jetA and 
00604       // the new jet that is added as jetB is inserted in a position that
00605       // has a future!
00606       if (jetA < jetB) {std::swap(jetA,jetB);}
00607 
00608       int nn; // new jet index
00609       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00610       
00611       //OBS// get the two history indices
00612       //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
00613       //OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index();
00614       //OBS// create the recombined jet
00615       //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00616       //OBSint nn = _jets.size() - 1;
00617       //OBS_jets[nn].set_cluster_hist_index(history_location);
00618       //OBS// update history
00619       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
00620       //OBS_add_step_to_history(history_location, 
00621       //OBS                min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b),
00622       //OBS                        nn, diJ_min);
00623       // what was jetB will now become the new jet
00624       _bj_remove_from_tiles(jetA);
00625       oldB = * jetB;  // take a copy because we will need it...
00626       _bj_remove_from_tiles(jetB);
00627       _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
00628                                  // (also registers the jet in the tiling)
00629     } else {
00630       // jet-beam recombination
00631       // get the hist_index
00632       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00633       //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
00634       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
00635       //OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min); 
00636       _bj_remove_from_tiles(jetA);
00637     }
00638 
00639     // first establish the set of tiles over which we are going to
00640     // have to run searches for updated and new nearest-neighbours --
00641     // basically a combination of vicinity of the tiles of the two old
00642     // and one new jet.
00643     int n_near_tiles = 0;
00644     _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
00645                                            tile_union, n_near_tiles);
00646     if (jetB != NULL) {
00647       if (jetB->tile_index != jetA->tile_index) {
00648         _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00649                                                tile_union,n_near_tiles);
00650       }
00651       if (oldB.tile_index != jetA->tile_index && 
00652           oldB.tile_index != jetB->tile_index) {
00653         _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00654                                                tile_union,n_near_tiles);
00655       }
00656     }
00657 
00658     // now update our nearest neighbour info and diJ table
00659     // first reduce size of table
00660     n--;
00661     // then compactify the diJ by taking the last of the diJ and copying
00662     // it to the position occupied by the diJ for jetA
00663     diJ[n].jet->diJ_posn = jetA->diJ_posn;
00664     diJ[jetA->diJ_posn] = diJ[n];
00665 
00666     // Initialise jetB's NN distance as well as updating it for 
00667     // other particles.
00668     // Run over all tiles in our union 
00669     for (int itile = 0; itile < n_near_tiles; itile++) {
00670       Tile * tile_ptr = &_tiles[tile_union[itile]];
00671       tile_ptr->tagged = false; // reset tag, since we're done with unions
00672       // run over all jets in the current tile
00673       for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
00674         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00675         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00676           jetI->NN_dist = _R2;
00677           jetI->NN      = NULL;
00678           // now go over tiles that are neighbours of I (include own tile)
00679           for (Tile ** near_tile  = tile_ptr->begin_tiles; 
00680                        near_tile != tile_ptr->end_tiles; near_tile++) {
00681             // and then over the contents of that tile
00682             for (TiledJet * jetJ  = (*near_tile)->head; 
00683                             jetJ != NULL; jetJ = jetJ->next) {
00684               double dist = _bj_dist(jetI,jetJ);
00685               if (dist < jetI->NN_dist && jetJ != jetI) {
00686                 jetI->NN_dist = dist; jetI->NN = jetJ;
00687               }
00688             }
00689           }
00690           diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
00691         }
00692         // check whether new jetB is closer than jetI's current NN and
00693         // if jetI is closer than jetB's current (evolving) nearest
00694         // neighbour. Where relevant update things
00695         if (jetB != NULL) {
00696           double dist = _bj_dist(jetI,jetB);
00697           if (dist < jetI->NN_dist) {
00698             if (jetI != jetB) {
00699               jetI->NN_dist = dist;
00700               jetI->NN = jetB;
00701               diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
00702             }
00703           }
00704           if (dist < jetB->NN_dist) {
00705             if (jetI != jetB) {
00706               jetB->NN_dist = dist;
00707               jetB->NN      = jetI;}
00708           }
00709         }
00710       }
00711     }
00712 
00713     // finally, register the updated kt distance for B
00714     if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
00715 
00716   }
00717 
00718   // final cleaning up;
00719   delete[] diJ;
00720   delete[] briefjets;
00721 }
00722 
00723 
00724 
00725 //----------------------------------------------------------------------
00726 /// run a tiled clustering, with our minheap for keeping track of the
00727 /// smallest dij
00728 void ClusterSequence::_minheap_faster_tiled_N2_cluster() {
00729 
00730   _initialise_tiles();
00731 
00732   int n = _jets.size();
00733   TiledJet * briefjets = new TiledJet[n];
00734   TiledJet * jetA = briefjets, * jetB;
00735   TiledJet oldB;
00736   
00737 
00738   // will be used quite deep inside loops, but declare it here so that
00739   // memory (de)allocation gets done only once
00740   vector<int> tile_union(3*n_tile_neighbours);
00741   
00742   // initialise the basic jet info 
00743   for (int i = 0; i< n; i++) {
00744     _tj_set_jetinfo(jetA, i);
00745     //cout << i<<": "<<jetA->tile_index<<"\n";
00746     jetA++; // move on to next entry of briefjets
00747   }
00748   TiledJet * head = briefjets; // a nicer way of naming start
00749 
00750   // set up the initial nearest neighbour information
00751   vector<Tile>::const_iterator tile;
00752   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00753     // first do it on this tile
00754     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00755       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00756         double dist = _bj_dist(jetA,jetB);
00757         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00758         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00759       }
00760     }
00761     // then do it for RH tiles
00762     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00763       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00764         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00765           double dist = _bj_dist(jetA,jetB);
00766           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00767           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00768         }
00769       }
00770     }
00771     // no need to do it for LH tiles, since they are implicitly done
00772     // when we set NN for both jetA and jetB on the RH tiles.
00773   }
00774 
00775   
00776   //// now create the diJ (where J is i's NN) table -- remember that 
00777   //// we differ from standard normalisation here by a factor of R2
00778   //// (corrected for at the end). 
00779   //struct diJ_plus_link {
00780   //  double     diJ; // the distance
00781   //  TiledJet * jet; // the jet (i) for which we've found this distance
00782   //                  // (whose NN will the J).
00783   //};
00784   //diJ_plus_link * diJ = new diJ_plus_link[n];
00785   //jetA = head;
00786   //for (int i = 0; i < n; i++) {
00787   //  diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
00788   //  diJ[i].jet = jetA;  // our compact diJ table will not be in            
00789   //  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
00790   //                      // so set up bi-directional correspondence here.
00791   //  jetA++; // have jetA follow i 
00792   //}
00793 
00794   vector<double> diJs(n);
00795   for (int i = 0; i < n; i++) {
00796     diJs[i] = _bj_diJ(&briefjets[i]);
00797     briefjets[i].label_minheap_update_done();
00798   }
00799   MinHeap minheap(diJs);
00800   // have a stack telling us which jets we'll have to update on the heap
00801   vector<TiledJet *> jets_for_minheap;
00802   jets_for_minheap.reserve(n); 
00803 
00804   // now run the recombination loop
00805   int history_location = n-1;
00806   while (n > 0) {
00807 
00808     double diJ_min = minheap.minval() *_invR2;
00809     jetA = head + minheap.minloc();
00810 
00811     // do the recombination between A and B
00812     history_location++;
00813     jetB = jetA->NN;
00814 
00815     if (jetB != NULL) {
00816       // jet-jet recombination
00817       // If necessary relabel A & B to ensure jetB < jetA, that way if
00818       // the larger of them == newtail then that ends up being jetA and 
00819       // the new jet that is added as jetB is inserted in a position that
00820       // has a future!
00821       if (jetA < jetB) {std::swap(jetA,jetB);}
00822 
00823       int nn; // new jet index
00824       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00825       
00826       // what was jetB will now become the new jet
00827       _bj_remove_from_tiles(jetA);
00828       oldB = * jetB;  // take a copy because we will need it...
00829       _bj_remove_from_tiles(jetB);
00830       _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
00831                                  // (also registers the jet in the tiling)
00832     } else {
00833       // jet-beam recombination
00834       // get the hist_index
00835       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00836       _bj_remove_from_tiles(jetA);
00837     }
00838 
00839     // remove the minheap entry for jetA
00840     minheap.remove(jetA-head);
00841 
00842     // first establish the set of tiles over which we are going to
00843     // have to run searches for updated and new nearest-neighbours --
00844     // basically a combination of vicinity of the tiles of the two old
00845     // and one new jet.
00846     int n_near_tiles = 0;
00847     _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
00848                                            tile_union, n_near_tiles);
00849     if (jetB != NULL) {
00850       if (jetB->tile_index != jetA->tile_index) {
00851         _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00852                                                tile_union,n_near_tiles);
00853       }
00854       if (oldB.tile_index != jetA->tile_index && 
00855           oldB.tile_index != jetB->tile_index) {
00856         // GS: the line below generates a warning that oldB.tile_index
00857         // may be used uninitialised. However, to reach this point, we
00858         // ned jetB != NULL (see test a few lines above) and is jetB
00859         // !=NULL, one would have gone through "oldB = *jetB before
00860         // (see piece of code ~20 line above), so the index is
00861         // initialised. We do not do anything to avoid the warning to
00862         // avoid any potential speed impact.
00863         _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00864                                                tile_union,n_near_tiles);
00865       }
00866       // indicate that we'll have to update jetB in the minheap
00867       jetB->label_minheap_update_needed();
00868       jets_for_minheap.push_back(jetB);
00869     }
00870 
00871 
00872     // Initialise jetB's NN distance as well as updating it for 
00873     // other particles.
00874     // Run over all tiles in our union 
00875     for (int itile = 0; itile < n_near_tiles; itile++) {
00876       Tile * tile_ptr = &_tiles[tile_union[itile]];
00877       tile_ptr->tagged = false; // reset tag, since we're done with unions
00878       // run over all jets in the current tile
00879       for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
00880         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00881         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00882           jetI->NN_dist = _R2;
00883           jetI->NN      = NULL;
00884           // label jetI as needing heap action...
00885           if (!jetI->minheap_update_needed()) {
00886             jetI->label_minheap_update_needed();
00887             jets_for_minheap.push_back(jetI);}
00888           // now go over tiles that are neighbours of I (include own tile)
00889           for (Tile ** near_tile  = tile_ptr->begin_tiles; 
00890                        near_tile != tile_ptr->end_tiles; near_tile++) {
00891             // and then over the contents of that tile
00892             for (TiledJet * jetJ  = (*near_tile)->head; 
00893                             jetJ != NULL; jetJ = jetJ->next) {
00894               double dist = _bj_dist(jetI,jetJ);
00895               if (dist < jetI->NN_dist && jetJ != jetI) {
00896                 jetI->NN_dist = dist; jetI->NN = jetJ;
00897               }
00898             }
00899           }
00900         }
00901         // check whether new jetB is closer than jetI's current NN and
00902         // if jetI is closer than jetB's current (evolving) nearest
00903         // neighbour. Where relevant update things
00904         if (jetB != NULL) {
00905           double dist = _bj_dist(jetI,jetB);
00906           if (dist < jetI->NN_dist) {
00907             if (jetI != jetB) {
00908               jetI->NN_dist = dist;
00909               jetI->NN = jetB;
00910               // label jetI as needing heap action...
00911               if (!jetI->minheap_update_needed()) {
00912                 jetI->label_minheap_update_needed();
00913                 jets_for_minheap.push_back(jetI);}
00914             }
00915           }
00916           if (dist < jetB->NN_dist) {
00917             if (jetI != jetB) {
00918               jetB->NN_dist = dist;
00919               jetB->NN      = jetI;}
00920           }
00921         }
00922       }
00923     }
00924 
00925     // deal with jets whose minheap entry needs updating
00926     while (jets_for_minheap.size() > 0) {
00927       TiledJet * jetI = jets_for_minheap.back(); 
00928       jets_for_minheap.pop_back();
00929       minheap.update(jetI-head, _bj_diJ(jetI));
00930       jetI->label_minheap_update_done();
00931     }
00932     n--;
00933   }
00934 
00935   // final cleaning up;
00936   delete[] briefjets;
00937 }
00938 
00939 
00940 FASTJET_END_NAMESPACE
00941 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends