00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include<iostream>
00037 #include<vector>
00038 #include<cmath>
00039 #include<algorithm>
00040 #include "fastjet/PseudoJet.hh"
00041 #include "fastjet/ClusterSequence.hh"
00042 #include "fastjet/internal/MinHeap.hh"
00043
00044 FASTJET_BEGIN_NAMESPACE
00045
00046 using namespace std;
00047
00048
00049
00050 void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) {
00051 Tile * tile = & _tiles[jet->tile_index];
00052
00053 if (jet->previous == NULL) {
00054
00055
00056 tile->head = jet->next;
00057 } else {
00058
00059 jet->previous->next = jet->next;
00060 }
00061 if (jet->next != NULL) {
00062
00063 jet->next->previous = jet->previous;
00064 }
00065 }
00066
00067
00086 void ClusterSequence::_initialise_tiles() {
00087
00088
00089
00090 double default_size = max(0.1,_Rparam);
00091 _tile_size_eta = default_size;
00092 _n_tiles_phi = int(floor(twopi/default_size));
00093 _tile_size_phi = twopi / _n_tiles_phi;
00094
00095
00096 _tiles_eta_min = 0.0;
00097 _tiles_eta_max = 0.0;
00098
00099 const double maxrap = 7.0;
00100
00101
00102 for(unsigned int i = 0; i < _jets.size(); i++) {
00103 double eta = _jets[i].rap();
00104
00105
00106 if (abs(eta) < maxrap) {
00107 if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
00108 if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
00109 }
00110 }
00111
00112
00113 _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
00114 _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
00115 _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
00116 _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
00117
00118
00119 _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
00120
00121
00122 for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
00123 for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
00124 Tile * tile = & _tiles[_tile_index(ieta,iphi)];
00125
00126 tile->head = NULL;
00127 tile->begin_tiles[0] = tile;
00128 Tile ** pptile = & (tile->begin_tiles[0]);
00129 pptile++;
00130
00131
00132 tile->surrounding_tiles = pptile;
00133 if (ieta > _tiles_ieta_min) {
00134
00135
00136
00137 for (int idphi = -1; idphi <=+1; idphi++) {
00138 *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
00139 pptile++;
00140 }
00141 }
00142
00143 *pptile = & _tiles[_tile_index(ieta,iphi-1)];
00144 pptile++;
00145
00146 tile->RH_tiles = pptile;
00147 *pptile = & _tiles[_tile_index(ieta,iphi+1)];
00148 pptile++;
00149
00150 if (ieta < _tiles_ieta_max) {
00151 for (int idphi = -1; idphi <= +1; idphi++) {
00152 *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
00153 pptile++;
00154 }
00155 }
00156
00157 tile->end_tiles = pptile;
00158
00159 tile->tagged = false;
00160 }
00161 }
00162
00163 }
00164
00165
00166
00168 int ClusterSequence::_tile_index(const double & eta, const double & phi) const {
00169 int ieta, iphi;
00170 if (eta <= _tiles_eta_min) {ieta = 0;}
00171 else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
00172 else {
00173
00174 ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
00175
00176 if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
00177 ieta = _tiles_ieta_max-_tiles_ieta_min;}
00178 }
00179
00180
00181
00182
00183 iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
00184 return (iphi + ieta * _n_tiles_phi);
00185 }
00186
00187
00188
00189
00190
00191 inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
00192 const int _jets_index) {
00193
00194 _bj_set_jetinfo<>(jet, _jets_index);
00195
00196
00197
00198
00199 jet->tile_index = _tile_index(jet->eta, jet->phi);
00200
00201
00202 Tile * tile = &_tiles[jet->tile_index];
00203 jet->previous = NULL;
00204 jet->next = tile->head;
00205 if (jet->next != NULL) {jet->next->previous = jet;}
00206 tile->head = jet;
00207 }
00208
00209
00210
00212 void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
00213 for (vector<Tile>::const_iterator tile = _tiles.begin();
00214 tile < _tiles.end(); tile++) {
00215 cout << "Tile " << tile - _tiles.begin()<<" = ";
00216 vector<int> list;
00217 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00218 list.push_back(jetI-briefjets);
00219
00220 }
00221 sort(list.begin(),list.end());
00222 for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
00223 cout <<"\n";
00224 }
00225 }
00226
00227
00228
00235 void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index,
00236 vector<int> & tile_union, int & n_near_tiles) const {
00237 for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
00238 near_tile != _tiles[tile_index].end_tiles; near_tile++){
00239
00240 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
00241 n_near_tiles++;
00242 }
00243 }
00244
00245
00246
00250 inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
00251 const int tile_index,
00252 vector<int> & tile_union, int & n_near_tiles) {
00253 for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
00254 near_tile != _tiles[tile_index].end_tiles; near_tile++){
00255 if (! (*near_tile)->tagged) {
00256 (*near_tile)->tagged = true;
00257
00258 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
00259 n_near_tiles++;
00260 }
00261 }
00262 }
00263
00264
00265
00267 void ClusterSequence::_tiled_N2_cluster() {
00268
00269 _initialise_tiles();
00270
00271 int n = _jets.size();
00272 TiledJet * briefjets = new TiledJet[n];
00273 TiledJet * jetA = briefjets, * jetB;
00274 TiledJet oldB;
00275
00276
00277
00278
00279 vector<int> tile_union(3*n_tile_neighbours);
00280
00281
00282 for (int i = 0; i< n; i++) {
00283 _tj_set_jetinfo(jetA, i);
00284
00285 jetA++;
00286 }
00287 TiledJet * tail = jetA;
00288 TiledJet * head = briefjets;
00289
00290
00291 vector<Tile>::const_iterator tile;
00292 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00293
00294 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00295 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00296 double dist = _bj_dist(jetA,jetB);
00297 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00298 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00299 }
00300 }
00301
00302 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00303 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00304 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00305 double dist = _bj_dist(jetA,jetB);
00306 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00307 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00308 }
00309 }
00310 }
00311 }
00312
00313
00314
00315 double * diJ = new double[n];
00316 jetA = head;
00317 for (int i = 0; i < n; i++) {
00318 diJ[i] = _bj_diJ(jetA);
00319 jetA++;
00320 }
00321
00322
00323 int history_location = n-1;
00324 while (tail != head) {
00325
00326
00327 double diJ_min = diJ[0];
00328 int diJ_min_jet = 0;
00329 for (int i = 1; i < n; i++) {
00330 if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
00331 }
00332
00333
00334 history_location++;
00335 jetA = & briefjets[diJ_min_jet];
00336 jetB = jetA->NN;
00337
00338 diJ_min *= _invR2;
00339
00340
00341
00342
00343
00344 if (jetB != NULL) {
00345
00346
00347
00348
00349
00350 if (jetA < jetB) {swap(jetA,jetB);}
00351
00352 int nn;
00353 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369 _bj_remove_from_tiles(jetA);
00370 oldB = * jetB;
00371 _bj_remove_from_tiles(jetB);
00372 _tj_set_jetinfo(jetB, nn);
00373 } else {
00374
00375 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00376
00377
00378
00379
00380
00381 _bj_remove_from_tiles(jetA);
00382 }
00383
00384
00385
00386 int n_near_tiles = 0;
00387 _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
00388 if (jetB != NULL) {
00389 bool sort_it = false;
00390 if (jetB->tile_index != jetA->tile_index) {
00391 sort_it = true;
00392 _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
00393 }
00394 if (oldB.tile_index != jetA->tile_index &&
00395 oldB.tile_index != jetB->tile_index) {
00396 sort_it = true;
00397 _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
00398 }
00399
00400 if (sort_it) {
00401
00402 sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
00403
00404 int nnn = 1;
00405 for (int i = 1; i < n_near_tiles; i++) {
00406 if (tile_union[i] != tile_union[nnn-1]) {
00407 tile_union[nnn] = tile_union[i];
00408 nnn++;
00409 }
00410 }
00411 n_near_tiles = nnn;
00412 }
00413 }
00414
00415
00416
00417 tail--; n--;
00418 if (jetA == tail) {
00419
00420 } else {
00421
00422 *jetA = *tail;
00423 diJ[jetA - head] = diJ[tail-head];
00424
00425
00426
00427 if (jetA->previous == NULL) {
00428 _tiles[jetA->tile_index].head = jetA;
00429 } else {
00430 jetA->previous->next = jetA;
00431 }
00432 if (jetA->next != NULL) {jetA->next->previous = jetA;}
00433 }
00434
00435
00436
00437 for (int itile = 0; itile < n_near_tiles; itile++) {
00438 Tile * tile = &_tiles[tile_union[itile]];
00439 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00440
00441 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00442 jetI->NN_dist = _R2;
00443 jetI->NN = NULL;
00444
00445 for (Tile ** near_tile = tile->begin_tiles;
00446 near_tile != tile->end_tiles; near_tile++) {
00447
00448 for (TiledJet * jetJ = (*near_tile)->head;
00449 jetJ != NULL; jetJ = jetJ->next) {
00450 double dist = _bj_dist(jetI,jetJ);
00451 if (dist < jetI->NN_dist && jetJ != jetI) {
00452 jetI->NN_dist = dist; jetI->NN = jetJ;
00453 }
00454 }
00455 }
00456 diJ[jetI-head] = _bj_diJ(jetI);
00457 }
00458
00459
00460 if (jetB != NULL) {
00461 double dist = _bj_dist(jetI,jetB);
00462 if (dist < jetI->NN_dist) {
00463 if (jetI != jetB) {
00464 jetI->NN_dist = dist;
00465 jetI->NN = jetB;
00466 diJ[jetI-head] = _bj_diJ(jetI);
00467 }
00468 }
00469 if (dist < jetB->NN_dist) {
00470 if (jetI != jetB) {
00471 jetB->NN_dist = dist;
00472 jetB->NN = jetI;}
00473 }
00474 }
00475 }
00476 }
00477
00478
00479 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00480
00481
00482
00483 for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
00484 near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
00485
00486 for (TiledJet * jetJ = (*near_tile)->head;
00487 jetJ != NULL; jetJ = jetJ->next) {
00488 if (jetJ->NN == tail) {jetJ->NN = jetA;}
00489 }
00490 }
00491
00492
00493
00494
00495
00496
00497 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00498
00499
00500 }
00501
00502
00503 delete[] diJ;
00504 delete[] briefjets;
00505 }
00506
00507
00508
00510 void ClusterSequence::_faster_tiled_N2_cluster() {
00511
00512 _initialise_tiles();
00513
00514 int n = _jets.size();
00515 TiledJet * briefjets = new TiledJet[n];
00516 TiledJet * jetA = briefjets, * jetB;
00517 TiledJet oldB;
00518
00519
00520
00521
00522 vector<int> tile_union(3*n_tile_neighbours);
00523
00524
00525 for (int i = 0; i< n; i++) {
00526 _tj_set_jetinfo(jetA, i);
00527
00528 jetA++;
00529 }
00530 TiledJet * head = briefjets;
00531
00532
00533 vector<Tile>::const_iterator tile;
00534 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00535
00536 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00537 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00538 double dist = _bj_dist(jetA,jetB);
00539 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00540 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00541 }
00542 }
00543
00544 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00545 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00546 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00547 double dist = _bj_dist(jetA,jetB);
00548 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00549 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00550 }
00551 }
00552 }
00553
00554
00555 }
00556
00557
00558
00559
00560
00561 struct diJ_plus_link {
00562 double diJ;
00563 TiledJet * jet;
00564
00565 };
00566 diJ_plus_link * diJ = new diJ_plus_link[n];
00567 jetA = head;
00568 for (int i = 0; i < n; i++) {
00569 diJ[i].diJ = _bj_diJ(jetA);
00570 diJ[i].jet = jetA;
00571 jetA->diJ_posn = i;
00572
00573 jetA++;
00574 }
00575
00576
00577 int history_location = n-1;
00578 while (n > 0) {
00579
00580
00581 diJ_plus_link * best, *stop;
00582
00583
00584
00585 double diJ_min = diJ[0].diJ;
00586 best = diJ;
00587 stop = diJ+n;
00588 for (diJ_plus_link * here = diJ+1; here != stop; here++) {
00589 if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
00590 }
00591
00592
00593 history_location++;
00594 jetA = best->jet;
00595 jetB = jetA->NN;
00596
00597 diJ_min *= _invR2;
00598
00599 if (jetB != NULL) {
00600
00601
00602
00603
00604
00605 if (jetA < jetB) {swap(jetA,jetB);}
00606
00607 int nn;
00608 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623 _bj_remove_from_tiles(jetA);
00624 oldB = * jetB;
00625 _bj_remove_from_tiles(jetB);
00626 _tj_set_jetinfo(jetB, nn);
00627
00628 } else {
00629
00630
00631 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00632
00633
00634
00635 _bj_remove_from_tiles(jetA);
00636 }
00637
00638
00639
00640
00641
00642 int n_near_tiles = 0;
00643 _add_untagged_neighbours_to_tile_union(jetA->tile_index,
00644 tile_union, n_near_tiles);
00645 if (jetB != NULL) {
00646 if (jetB->tile_index != jetA->tile_index) {
00647 _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00648 tile_union,n_near_tiles);
00649 }
00650 if (oldB.tile_index != jetA->tile_index &&
00651 oldB.tile_index != jetB->tile_index) {
00652 _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00653 tile_union,n_near_tiles);
00654 }
00655 }
00656
00657
00658
00659 n--;
00660
00661
00662 diJ[n].jet->diJ_posn = jetA->diJ_posn;
00663 diJ[jetA->diJ_posn] = diJ[n];
00664
00665
00666
00667
00668 for (int itile = 0; itile < n_near_tiles; itile++) {
00669 Tile * tile = &_tiles[tile_union[itile]];
00670 tile->tagged = false;
00671
00672 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00673
00674 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00675 jetI->NN_dist = _R2;
00676 jetI->NN = NULL;
00677
00678 for (Tile ** near_tile = tile->begin_tiles;
00679 near_tile != tile->end_tiles; near_tile++) {
00680
00681 for (TiledJet * jetJ = (*near_tile)->head;
00682 jetJ != NULL; jetJ = jetJ->next) {
00683 double dist = _bj_dist(jetI,jetJ);
00684 if (dist < jetI->NN_dist && jetJ != jetI) {
00685 jetI->NN_dist = dist; jetI->NN = jetJ;
00686 }
00687 }
00688 }
00689 diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI);
00690 }
00691
00692
00693
00694 if (jetB != NULL) {
00695 double dist = _bj_dist(jetI,jetB);
00696 if (dist < jetI->NN_dist) {
00697 if (jetI != jetB) {
00698 jetI->NN_dist = dist;
00699 jetI->NN = jetB;
00700 diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI);
00701 }
00702 }
00703 if (dist < jetB->NN_dist) {
00704 if (jetI != jetB) {
00705 jetB->NN_dist = dist;
00706 jetB->NN = jetI;}
00707 }
00708 }
00709 }
00710 }
00711
00712
00713 if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
00714
00715 }
00716
00717
00718 delete[] diJ;
00719 delete[] briefjets;
00720 }
00721
00722
00723
00724
00727 void ClusterSequence::_minheap_faster_tiled_N2_cluster() {
00728
00729 _initialise_tiles();
00730
00731 int n = _jets.size();
00732 TiledJet * briefjets = new TiledJet[n];
00733 TiledJet * jetA = briefjets, * jetB;
00734 TiledJet oldB;
00735
00736
00737
00738
00739 vector<int> tile_union(3*n_tile_neighbours);
00740
00741
00742 for (int i = 0; i< n; i++) {
00743 _tj_set_jetinfo(jetA, i);
00744
00745 jetA++;
00746 }
00747 TiledJet * head = briefjets;
00748
00749
00750 vector<Tile>::const_iterator tile;
00751 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00752
00753 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00754 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00755 double dist = _bj_dist(jetA,jetB);
00756 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00757 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00758 }
00759 }
00760
00761 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00762 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00763 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00764 double dist = _bj_dist(jetA,jetB);
00765 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00766 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00767 }
00768 }
00769 }
00770
00771
00772 }
00773
00774
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793 vector<double> diJs(n);
00794 for (int i = 0; i < n; i++) {
00795 diJs[i] = _bj_diJ(&briefjets[i]);
00796 briefjets[i].label_minheap_update_done();
00797 }
00798 MinHeap minheap(diJs);
00799
00800 vector<TiledJet *> jets_for_minheap;
00801 jets_for_minheap.reserve(n);
00802
00803
00804 int history_location = n-1;
00805 while (n > 0) {
00806
00807 double diJ_min = minheap.minval() *_invR2;
00808 jetA = head + minheap.minloc();
00809
00810
00811 history_location++;
00812 jetB = jetA->NN;
00813
00814 if (jetB != NULL) {
00815
00816
00817
00818
00819
00820 if (jetA < jetB) {swap(jetA,jetB);}
00821
00822 int nn;
00823 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00824
00825
00826 _bj_remove_from_tiles(jetA);
00827 oldB = * jetB;
00828 _bj_remove_from_tiles(jetB);
00829 _tj_set_jetinfo(jetB, nn);
00830
00831 } else {
00832
00833
00834 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00835 _bj_remove_from_tiles(jetA);
00836 }
00837
00838
00839 minheap.remove(jetA-head);
00840
00841
00842
00843
00844
00845 int n_near_tiles = 0;
00846 _add_untagged_neighbours_to_tile_union(jetA->tile_index,
00847 tile_union, n_near_tiles);
00848 if (jetB != NULL) {
00849 if (jetB->tile_index != jetA->tile_index) {
00850 _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00851 tile_union,n_near_tiles);
00852 }
00853 if (oldB.tile_index != jetA->tile_index &&
00854 oldB.tile_index != jetB->tile_index) {
00855 _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00856 tile_union,n_near_tiles);
00857 }
00858
00859 jetB->label_minheap_update_needed();
00860 jets_for_minheap.push_back(jetB);
00861 }
00862
00863
00864
00865
00866
00867 for (int itile = 0; itile < n_near_tiles; itile++) {
00868 Tile * tile = &_tiles[tile_union[itile]];
00869 tile->tagged = false;
00870
00871 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00872
00873 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00874 jetI->NN_dist = _R2;
00875 jetI->NN = NULL;
00876
00877 if (!jetI->minheap_update_needed()) {
00878 jetI->label_minheap_update_needed();
00879 jets_for_minheap.push_back(jetI);}
00880
00881 for (Tile ** near_tile = tile->begin_tiles;
00882 near_tile != tile->end_tiles; near_tile++) {
00883
00884 for (TiledJet * jetJ = (*near_tile)->head;
00885 jetJ != NULL; jetJ = jetJ->next) {
00886 double dist = _bj_dist(jetI,jetJ);
00887 if (dist < jetI->NN_dist && jetJ != jetI) {
00888 jetI->NN_dist = dist; jetI->NN = jetJ;
00889 }
00890 }
00891 }
00892 }
00893
00894
00895
00896 if (jetB != NULL) {
00897 double dist = _bj_dist(jetI,jetB);
00898 if (dist < jetI->NN_dist) {
00899 if (jetI != jetB) {
00900 jetI->NN_dist = dist;
00901 jetI->NN = jetB;
00902
00903 if (!jetI->minheap_update_needed()) {
00904 jetI->label_minheap_update_needed();
00905 jets_for_minheap.push_back(jetI);}
00906 }
00907 }
00908 if (dist < jetB->NN_dist) {
00909 if (jetI != jetB) {
00910 jetB->NN_dist = dist;
00911 jetB->NN = jetI;}
00912 }
00913 }
00914 }
00915 }
00916
00917
00918 while (jets_for_minheap.size() > 0) {
00919 TiledJet * jetI = jets_for_minheap.back();
00920 jets_for_minheap.pop_back();
00921 minheap.update(jetI-head, _bj_diJ(jetI));
00922 jetI->label_minheap_update_done();
00923 }
00924 n--;
00925 }
00926
00927
00928 delete[] briefjets;
00929 }
00930
00931
00932 FASTJET_END_NAMESPACE
00933