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