FastJet  3.3.0
ClusterSequence_TiledN2.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence_TiledN2.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 
32 // The tiled N^2 part of the ClusterSequence class -- separated out
33 // from the rest of the class implementation so as to speed up
34 // compilation of this particular part while it is under test.
35 
36 #include<iostream>
37 #include<vector>
38 #include<cmath>
39 #include<algorithm>
40 #include "fastjet/PseudoJet.hh"
41 #include "fastjet/ClusterSequence.hh"
42 #include "fastjet/internal/MinHeap.hh"
43 #include "fastjet/internal/TilingExtent.hh"
44 
45 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
46 
47 using namespace std;
48 
49 
50 //----------------------------------------------------------------------
51 void ClusterSequence::_bj_remove_from_tiles(TiledJet * const jet) {
52  Tile * tile = & _tiles[jet->tile_index];
53 
54  if (jet->previous == NULL) {
55  // we are at head of the tile, so reset it.
56  // If this was the only jet on the tile then tile->head will now be NULL
57  tile->head = jet->next;
58  } else {
59  // adjust link from previous jet in this tile
60  jet->previous->next = jet->next;
61  }
62  if (jet->next != NULL) {
63  // adjust backwards-link from next jet in this tile
64  jet->next->previous = jet->previous;
65  }
66 }
67 
68 //----------------------------------------------------------------------
69 /// Set up the tiles:
70 /// - decide the range in eta
71 /// - allocate the tiles
72 /// - set up the cross-referencing info between tiles
73 ///
74 /// The neighbourhood of a tile is set up as follows
75 ///
76 /// LRR
77 /// LXR
78 /// LLR
79 ///
80 /// such that tiles is an array containing XLLLLRRRR with pointers
81 /// | \ RH_tiles
82 /// \ surrounding_tiles
83 ///
84 /// with appropriate precautions when close to the edge of the tiled
85 /// region.
86 ///
87 void ClusterSequence::_initialise_tiles() {
88 
89  // first decide tile sizes (with a lower bound to avoid huge memory use with
90  // very small R)
91  double default_size = max(0.1,_Rparam);
92  _tile_size_eta = default_size;
93  // it makes no sense to go below 3 tiles in phi -- 3 tiles is
94  // sufficient to make sure all pair-wise combinations up to pi in
95  // phi are possible
96  _n_tiles_phi = max(3,int(floor(twopi/default_size)));
97  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
98 
99  TilingExtent tiling_analysis(*this);
100  _tiles_eta_min = tiling_analysis.minrap();
101  _tiles_eta_max = tiling_analysis.maxrap();
102 
103  // // always include zero rapidity in the tiling region
104  // _tiles_eta_min = 0.0;
105  // _tiles_eta_max = 0.0;
106  // // but go no further than following
107  // const double maxrap = 7.0;
108  //
109  // // and find out how much further one should go
110  // for(unsigned int i = 0; i < _jets.size(); i++) {
111  // double eta = _jets[i].rap();
112  // // first check if eta is in range -- to avoid taking into account
113  // // very spurious rapidities due to particles with near-zero kt.
114  // if (abs(eta) < maxrap) {
115  // if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
116  // if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
117  // }
118  // }
119 
120  // now adjust the values
121  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
122  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
123  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
124  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
125 
126  // allocate the tiles
127  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
128 
129  // now set up the cross-referencing between tiles
130  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
131  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
132  Tile * tile = & _tiles[_tile_index(ieta,iphi)];
133  // no jets in this tile yet
134  tile->head = NULL; // first element of tiles points to itself
135  tile->begin_tiles[0] = tile;
136  Tile ** pptile = & (tile->begin_tiles[0]);
137  pptile++;
138  //
139  // set up L's in column to the left of X
140  tile->surrounding_tiles = pptile;
141  if (ieta > _tiles_ieta_min) {
142  // with the itile subroutine, we can safely run tiles from
143  // idphi=-1 to idphi=+1, because it takes care of
144  // negative and positive boundaries
145  for (int idphi = -1; idphi <=+1; idphi++) {
146  *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
147  pptile++;
148  }
149  }
150  // now set up last L (below X)
151  *pptile = & _tiles[_tile_index(ieta,iphi-1)];
152  pptile++;
153  // set up first R (above X)
154  tile->RH_tiles = pptile;
155  *pptile = & _tiles[_tile_index(ieta,iphi+1)];
156  pptile++;
157  // set up remaining R's, to the right of X
158  if (ieta < _tiles_ieta_max) {
159  for (int idphi = -1; idphi <= +1; idphi++) {
160  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
161  pptile++;
162  }
163  }
164  // now put semaphore for end tile
165  tile->end_tiles = pptile;
166  // finally make sure tiles are untagged
167  tile->tagged = false;
168  }
169  }
170 
171 }
172 
173 
174 //----------------------------------------------------------------------
175 /// return the tile index corresponding to the given eta,phi point
176 int ClusterSequence::_tile_index(const double eta, const double phi) const {
177  int ieta, iphi;
178  if (eta <= _tiles_eta_min) {ieta = 0;}
179  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
180  else {
181  //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
182  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
183  // following needed in case of rare but nasty rounding errors
184  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
185  ieta = _tiles_ieta_max-_tiles_ieta_min;}
186  }
187  // allow for some extent of being beyond range in calculation of phi
188  // as well
189  //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
190  // with just int and no floor, things run faster but beware
191  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
192  return (iphi + ieta * _n_tiles_phi);
193 }
194 
195 
196 //----------------------------------------------------------------------
197 // overloaded version which additionally sets up information regarding the
198 // tiling
199 inline void ClusterSequence::_tj_set_jetinfo( TiledJet * const jet,
200  const int _jets_index) {
201  // first call the generic setup
202  _bj_set_jetinfo<>(jet, _jets_index);
203 
204  // Then do the setup specific to the tiled case.
205 
206  // Find out which tile it belonds to
207  jet->tile_index = _tile_index(jet->eta, jet->phi);
208 
209  // Insert it into the tile's linked list of jets
210  Tile * tile = &_tiles[jet->tile_index];
211  jet->previous = NULL;
212  jet->next = tile->head;
213  if (jet->next != NULL) {jet->next->previous = jet;}
214  tile->head = jet;
215 }
216 
217 
218 //----------------------------------------------------------------------
219 /// output the contents of the tiles
220 void ClusterSequence::_print_tiles(TiledJet * briefjets ) const {
221  for (vector<Tile>::const_iterator tile = _tiles.begin();
222  tile < _tiles.end(); tile++) {
223  cout << "Tile " << tile - _tiles.begin()<<" = ";
224  vector<int> list;
225  for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
226  list.push_back(jetI-briefjets);
227  //cout <<" "<<jetI-briefjets;
228  }
229  sort(list.begin(),list.end());
230  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
231  cout <<"\n";
232  }
233 }
234 
235 
236 //----------------------------------------------------------------------
237 /// Add to the vector tile_union the tiles that are in the neighbourhood
238 /// of the specified tile_index, including itself -- start adding
239 /// from position n_near_tiles-1, and increase n_near_tiles as
240 /// you go along (could have done it more C++ like with vector with reserved
241 /// space, but fear is that it would have been slower, e.g. checking
242 /// for end of vector at each stage to decide whether to resize it)
243 void ClusterSequence::_add_neighbours_to_tile_union(const int tile_index,
244  vector<int> & tile_union, int & n_near_tiles) const {
245  for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
246  near_tile != _tiles[tile_index].end_tiles; near_tile++){
247  // get the tile number
248  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
249  n_near_tiles++;
250  }
251 }
252 
253 
254 //----------------------------------------------------------------------
255 /// Like _add_neighbours_to_tile_union, but only adds neighbours if
256 /// their "tagged" status is false; when a neighbour is added its
257 /// tagged status is set to true.
258 ///
259 /// Note that with a high level of warnings (-pedantic -Wextra -ansi,
260 /// gcc complains about tile_index maybe being used uninitialised for
261 /// oldB in ClusterSequence::_minheap_faster_tiled_N2_cluster(). We
262 /// have explicitly checked that it was harmless so we could disable
263 /// the gcc warning by hand using the construct below
264 ///
265 /// #pragma GCC diagnostic push
266 /// #pragma GCC diagnostic ignored "-Wpragmas"
267 /// #pragma GCC diagnostic ignored "-Wuninitialized"
268 /// #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
269 /// ...
270 /// #pragma GCC diagnostic pop
271 ///
272 /// the @GCC diagnostic push/pop directive was only introduced in
273 /// gcc-4.6, so for broader usage, we'd need to insert #pragma GCC
274 /// diagnostic ignored "-Wpragmas" at the top of this file
275 inline void ClusterSequence::_add_untagged_neighbours_to_tile_union(
276  const int tile_index,
277  vector<int> & tile_union, int & n_near_tiles) {
278  for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
279  near_tile != _tiles[tile_index].end_tiles; near_tile++){
280  if (! (*near_tile)->tagged) {
281  (*near_tile)->tagged = true;
282  // get the tile number
283  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
284  n_near_tiles++;
285  }
286  }
287 }
288 
289 
290 //----------------------------------------------------------------------
291 /// run a tiled clustering
292 void ClusterSequence::_tiled_N2_cluster() {
293 
294  _initialise_tiles();
295 
296  int n = _jets.size();
297  TiledJet * briefjets = new TiledJet[n];
298  TiledJet * jetA = briefjets, * jetB;
299  TiledJet oldB;
300  oldB.tile_index=0; // prevents a gcc warning
301 
302  // will be used quite deep inside loops, but declare it here so that
303  // memory (de)allocation gets done only once
304  vector<int> tile_union(3*n_tile_neighbours);
305 
306  // initialise the basic jet info
307  for (int i = 0; i< n; i++) {
308  _tj_set_jetinfo(jetA, i);
309  //cout << i<<": "<<jetA->tile_index<<"\n";
310  jetA++; // move on to next entry of briefjets
311  }
312  TiledJet * tail = jetA; // a semaphore for the end of briefjets
313  TiledJet * head = briefjets; // a nicer way of naming start
314 
315  // set up the initial nearest neighbour information
316  vector<Tile>::const_iterator tile;
317  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
318  // first do it on this tile
319  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
320  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
321  double dist = _bj_dist(jetA,jetB);
322  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
323  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
324  }
325  }
326  // then do it for RH tiles
327  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
328  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
329  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
330  double dist = _bj_dist(jetA,jetB);
331  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
332  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
333  }
334  }
335  }
336  }
337 
338  // now create the diJ (where J is i's NN) table -- remember that
339  // we differ from standard normalisation here by a factor of R2
340  double * diJ = new double[n];
341  jetA = head;
342  for (int i = 0; i < n; i++) {
343  diJ[i] = _bj_diJ(jetA);
344  jetA++; // have jetA follow i
345  }
346 
347  // now run the recombination loop
348  int history_location = n-1;
349  while (tail != head) {
350 
351  // find the minimum of the diJ on this round
352  double diJ_min = diJ[0];
353  int diJ_min_jet = 0;
354  for (int i = 1; i < n; i++) {
355  if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
356  }
357 
358  // do the recombination between A and B
359  history_location++;
360  jetA = & briefjets[diJ_min_jet];
361  jetB = jetA->NN;
362  // put the normalisation back in
363  diJ_min *= _invR2;
364 
365  //if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";}
366 
367  //cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n";
368 
369  if (jetB != NULL) {
370  // jet-jet recombination
371  // If necessary relabel A & B to ensure jetB < jetA, that way if
372  // the larger of them == newtail then that ends up being jetA and
373  // the new jet that is added as jetB is inserted in a position that
374  // has a future!
375  if (jetA < jetB) {std::swap(jetA,jetB);}
376 
377  int nn; // new jet index
378  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
379 
380  //OBS// get the two history indices
381  //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
382  //OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index();
383  //OBS// create the recombined jet
384  //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
385  //OBSint nn = _jets.size() - 1;
386  //OBS_jets[nn].set_cluster_hist_index(history_location);
387  //OBS// update history
388  //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
389  //OBS_add_step_to_history(history_location,
390  //OBS min(hist_a,hist_b),max(hist_a,hist_b),
391  //OBS nn, diJ_min);
392 
393  // what was jetB will now become the new jet
394  _bj_remove_from_tiles(jetA);
395  oldB = * jetB; // take a copy because we will need it...
396  _bj_remove_from_tiles(jetB);
397  _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
398  } else {
399  // jet-beam recombination
400  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
401 
402  //OBS// get the hist_index
403  //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
404  //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
405  //OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min);
406  _bj_remove_from_tiles(jetA);
407  }
408 
409  // first establish the set of tiles over which we are going to
410  // have to run searches for updated and new nearest-neighbours
411  int n_near_tiles = 0;
412  _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
413  if (jetB != NULL) {
414  bool sort_it = false;
415  if (jetB->tile_index != jetA->tile_index) {
416  sort_it = true;
417  _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
418  }
419  if (oldB.tile_index != jetA->tile_index &&
420  oldB.tile_index != jetB->tile_index) {
421  sort_it = true;
422  _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
423  }
424 
425  if (sort_it) {
426  // sort the tiles before then compressing the list
427  sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
428  // and now condense the list
429  int nnn = 1;
430  for (int i = 1; i < n_near_tiles; i++) {
431  if (tile_union[i] != tile_union[nnn-1]) {
432  tile_union[nnn] = tile_union[i];
433  nnn++;
434  }
435  }
436  n_near_tiles = nnn;
437  }
438  }
439 
440  // now update our nearest neighbour info and diJ table
441  // first reduce size of table
442  tail--; n--;
443  if (jetA == tail) {
444  // there is nothing to be done
445  } else {
446  // Copy last jet contents and diJ info into position of jetA
447  *jetA = *tail;
448  diJ[jetA - head] = diJ[tail-head];
449  // IN the tiling fix pointers to tail and turn them into
450  // pointers to jetA (from predecessors, successors and the tile
451  // head if need be)
452  if (jetA->previous == NULL) {
453  _tiles[jetA->tile_index].head = jetA;
454  } else {
455  jetA->previous->next = jetA;
456  }
457  if (jetA->next != NULL) {jetA->next->previous = jetA;}
458  }
459 
460  // Initialise jetB's NN distance as well as updating it for
461  // other particles.
462  for (int itile = 0; itile < n_near_tiles; itile++) {
463  Tile * tile_ptr = &_tiles[tile_union[itile]];
464  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
465  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
466  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
467  jetI->NN_dist = _R2;
468  jetI->NN = NULL;
469  // now go over tiles that are neighbours of I (include own tile)
470  for (Tile ** near_tile = tile_ptr->begin_tiles;
471  near_tile != tile_ptr->end_tiles; near_tile++) {
472  // and then over the contents of that tile
473  for (TiledJet * jetJ = (*near_tile)->head;
474  jetJ != NULL; jetJ = jetJ->next) {
475  double dist = _bj_dist(jetI,jetJ);
476  if (dist < jetI->NN_dist && jetJ != jetI) {
477  jetI->NN_dist = dist; jetI->NN = jetJ;
478  }
479  }
480  }
481  diJ[jetI-head] = _bj_diJ(jetI); // update diJ
482  }
483  // check whether new jetB is closer than jetI's current NN and
484  // if need to update things
485  if (jetB != NULL) {
486  double dist = _bj_dist(jetI,jetB);
487  if (dist < jetI->NN_dist) {
488  if (jetI != jetB) {
489  jetI->NN_dist = dist;
490  jetI->NN = jetB;
491  diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
492  }
493  }
494  if (dist < jetB->NN_dist) {
495  if (jetI != jetB) {
496  jetB->NN_dist = dist;
497  jetB->NN = jetI;}
498  }
499  }
500  }
501  }
502 
503 
504  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
505  //cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
506 
507  // remember to update pointers to tail
508  for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
509  near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
510  // and then the contents of that tile
511  for (TiledJet * jetJ = (*near_tile)->head;
512  jetJ != NULL; jetJ = jetJ->next) {
513  if (jetJ->NN == tail) {jetJ->NN = jetA;}
514  }
515  }
516 
517  //for (int i = 0; i < n; i++) {
518  // 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";}
519  //}
520 
521 
522  if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
523  //cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
524 
525  }
526 
527  // final cleaning up;
528  delete[] diJ;
529  delete[] briefjets;
530 }
531 
532 
533 //----------------------------------------------------------------------
534 /// run a tiled clustering
535 void ClusterSequence::_faster_tiled_N2_cluster() {
536 
537  _initialise_tiles();
538 
539  int n = _jets.size();
540  TiledJet * briefjets = new TiledJet[n];
541  TiledJet * jetA = briefjets, * jetB;
542  TiledJet oldB;
543  oldB.tile_index=0; // prevents a gcc warning
544 
545  // will be used quite deep inside loops, but declare it here so that
546  // memory (de)allocation gets done only once
547  vector<int> tile_union(3*n_tile_neighbours);
548 
549  // initialise the basic jet info
550  for (int i = 0; i< n; i++) {
551  _tj_set_jetinfo(jetA, i);
552  //cout << i<<": "<<jetA->tile_index<<"\n";
553  jetA++; // move on to next entry of briefjets
554  }
555  TiledJet * head = briefjets; // a nicer way of naming start
556 
557  // set up the initial nearest neighbour information
558  vector<Tile>::const_iterator tile;
559  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
560  // first do it on this tile
561  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
562  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
563  double dist = _bj_dist(jetA,jetB);
564  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
565  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
566  }
567  }
568  // then do it for RH tiles
569  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
570  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
571  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
572  double dist = _bj_dist(jetA,jetB);
573  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
574  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
575  }
576  }
577  }
578  // no need to do it for LH tiles, since they are implicitly done
579  // when we set NN for both jetA and jetB on the RH tiles.
580  }
581 
582  // now create the diJ (where J is i's NN) table -- remember that
583  // we differ from standard normalisation here by a factor of R2
584  // (corrected for at the end).
585  struct diJ_plus_link {
586  double diJ; // the distance
587  TiledJet * jet; // the jet (i) for which we've found this distance
588  // (whose NN will the J).
589  };
590  diJ_plus_link * diJ = new diJ_plus_link[n];
591  jetA = head;
592  for (int i = 0; i < n; i++) {
593  diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
594  diJ[i].jet = jetA; // our compact diJ table will not be in
595  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
596  // so set up bi-directional correspondence here.
597  jetA++; // have jetA follow i
598  }
599 
600  // now run the recombination loop
601  int history_location = n-1;
602  while (n > 0) {
603 
604  // find the minimum of the diJ on this round
605  diJ_plus_link * best, *stop; // pointers a bit faster than indices
606  // could use best to keep track of diJ min, but it turns out to be
607  // marginally faster to have a separate variable (avoids n
608  // dereferences at the expense of n/2 assignments).
609  double diJ_min = diJ[0].diJ; // initialise the best one here.
610  best = diJ; // and here
611  stop = diJ+n;
612  for (diJ_plus_link * here = diJ+1; here != stop; here++) {
613  if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
614  }
615 
616  // do the recombination between A and B
617  history_location++;
618  jetA = best->jet;
619  jetB = jetA->NN;
620  // put the normalisation back in
621  diJ_min *= _invR2;
622 
623  if (jetB != NULL) {
624  // jet-jet recombination
625  // If necessary relabel A & B to ensure jetB < jetA, that way if
626  // the larger of them == newtail then that ends up being jetA and
627  // the new jet that is added as jetB is inserted in a position that
628  // has a future!
629  if (jetA < jetB) {std::swap(jetA,jetB);}
630 
631  int nn; // new jet index
632  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
633 
634  //OBS// get the two history indices
635  //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
636  //OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index();
637  //OBS// create the recombined jet
638  //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
639  //OBSint nn = _jets.size() - 1;
640  //OBS_jets[nn].set_cluster_hist_index(history_location);
641  //OBS// update history
642  //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
643  //OBS_add_step_to_history(history_location,
644  //OBS min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b),
645  //OBS nn, diJ_min);
646  // what was jetB will now become the new jet
647  _bj_remove_from_tiles(jetA);
648  oldB = * jetB; // take a copy because we will need it...
649  _bj_remove_from_tiles(jetB);
650  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
651  // (also registers the jet in the tiling)
652  } else {
653  // jet-beam recombination
654  // get the hist_index
655  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
656  //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
657  //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
658  //OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min);
659  _bj_remove_from_tiles(jetA);
660  }
661 
662  // first establish the set of tiles over which we are going to
663  // have to run searches for updated and new nearest-neighbours --
664  // basically a combination of vicinity of the tiles of the two old
665  // and one new jet.
666  int n_near_tiles = 0;
667  _add_untagged_neighbours_to_tile_union(jetA->tile_index,
668  tile_union, n_near_tiles);
669  if (jetB != NULL) {
670  if (jetB->tile_index != jetA->tile_index) {
671  _add_untagged_neighbours_to_tile_union(jetB->tile_index,
672  tile_union,n_near_tiles);
673  }
674  if (oldB.tile_index != jetA->tile_index &&
675  oldB.tile_index != jetB->tile_index) {
676  _add_untagged_neighbours_to_tile_union(oldB.tile_index,
677  tile_union,n_near_tiles);
678  }
679  }
680 
681  // now update our nearest neighbour info and diJ table
682  // first reduce size of table
683  n--;
684  // then compactify the diJ by taking the last of the diJ and copying
685  // it to the position occupied by the diJ for jetA
686  diJ[n].jet->diJ_posn = jetA->diJ_posn;
687  diJ[jetA->diJ_posn] = diJ[n];
688 
689  // Initialise jetB's NN distance as well as updating it for
690  // other particles.
691  // Run over all tiles in our union
692  for (int itile = 0; itile < n_near_tiles; itile++) {
693  Tile * tile_ptr = &_tiles[tile_union[itile]];
694  tile_ptr->tagged = false; // reset tag, since we're done with unions
695  // run over all jets in the current tile
696  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
697  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
698  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
699  jetI->NN_dist = _R2;
700  jetI->NN = NULL;
701  // now go over tiles that are neighbours of I (include own tile)
702  for (Tile ** near_tile = tile_ptr->begin_tiles;
703  near_tile != tile_ptr->end_tiles; near_tile++) {
704  // and then over the contents of that tile
705  for (TiledJet * jetJ = (*near_tile)->head;
706  jetJ != NULL; jetJ = jetJ->next) {
707  double dist = _bj_dist(jetI,jetJ);
708  if (dist < jetI->NN_dist && jetJ != jetI) {
709  jetI->NN_dist = dist; jetI->NN = jetJ;
710  }
711  }
712  }
713  diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
714  }
715  // check whether new jetB is closer than jetI's current NN and
716  // if jetI is closer than jetB's current (evolving) nearest
717  // neighbour. Where relevant update things
718  if (jetB != NULL) {
719  double dist = _bj_dist(jetI,jetB);
720  if (dist < jetI->NN_dist) {
721  if (jetI != jetB) {
722  jetI->NN_dist = dist;
723  jetI->NN = jetB;
724  diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
725  }
726  }
727  if (dist < jetB->NN_dist) {
728  if (jetI != jetB) {
729  jetB->NN_dist = dist;
730  jetB->NN = jetI;}
731  }
732  }
733  }
734  }
735 
736  // finally, register the updated kt distance for B
737  if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
738 
739  }
740 
741  // final cleaning up;
742  delete[] diJ;
743  delete[] briefjets;
744 }
745 
746 //----------------------------------------------------------------------
747 /// run a tiled clustering, with our minheap for keeping track of the
748 /// smallest dij
749 void ClusterSequence::_minheap_faster_tiled_N2_cluster() {
750 
751  _initialise_tiles();
752 
753  int n = _jets.size();
754  TiledJet * briefjets = new TiledJet[n];
755  TiledJet * jetA = briefjets, * jetB;
756  TiledJet oldB;
757  oldB.tile_index=0; // prevents a gcc warning
758 
759  // will be used quite deep inside loops, but declare it here so that
760  // memory (de)allocation gets done only once
761  vector<int> tile_union(3*n_tile_neighbours);
762 
763  // initialise the basic jet info
764  for (int i = 0; i< n; i++) {
765  _tj_set_jetinfo(jetA, i);
766  //cout << i<<": "<<jetA->tile_index<<"\n";
767  jetA++; // move on to next entry of briefjets
768  }
769  TiledJet * head = briefjets; // a nicer way of naming start
770 
771  // set up the initial nearest neighbour information
772  vector<Tile>::const_iterator tile;
773  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
774  // first do it on this tile
775  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
776  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
777  double dist = _bj_dist(jetA,jetB);
778  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
779  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
780  }
781  }
782  // then do it for RH tiles
783  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
784  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
785  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
786  double dist = _bj_dist(jetA,jetB);
787  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
788  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
789  }
790  }
791  }
792  // no need to do it for LH tiles, since they are implicitly done
793  // when we set NN for both jetA and jetB on the RH tiles.
794  }
795 
796 
797  //// now create the diJ (where J is i's NN) table -- remember that
798  //// we differ from standard normalisation here by a factor of R2
799  //// (corrected for at the end).
800  //struct diJ_plus_link {
801  // double diJ; // the distance
802  // TiledJet * jet; // the jet (i) for which we've found this distance
803  // // (whose NN will the J).
804  //};
805  //diJ_plus_link * diJ = new diJ_plus_link[n];
806  //jetA = head;
807  //for (int i = 0; i < n; i++) {
808  // diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
809  // diJ[i].jet = jetA; // our compact diJ table will not be in
810  // jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
811  // // so set up bi-directional correspondence here.
812  // jetA++; // have jetA follow i
813  //}
814 
815  vector<double> diJs(n);
816  for (int i = 0; i < n; i++) {
817  diJs[i] = _bj_diJ(&briefjets[i]);
818  briefjets[i].label_minheap_update_done();
819  }
820  MinHeap minheap(diJs);
821  // have a stack telling us which jets we'll have to update on the heap
822  vector<TiledJet *> jets_for_minheap;
823  jets_for_minheap.reserve(n);
824 
825  // now run the recombination loop
826  int history_location = n-1;
827  while (n > 0) {
828 
829  double diJ_min = minheap.minval() *_invR2;
830  jetA = head + minheap.minloc();
831 
832  // do the recombination between A and B
833  history_location++;
834  jetB = jetA->NN;
835 
836  if (jetB != NULL) {
837  // jet-jet recombination
838  // If necessary relabel A & B to ensure jetB < jetA, that way if
839  // the larger of them == newtail then that ends up being jetA and
840  // the new jet that is added as jetB is inserted in a position that
841  // has a future!
842  if (jetA < jetB) {std::swap(jetA,jetB);}
843 
844  int nn; // new jet index
845  _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
846 
847  // what was jetB will now become the new jet
848  _bj_remove_from_tiles(jetA);
849  oldB = * jetB; // take a copy because we will need it...
850  _bj_remove_from_tiles(jetB);
851  _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
852  // (also registers the jet in the tiling)
853  } else {
854  // jet-beam recombination
855  // get the hist_index
856  _do_iB_recombination_step(jetA->_jets_index, diJ_min);
857  _bj_remove_from_tiles(jetA);
858  }
859 
860  // remove the minheap entry for jetA
861  minheap.remove(jetA-head);
862 
863  // first establish the set of tiles over which we are going to
864  // have to run searches for updated and new nearest-neighbours --
865  // basically a combination of vicinity of the tiles of the two old
866  // and one new jet.
867  int n_near_tiles = 0;
868  _add_untagged_neighbours_to_tile_union(jetA->tile_index,
869  tile_union, n_near_tiles);
870  if (jetB != NULL) {
871  if (jetB->tile_index != jetA->tile_index) {
872  _add_untagged_neighbours_to_tile_union(jetB->tile_index,
873  tile_union,n_near_tiles);
874  }
875  if (oldB.tile_index != jetA->tile_index &&
876  oldB.tile_index != jetB->tile_index) {
877  // GS: the line below generates a warning that oldB.tile_index
878  // may be used uninitialised. However, to reach this point, we
879  // ned jetB != NULL (see test a few lines above) and is jetB
880  // !=NULL, one would have gone through "oldB = *jetB before
881  // (see piece of code ~20 line above), so the index is
882  // initialised. We do not do anything to avoid the warning to
883  // avoid any potential speed impact.
884  _add_untagged_neighbours_to_tile_union(oldB.tile_index,
885  tile_union,n_near_tiles);
886  }
887  // indicate that we'll have to update jetB in the minheap
888  jetB->label_minheap_update_needed();
889  jets_for_minheap.push_back(jetB);
890  }
891 
892 
893  // Initialise jetB's NN distance as well as updating it for
894  // other particles.
895  // Run over all tiles in our union
896  for (int itile = 0; itile < n_near_tiles; itile++) {
897  Tile * tile_ptr = &_tiles[tile_union[itile]];
898  tile_ptr->tagged = false; // reset tag, since we're done with unions
899  // run over all jets in the current tile
900  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
901  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
902  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
903  jetI->NN_dist = _R2;
904  jetI->NN = NULL;
905  // label jetI as needing heap action...
906  if (!jetI->minheap_update_needed()) {
907  jetI->label_minheap_update_needed();
908  jets_for_minheap.push_back(jetI);}
909  // now go over tiles that are neighbours of I (include own tile)
910  for (Tile ** near_tile = tile_ptr->begin_tiles;
911  near_tile != tile_ptr->end_tiles; near_tile++) {
912  // and then over the contents of that tile
913  for (TiledJet * jetJ = (*near_tile)->head;
914  jetJ != NULL; jetJ = jetJ->next) {
915  double dist = _bj_dist(jetI,jetJ);
916  if (dist < jetI->NN_dist && jetJ != jetI) {
917  jetI->NN_dist = dist; jetI->NN = jetJ;
918  }
919  }
920  }
921  }
922  // check whether new jetB is closer than jetI's current NN and
923  // if jetI is closer than jetB's current (evolving) nearest
924  // neighbour. Where relevant update things
925  if (jetB != NULL) {
926  double dist = _bj_dist(jetI,jetB);
927  if (dist < jetI->NN_dist) {
928  if (jetI != jetB) {
929  jetI->NN_dist = dist;
930  jetI->NN = jetB;
931  // label jetI as needing heap action...
932  if (!jetI->minheap_update_needed()) {
933  jetI->label_minheap_update_needed();
934  jets_for_minheap.push_back(jetI);}
935  }
936  }
937  if (dist < jetB->NN_dist) {
938  if (jetI != jetB) {
939  jetB->NN_dist = dist;
940  jetB->NN = jetI;}
941  }
942  }
943  }
944  }
945 
946  // deal with jets whose minheap entry needs updating
947  while (jets_for_minheap.size() > 0) {
948  TiledJet * jetI = jets_for_minheap.back();
949  jets_for_minheap.pop_back();
950  minheap.update(jetI-head, _bj_diJ(jetI));
951  jetI->label_minheap_update_done();
952  }
953  n--;
954  }
955 
956  // final cleaning up;
957  delete[] briefjets;
958 }
959 
960 
961 FASTJET_END_NAMESPACE
962