38 #include "fastjet/PseudoJet.hh"
39 #include "fastjet/ClusterSequence.hh"
40 #include "fastjet/internal/MinHeap.hh"
42 FASTJET_BEGIN_NAMESPACE
48 void ClusterSequence::_bj_remove_from_tiles(TiledJet *
const jet) {
49 Tile * tile = & _tiles[jet->tile_index];
51 if (jet->previous == NULL) {
54 tile->head = jet->next;
57 jet->previous->next = jet->next;
59 if (jet->next != NULL) {
61 jet->next->previous = jet->previous;
84 void ClusterSequence::_initialise_tiles() {
88 double default_size = max(0.1,_Rparam);
89 _tile_size_eta = default_size;
93 _n_tiles_phi = max(3,
int(floor(twopi/default_size)));
94 _tile_size_phi = twopi / _n_tiles_phi;
100 const double maxrap = 7.0;
103 for(
unsigned int i = 0; i < _jets.size(); i++) {
104 double eta = _jets[i].rap();
107 if (abs(eta) < maxrap) {
108 if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
109 if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
114 _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
115 _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
116 _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
117 _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
120 _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
123 for (
int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
124 for (
int iphi = 0; iphi < _n_tiles_phi; iphi++) {
125 Tile * tile = & _tiles[_tile_index(ieta,iphi)];
128 tile->begin_tiles[0] = tile;
129 Tile ** pptile = & (tile->begin_tiles[0]);
133 tile->surrounding_tiles = pptile;
134 if (ieta > _tiles_ieta_min) {
138 for (
int idphi = -1; idphi <=+1; idphi++) {
139 *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
144 *pptile = & _tiles[_tile_index(ieta,iphi-1)];
147 tile->RH_tiles = pptile;
148 *pptile = & _tiles[_tile_index(ieta,iphi+1)];
151 if (ieta < _tiles_ieta_max) {
152 for (
int idphi = -1; idphi <= +1; idphi++) {
153 *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
158 tile->end_tiles = pptile;
160 tile->tagged =
false;
169 int ClusterSequence::_tile_index(
const double & eta,
const double & phi)
const {
171 if (eta <= _tiles_eta_min) {ieta = 0;}
172 else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
175 ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
177 if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
178 ieta = _tiles_ieta_max-_tiles_ieta_min;}
184 iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
185 return (iphi + ieta * _n_tiles_phi);
192 inline void ClusterSequence::_tj_set_jetinfo( TiledJet *
const jet,
193 const int _jets_index) {
195 _bj_set_jetinfo<>(jet, _jets_index);
200 jet->tile_index = _tile_index(jet->eta, jet->phi);
203 Tile * tile = &_tiles[jet->tile_index];
204 jet->previous = NULL;
205 jet->next = tile->head;
206 if (jet->next != NULL) {jet->next->previous = jet;}
213 void ClusterSequence::_print_tiles(TiledJet * briefjets )
const {
214 for (vector<Tile>::const_iterator tile = _tiles.begin();
215 tile < _tiles.end(); tile++) {
216 cout <<
"Tile " << tile - _tiles.begin()<<
" = ";
218 for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
219 list.push_back(jetI-briefjets);
222 sort(list.begin(),list.end());
223 for (
unsigned int i = 0; i < list.size(); i++) {cout <<
" "<<list[i];}
236 void ClusterSequence::_add_neighbours_to_tile_union(
const int tile_index,
237 vector<int> & tile_union,
int & n_near_tiles)
const {
238 for (Tile *
const * near_tile = _tiles[tile_index].begin_tiles;
239 near_tile != _tiles[tile_index].end_tiles; near_tile++){
241 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
251 inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
252 const int tile_index,
253 vector<int> & tile_union,
int & n_near_tiles) {
254 for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
255 near_tile != _tiles[tile_index].end_tiles; near_tile++){
256 if (! (*near_tile)->tagged) {
257 (*near_tile)->tagged =
true;
259 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
268 void ClusterSequence::_tiled_N2_cluster() {
272 int n = _jets.size();
273 TiledJet * briefjets =
new TiledJet[n];
274 TiledJet * jetA = briefjets, * jetB;
280 vector<int> tile_union(3*n_tile_neighbours);
283 for (
int i = 0; i< n; i++) {
284 _tj_set_jetinfo(jetA, i);
288 TiledJet * tail = jetA;
289 TiledJet * head = briefjets;
292 vector<Tile>::const_iterator tile;
293 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
295 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
296 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
297 double dist = _bj_dist(jetA,jetB);
298 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
299 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
303 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
304 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
305 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
306 double dist = _bj_dist(jetA,jetB);
307 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
308 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
316 double * diJ =
new double[n];
318 for (
int i = 0; i < n; i++) {
319 diJ[i] = _bj_diJ(jetA);
324 int history_location = n-1;
325 while (tail != head) {
328 double diJ_min = diJ[0];
330 for (
int i = 1; i < n; i++) {
331 if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
336 jetA = & briefjets[diJ_min_jet];
351 if (jetA < jetB) {std::swap(jetA,jetB);}
354 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
370 _bj_remove_from_tiles(jetA);
372 _bj_remove_from_tiles(jetB);
373 _tj_set_jetinfo(jetB, nn);
376 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
382 _bj_remove_from_tiles(jetA);
387 int n_near_tiles = 0;
388 _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
390 bool sort_it =
false;
391 if (jetB->tile_index != jetA->tile_index) {
393 _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
395 if (oldB.tile_index != jetA->tile_index &&
396 oldB.tile_index != jetB->tile_index) {
398 _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
403 sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
406 for (
int i = 1; i < n_near_tiles; i++) {
407 if (tile_union[i] != tile_union[nnn-1]) {
408 tile_union[nnn] = tile_union[i];
424 diJ[jetA - head] = diJ[tail-head];
428 if (jetA->previous == NULL) {
429 _tiles[jetA->tile_index].head = jetA;
431 jetA->previous->next = jetA;
433 if (jetA->next != NULL) {jetA->next->previous = jetA;}
438 for (
int itile = 0; itile < n_near_tiles; itile++) {
439 Tile * tile_ptr = &_tiles[tile_union[itile]];
440 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
442 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
446 for (Tile ** near_tile = tile_ptr->begin_tiles;
447 near_tile != tile_ptr->end_tiles; near_tile++) {
449 for (TiledJet * jetJ = (*near_tile)->head;
450 jetJ != NULL; jetJ = jetJ->next) {
451 double dist = _bj_dist(jetI,jetJ);
452 if (dist < jetI->NN_dist && jetJ != jetI) {
453 jetI->NN_dist = dist; jetI->NN = jetJ;
457 diJ[jetI-head] = _bj_diJ(jetI);
462 double dist = _bj_dist(jetI,jetB);
463 if (dist < jetI->NN_dist) {
465 jetI->NN_dist = dist;
467 diJ[jetI-head] = _bj_diJ(jetI);
470 if (dist < jetB->NN_dist) {
472 jetB->NN_dist = dist;
480 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
484 for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
485 near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
487 for (TiledJet * jetJ = (*near_tile)->head;
488 jetJ != NULL; jetJ = jetJ->next) {
489 if (jetJ->NN == tail) {jetJ->NN = jetA;}
498 if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
511 void ClusterSequence::_faster_tiled_N2_cluster() {
515 int n = _jets.size();
516 TiledJet * briefjets =
new TiledJet[n];
517 TiledJet * jetA = briefjets, * jetB;
523 vector<int> tile_union(3*n_tile_neighbours);
526 for (
int i = 0; i< n; i++) {
527 _tj_set_jetinfo(jetA, i);
531 TiledJet * head = briefjets;
534 vector<Tile>::const_iterator tile;
535 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
537 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
538 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
539 double dist = _bj_dist(jetA,jetB);
540 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
541 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
545 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
546 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
547 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
548 double dist = _bj_dist(jetA,jetB);
549 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
550 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
562 struct diJ_plus_link {
567 diJ_plus_link * diJ =
new diJ_plus_link[n];
569 for (
int i = 0; i < n; i++) {
570 diJ[i].diJ = _bj_diJ(jetA);
578 int history_location = n-1;
582 diJ_plus_link * best, *stop;
586 double diJ_min = diJ[0].diJ;
589 for (diJ_plus_link * here = diJ+1; here != stop; here++) {
590 if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
606 if (jetA < jetB) {std::swap(jetA,jetB);}
609 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
624 _bj_remove_from_tiles(jetA);
626 _bj_remove_from_tiles(jetB);
627 _tj_set_jetinfo(jetB, nn);
632 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
636 _bj_remove_from_tiles(jetA);
643 int n_near_tiles = 0;
644 _add_untagged_neighbours_to_tile_union(jetA->tile_index,
645 tile_union, n_near_tiles);
647 if (jetB->tile_index != jetA->tile_index) {
648 _add_untagged_neighbours_to_tile_union(jetB->tile_index,
649 tile_union,n_near_tiles);
651 if (oldB.tile_index != jetA->tile_index &&
652 oldB.tile_index != jetB->tile_index) {
653 _add_untagged_neighbours_to_tile_union(oldB.tile_index,
654 tile_union,n_near_tiles);
663 diJ[n].jet->diJ_posn = jetA->diJ_posn;
664 diJ[jetA->diJ_posn] = diJ[n];
669 for (
int itile = 0; itile < n_near_tiles; itile++) {
670 Tile * tile_ptr = &_tiles[tile_union[itile]];
671 tile_ptr->tagged =
false;
673 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
675 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
679 for (Tile ** near_tile = tile_ptr->begin_tiles;
680 near_tile != tile_ptr->end_tiles; near_tile++) {
682 for (TiledJet * jetJ = (*near_tile)->head;
683 jetJ != NULL; jetJ = jetJ->next) {
684 double dist = _bj_dist(jetI,jetJ);
685 if (dist < jetI->NN_dist && jetJ != jetI) {
686 jetI->NN_dist = dist; jetI->NN = jetJ;
690 diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI);
696 double dist = _bj_dist(jetI,jetB);
697 if (dist < jetI->NN_dist) {
699 jetI->NN_dist = dist;
701 diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI);
704 if (dist < jetB->NN_dist) {
706 jetB->NN_dist = dist;
714 if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
728 void ClusterSequence::_minheap_faster_tiled_N2_cluster() {
732 int n = _jets.size();
733 TiledJet * briefjets =
new TiledJet[n];
734 TiledJet * jetA = briefjets, * jetB;
740 vector<int> tile_union(3*n_tile_neighbours);
743 for (
int i = 0; i< n; i++) {
744 _tj_set_jetinfo(jetA, i);
748 TiledJet * head = briefjets;
751 vector<Tile>::const_iterator tile;
752 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
754 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
755 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
756 double dist = _bj_dist(jetA,jetB);
757 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
758 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
762 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
763 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
764 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
765 double dist = _bj_dist(jetA,jetB);
766 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
767 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
794 vector<double> diJs(n);
795 for (
int i = 0; i < n; i++) {
796 diJs[i] = _bj_diJ(&briefjets[i]);
797 briefjets[i].label_minheap_update_done();
799 MinHeap minheap(diJs);
801 vector<TiledJet *> jets_for_minheap;
802 jets_for_minheap.reserve(n);
805 int history_location = n-1;
808 double diJ_min = minheap.minval() *_invR2;
809 jetA = head + minheap.minloc();
821 if (jetA < jetB) {std::swap(jetA,jetB);}
824 _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
827 _bj_remove_from_tiles(jetA);
829 _bj_remove_from_tiles(jetB);
830 _tj_set_jetinfo(jetB, nn);
835 _do_iB_recombination_step(jetA->_jets_index, diJ_min);
836 _bj_remove_from_tiles(jetA);
840 minheap.remove(jetA-head);
846 int n_near_tiles = 0;
847 _add_untagged_neighbours_to_tile_union(jetA->tile_index,
848 tile_union, n_near_tiles);
850 if (jetB->tile_index != jetA->tile_index) {
851 _add_untagged_neighbours_to_tile_union(jetB->tile_index,
852 tile_union,n_near_tiles);
854 if (oldB.tile_index != jetA->tile_index &&
855 oldB.tile_index != jetB->tile_index) {
863 _add_untagged_neighbours_to_tile_union(oldB.tile_index,
864 tile_union,n_near_tiles);
867 jetB->label_minheap_update_needed();
868 jets_for_minheap.push_back(jetB);
875 for (
int itile = 0; itile < n_near_tiles; itile++) {
876 Tile * tile_ptr = &_tiles[tile_union[itile]];
877 tile_ptr->tagged =
false;
879 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
881 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
885 if (!jetI->minheap_update_needed()) {
886 jetI->label_minheap_update_needed();
887 jets_for_minheap.push_back(jetI);}
889 for (Tile ** near_tile = tile_ptr->begin_tiles;
890 near_tile != tile_ptr->end_tiles; near_tile++) {
892 for (TiledJet * jetJ = (*near_tile)->head;
893 jetJ != NULL; jetJ = jetJ->next) {
894 double dist = _bj_dist(jetI,jetJ);
895 if (dist < jetI->NN_dist && jetJ != jetI) {
896 jetI->NN_dist = dist; jetI->NN = jetJ;
905 double dist = _bj_dist(jetI,jetB);
906 if (dist < jetI->NN_dist) {
908 jetI->NN_dist = dist;
911 if (!jetI->minheap_update_needed()) {
912 jetI->label_minheap_update_needed();
913 jets_for_minheap.push_back(jetI);}
916 if (dist < jetB->NN_dist) {
918 jetB->NN_dist = dist;
926 while (jets_for_minheap.size() > 0) {
927 TiledJet * jetI = jets_for_minheap.back();
928 jets_for_minheap.pop_back();
929 minheap.update(jetI-head, _bj_diJ(jetI));
930 jetI->label_minheap_update_done();
940 FASTJET_END_NAMESPACE