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 "fastjet/PseudoJet.hh"
00040 #include "fastjet/ClusterSequence.hh"
00041 #include "fastjet/internal/MinHeap.hh"
00042
00043 FASTJET_BEGIN_NAMESPACE
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
00054
00055 tile->head = jet->next;
00056 } else {
00057
00058 jet->previous->next = jet->next;
00059 }
00060 if (jet->next != NULL) {
00061
00062 jet->next->previous = jet->previous;
00063 }
00064 }
00065
00066
00085 void ClusterSequence::_initialise_tiles() {
00086
00087
00088 _tile_size_eta = _Rparam;
00089 _n_tiles_phi = int(floor(twopi/_Rparam));
00090 _tile_size_phi = twopi / _n_tiles_phi;
00091
00092
00093 _tiles_eta_min = 0.0;
00094 _tiles_eta_max = 0.0;
00095
00096 const double maxrap = 7.0;
00097
00098
00099 for(unsigned int i = 0; i < _jets.size(); i++) {
00100 double eta = _jets[i].rap();
00101
00102
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
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
00116 _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
00117
00118
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
00123 tile->head = NULL;
00124 tile->begin_tiles[0] = tile;
00125 Tile ** pptile = & (tile->begin_tiles[0]);
00126 pptile++;
00127
00128
00129 tile->surrounding_tiles = pptile;
00130 if (ieta > _tiles_ieta_min) {
00131
00132
00133
00134 for (int idphi = -1; idphi <=+1; idphi++) {
00135 *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
00136 pptile++;
00137 }
00138 }
00139
00140 *pptile = & _tiles[_tile_index(ieta,iphi-1)];
00141 pptile++;
00142
00143 tile->RH_tiles = pptile;
00144 *pptile = & _tiles[_tile_index(ieta,iphi+1)];
00145 pptile++;
00146
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
00154 tile->end_tiles = pptile;
00155
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
00171 ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
00172
00173 if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
00174 ieta = _tiles_ieta_max-_tiles_ieta_min;}
00175 }
00176
00177
00178
00179
00180 iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
00181 return (iphi + ieta * _n_tiles_phi);
00182 }
00183
00184
00185
00186
00187
00188 inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
00189 const int _jets_index) {
00190
00191 _bj_set_jetinfo<>(jet, _jets_index);
00192
00193
00194
00195
00196 jet->tile_index = _tile_index(jet->eta, jet->phi);
00197
00198
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
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
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
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
00275
00276 vector<int> tile_union(3*n_tile_neighbours);
00277
00278
00279 for (int i = 0; i< n; i++) {
00280 _tj_set_jetinfo(jetA, i);
00281
00282 jetA++;
00283 }
00284 TiledJet * tail = jetA;
00285 TiledJet * head = briefjets;
00286
00287
00288 vector<Tile>::const_iterator tile;
00289 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00290
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
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
00311
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++;
00317 }
00318
00319
00320 int history_location = n-1;
00321 while (tail != head) {
00322
00323
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
00331 history_location++;
00332 jetA = & briefjets[diJ_min_jet];
00333 jetB = jetA->NN;
00334
00335 diJ_min *= _invR2;
00336
00337
00338
00339
00340
00341 if (jetB != NULL) {
00342
00343
00344
00345
00346
00347 if (jetA < jetB) {swap(jetA,jetB);}
00348
00349 int nn;
00350 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 _bj_remove_from_tiles(jetA);
00367 oldB = * jetB;
00368 _bj_remove_from_tiles(jetB);
00369 _tj_set_jetinfo(jetB, nn);
00370 } else {
00371
00372 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00373
00374
00375
00376
00377
00378 _bj_remove_from_tiles(jetA);
00379 }
00380
00381
00382
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
00399 sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
00400
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
00413
00414 tail--; n--;
00415 if (jetA == tail) {
00416
00417 } else {
00418
00419 *jetA = *tail;
00420 diJ[jetA - head] = diJ[tail-head];
00421
00422
00423
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
00433
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
00438 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00439 jetI->NN_dist = _R2;
00440 jetI->NN = NULL;
00441
00442 for (Tile ** near_tile = tile->begin_tiles;
00443 near_tile != tile->end_tiles; near_tile++) {
00444
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);
00454 }
00455
00456
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);
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
00478
00479
00480 for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
00481 near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
00482
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
00490
00491
00492
00493
00494 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00495
00496
00497 }
00498
00499
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
00518
00519 vector<int> tile_union(3*n_tile_neighbours);
00520
00521
00522 for (int i = 0; i< n; i++) {
00523 _tj_set_jetinfo(jetA, i);
00524
00525 jetA++;
00526 }
00527 TiledJet * head = briefjets;
00528
00529
00530 vector<Tile>::const_iterator tile;
00531 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00532
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
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
00551
00552 }
00553
00554
00555
00556
00557
00558 struct diJ_plus_link {
00559 double diJ;
00560 TiledJet * jet;
00561
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);
00567 diJ[i].jet = jetA;
00568 jetA->diJ_posn = i;
00569
00570 jetA++;
00571 }
00572
00573
00574 int history_location = n-1;
00575 while (n > 0) {
00576
00577
00578 diJ_plus_link * best, *stop;
00579
00580
00581
00582 double diJ_min = diJ[0].diJ;
00583 best = diJ;
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
00590 history_location++;
00591 jetA = best->jet;
00592 jetB = jetA->NN;
00593
00594 diJ_min *= _invR2;
00595
00596 if (jetB != NULL) {
00597
00598
00599
00600
00601
00602 if (jetA < jetB) {swap(jetA,jetB);}
00603
00604 int nn;
00605 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 _bj_remove_from_tiles(jetA);
00621 oldB = * jetB;
00622 _bj_remove_from_tiles(jetB);
00623 _tj_set_jetinfo(jetB, nn);
00624
00625 } else {
00626
00627
00628 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00629
00630
00631
00632 _bj_remove_from_tiles(jetA);
00633 }
00634
00635
00636
00637
00638
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
00655
00656 n--;
00657
00658
00659 diJ[n].jet->diJ_posn = jetA->diJ_posn;
00660 diJ[jetA->diJ_posn] = diJ[n];
00661
00662
00663
00664
00665 for (int itile = 0; itile < n_near_tiles; itile++) {
00666 Tile * tile = &_tiles[tile_union[itile]];
00667 tile->tagged = false;
00668
00669 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00670
00671 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00672 jetI->NN_dist = _R2;
00673 jetI->NN = NULL;
00674
00675 for (Tile ** near_tile = tile->begin_tiles;
00676 near_tile != tile->end_tiles; near_tile++) {
00677
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);
00687 }
00688
00689
00690
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);
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
00710 if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
00711
00712 }
00713
00714
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
00735
00736 vector<int> tile_union(3*n_tile_neighbours);
00737
00738
00739 for (int i = 0; i< n; i++) {
00740 _tj_set_jetinfo(jetA, i);
00741
00742 jetA++;
00743 }
00744 TiledJet * head = briefjets;
00745
00746
00747 vector<Tile>::const_iterator tile;
00748 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00749
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
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
00768
00769 }
00770
00771
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
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
00797 vector<TiledJet *> jets_for_minheap;
00798 jets_for_minheap.reserve(n);
00799
00800
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
00808 history_location++;
00809 jetB = jetA->NN;
00810
00811 if (jetB != NULL) {
00812
00813
00814
00815
00816
00817 if (jetA < jetB) {swap(jetA,jetB);}
00818
00819 int nn;
00820 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00821
00822
00823 _bj_remove_from_tiles(jetA);
00824 oldB = * jetB;
00825 _bj_remove_from_tiles(jetB);
00826 _tj_set_jetinfo(jetB, nn);
00827
00828 } else {
00829
00830
00831 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00832 _bj_remove_from_tiles(jetA);
00833 }
00834
00835
00836 minheap.remove(jetA-head);
00837
00838
00839
00840
00841
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
00856 jetB->label_minheap_update_needed();
00857 jets_for_minheap.push_back(jetB);
00858 }
00859
00860
00861
00862
00863
00864 for (int itile = 0; itile < n_near_tiles; itile++) {
00865 Tile * tile = &_tiles[tile_union[itile]];
00866 tile->tagged = false;
00867
00868 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00869
00870 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00871 jetI->NN_dist = _R2;
00872 jetI->NN = NULL;
00873
00874 if (!jetI->minheap_update_needed()) {
00875 jetI->label_minheap_update_needed();
00876 jets_for_minheap.push_back(jetI);}
00877
00878 for (Tile ** near_tile = tile->begin_tiles;
00879 near_tile != tile->end_tiles; near_tile++) {
00880
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
00891
00892
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
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
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
00925 delete[] briefjets;
00926 }
00927
00928
00929 FASTJET_END_NAMESPACE
00930