FastJet 3.0.2
|
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