ClusterSequence_TiledN2.cc

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

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2