FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LazyTiling9SeparateGhosts.cc
1 //FJSTARTHEADER
2 // $Id: LazyTiling9SeparateGhosts.cc 3596 2014-08-12 15:27:19Z soyez $
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 #include "fastjet/internal/LazyTiling9SeparateGhosts.hh"
32 #include "fastjet/internal/TilingExtent.hh"
33 #include <iomanip>
34 using namespace std;
35 
36 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37 
38 double LazyTiling9SeparateGhosts::ghost_pt2_threshold = 1e-100;
39 
40 LazyTiling9SeparateGhosts::LazyTiling9SeparateGhosts(ClusterSequence & cs) :
41  _cs(cs), _jets(cs.jets())
42  //, _minheap(_jets.size())
43 {
44  _Rparam = cs.jet_def().R();
45  _R2 = _Rparam * _Rparam;
46  _invR2 = 1.0 / _R2;
47  _initialise_tiles();
48 }
49 
50 
51 //----------------------------------------------------------------------
52 /// Set up the tiles:
53 /// - decide the range in eta
54 /// - allocate the tiles
55 /// - set up the cross-referencing info between tiles
56 ///
57 /// The neighbourhood of a tile is set up as follows
58 ///
59 /// LRR
60 /// LXR
61 /// LLR
62 ///
63 /// such that tiles is an array containing XLLLLRRRR with pointers
64 /// | \ RH_tiles
65 /// \ surrounding_tiles
66 ///
67 /// with appropriate precautions when close to the edge of the tiled
68 /// region.
69 ///
70 void LazyTiling9SeparateGhosts::_initialise_tiles() {
71 
72  // first decide tile sizes (with a lower bound to avoid huge memory use with
73  // very small R)
74  double default_size = max(0.1,_Rparam);
75  _tile_size_eta = default_size;
76  // it makes no sense to go below 3 tiles in phi -- 3 tiles is
77  // sufficient to make sure all pair-wise combinations up to pi in
78  // phi are possible
79  _n_tiles_phi = max(3,int(floor(twopi/default_size)));
80  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
81 
82  // always include zero rapidity in the tiling region
83  _tiles_eta_min = 0.0;
84  _tiles_eta_max = 0.0;
85  // but go no further than following
86  const double maxrap = 7.0;
87 
88  // and find out how much further one should go
89  for(unsigned int i = 0; i < _jets.size(); i++) {
90  double eta = _jets[i].rap();
91  // first check if eta is in range -- to avoid taking into account
92  // very spurious rapidities due to particles with near-zero kt.
93  if (abs(eta) < maxrap) {
94  if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
95  if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
96  }
97  }
98 
99  // now adjust the values
100  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
101  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
102  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
103  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
104 
105  _tile_half_size_eta = _tile_size_eta * 0.5;
106  _tile_half_size_phi = _tile_size_phi * 0.5;
107 
108  // allocate the tiles
109  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
110 
111  // now set up the cross-referencing between tiles
112  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
113  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
114  Tile3 * tile = & _tiles[_tile_index(ieta,iphi)];
115  // no jets in this tile yet
116  tile->head = NULL; // first element of tiles points to itself
117  tile->ghost_head = NULL; // first element of tiles points to itself
118  tile->begin_tiles[0] = tile;
119  Tile3 ** pptile = & (tile->begin_tiles[0]);
120  pptile++;
121  //
122  // set up L's in column to the left of X
123  tile->surrounding_tiles = pptile;
124  if (ieta > _tiles_ieta_min) {
125  // with the itile subroutine, we can safely run tiles from
126  // idphi=-1 to idphi=+1, because it takes care of
127  // negative and positive boundaries
128  for (int idphi = -1; idphi <=+1; idphi++) {
129  *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
130  pptile++;
131  }
132  }
133  // now set up last L (below X)
134  *pptile = & _tiles[_tile_index(ieta,iphi-1)];
135  pptile++;
136  // set up first R (above X)
137  tile->RH_tiles = pptile;
138  *pptile = & _tiles[_tile_index(ieta,iphi+1)];
139  pptile++;
140  // set up remaining R's, to the right of X
141  if (ieta < _tiles_ieta_max) {
142  for (int idphi = -1; idphi <= +1; idphi++) {
143  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
144  pptile++;
145  }
146  }
147  // now put semaphore for end tile
148  tile->end_tiles = pptile;
149  // finally make sure tiles are untagged
150  tile->tagged = false;
151  // and ensure max distance is sensibly initialised
152  tile->max_NN_dist = 0;
153  // and also position of centre of tile
154  tile->eta_centre = (ieta+0.5)*_tile_size_eta;
155  tile->phi_centre = (iphi+0.5)*_tile_size_phi;
156  }
157  }
158 
159 }
160 
161 //----------------------------------------------------------------------
162 /// return the tile index corresponding to the given eta,phi point
163 int LazyTiling9SeparateGhosts::_tile_index(const double eta, const double phi) const {
164  int ieta, iphi;
165  if (eta <= _tiles_eta_min) {ieta = 0;}
166  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
167  else {
168  //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
169  ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
170  // following needed in case of rare but nasty rounding errors
171  if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
172  ieta = _tiles_ieta_max-_tiles_ieta_min;}
173  }
174  // allow for some extent of being beyond range in calculation of phi
175  // as well
176  //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
177  // with just int and no floor, things run faster but beware
178  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
179  return (iphi + ieta * _n_tiles_phi);
180 }
181 
182 
183 //----------------------------------------------------------------------
184 // sets up information regarding the tiling of the given jet
185 inline void LazyTiling9SeparateGhosts::_tj_set_jetinfo( TiledJet3 * const jet,
186  const int _jets_index, bool is_ghost) {
187  // first call the generic setup
188  _bj_set_jetinfo<>(jet, _jets_index);
189 
190  // Then do the setup specific to the tiled case.
191  jet->is_ghost = is_ghost;
192 
193  // Find out which tile it belonds to
194  jet->tile_index = _tile_index(jet->eta, jet->phi);
195 
196  // Insert it into the tile's linked list of jets
197  Tile3 * tile = &_tiles[jet->tile_index];
198  jet->previous = NULL;
199  if (is_ghost) {
200  jet->next = tile->ghost_head;
201  tile->ghost_head = jet;
202  } else {
203  jet->next = tile->head;
204  tile->head = jet;
205  }
206  if (jet->next != NULL) {jet->next->previous = jet;}
207 }
208 
209 
210 //----------------------------------------------------------------------
211 void LazyTiling9SeparateGhosts::_bj_remove_from_tiles(TiledJet3 * const jet) {
212  Tile3 * tile = & _tiles[jet->tile_index];
213 
214  if (jet->previous == NULL) {
215  // we are at head of the tile, so reset it.
216  // If this was the only jet on the tile then tile->head will now be NULL
217  if (jet->is_ghost) {
218  tile->ghost_head = jet->next;
219  } else {
220  tile->head = jet->next;
221  }
222  } else {
223  // adjust link from previous jet in this tile
224  jet->previous->next = jet->next;
225  }
226  if (jet->next != NULL) {
227  // adjust backwards-link from next jet in this tile
228  jet->next->previous = jet->previous;
229  }
230 }
231 
232 
233 //----------------------------------------------------------------------
234 /// output the contents of the tiles
235 void LazyTiling9SeparateGhosts::_print_tiles(TiledJet3 * briefjets ) const {
236  for (vector<Tile3>::const_iterator tile = _tiles.begin();
237  tile < _tiles.end(); tile++) {
238  cout << "Tile " << tile - _tiles.begin()<<" = ";
239  vector<int> list;
240  for (TiledJet3 * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
241  list.push_back(jetI-briefjets);
242  //cout <<" "<<jetI-briefjets;
243  }
244  sort(list.begin(),list.end());
245  for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
246  cout <<"\n";
247  }
248 }
249 
250 
251 //----------------------------------------------------------------------
252 /// Add to the vector tile_union the tiles that are in the neighbourhood
253 /// of the specified tile_index, including itself -- start adding
254 /// from position n_near_tiles-1, and increase n_near_tiles as
255 /// you go along (could have done it more C++ like with vector with reserved
256 /// space, but fear is that it would have been slower, e.g. checking
257 /// for end of vector at each stage to decide whether to resize it)
258 void LazyTiling9SeparateGhosts::_add_neighbours_to_tile_union(const int tile_index,
259  vector<int> & tile_union, int & n_near_tiles) const {
260  for (Tile3 * const * near_tile = _tiles[tile_index].begin_tiles;
261  near_tile != _tiles[tile_index].end_tiles; near_tile++){
262  // get the tile number
263  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
264  n_near_tiles++;
265  }
266 }
267 
268 
269 //----------------------------------------------------------------------
270 /// Like _add_neighbours_to_tile_union, but only adds neighbours if
271 /// their "tagged" status is false; when a neighbour is added its
272 /// tagged status is set to true.
273 inline void LazyTiling9SeparateGhosts::_add_untagged_neighbours_to_tile_union(
274  const int tile_index,
275  vector<int> & tile_union, int & n_near_tiles) {
276  for (Tile3 ** near_tile = _tiles[tile_index].begin_tiles;
277  near_tile != _tiles[tile_index].end_tiles; near_tile++){
278  if (! (*near_tile)->tagged) {
279  (*near_tile)->tagged = true;
280  // get the tile number
281  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
282  n_near_tiles++;
283  }
284  }
285 }
286 
287 //----------------------------------------------------------------------
288 /// Like _add_neighbours_to_tile_union, but adds tiles that are
289 /// "neighbours" of a jet (rather than a tile) and only if a
290 /// neighbouring tile's max_NN_dist is >= the distance between the jet
291 /// and the nearest point on the tile. It ignores tiles that have
292 /// already been tagged.
293 inline void LazyTiling9SeparateGhosts::_add_untagged_neighbours_to_tile_union_using_max_info(
294  const TiledJet3 * jet,
295  vector<int> & tile_union, int & n_near_tiles) {
296  Tile3 & tile = _tiles[jet->tile_index];
297 
298  for (Tile3 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
299  if ((*near_tile)->tagged) continue;
300  double dist = _distance_to_tile(jet, *near_tile);
301  // cout << " max info looked at tile " << *near_tile - &_tiles[0]
302  // << ", dist = " << dist << " " << (*near_tile)->max_NN_dist
303  // << endl;
304  if (dist > (*near_tile)->max_NN_dist) continue;
305 
306  // cout << " max info tagged tile " << *near_tile - &_tiles[0] << endl;
307  (*near_tile)->tagged = true;
308  // get the tile number
309  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
310  n_near_tiles++;
311  }
312 }
313 
314 ////--------TMPTMPTMPTMPTMP-----GPS TEMP--------------------
315 //ostream & operator<<(ostream & ostr, const TiledJet3 & jet) {
316 // ostr << "j" << setw(3) << jet._jets_index << ":pt2,rap,phi=" ; ostr.flush();
317 // ostr << jet.kt2 << ","; ostr.flush();
318 // ostr << jet.eta << ","; ostr.flush();
319 // ostr << jet.phi; ostr.flush();
320 // ostr << ", tile=" << jet.tile_index; ostr.flush();
321 // return ostr;
322 //}
323 
324 
325 //----------------------------------------------------------------------
326 /// returns a particle's distance to the edge of the specified tile
327 inline double LazyTiling9SeparateGhosts::_distance_to_tile(const TiledJet3 * bj, const Tile3 * tile) const {
328 
329  // Note the careful way of checking the minimum potential deta:
330  // unlike the phi case below, we don't calculate the distance to the
331  // centre and subtract spacing/2. This is because of issue of
332  // boundary tiles, which can extend far beyond spacing/2 in eta.
333  // Using the positions of tile centers should instead be safe.
334  double deta;
335  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
336  //else deta = std::abs(bj->eta - tile->eta_centre) - 0.5*_tile_size_eta;
337  else deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
338  // ------
339  // |
340  // A | B
341  // ------
342  // |
343  // C | D
344  // ------
345 
346  double dphi = std::abs(bj->phi - tile->phi_centre);
347  if (dphi > pi) dphi = twopi-dphi;
348  dphi -= _tile_half_size_phi;
349  //dphi -= 0.5*_tile_size_phi;
350  if (dphi < 0) dphi = 0;
351 
352  return dphi*dphi + deta*deta;
353 }
354 
355 
356 
357 
358 //----------------------------------------------------------------------
359 /// looks at distance between jetX and jetI and updates the NN
360 /// information if relevant; also pushes identity of jetI onto
361 /// the vector of jets for minheap, to signal that it will have
362 /// to be handled later.
363 ///
364 inline void LazyTiling9SeparateGhosts::_update_jetX_jetI_NN(TiledJet3 * jetX, TiledJet3 * jetI, vector<TiledJet3 *> & jets_for_minheap) {
365  assert(! (jetX->is_ghost || jetI->is_ghost));
366  double dist = _bj_dist(jetI,jetX);
367  if (dist < jetI->NN_dist) {
368  if (jetI != jetX) {
369  jetI->NN_dist = dist;
370  jetI->NN = jetX;
371  // label jetI as needing heap action...
372  if (!jetI->minheap_update_needed()) {
373  jetI->label_minheap_update_needed();
374  jets_for_minheap.push_back(jetI);
375  }
376  }
377  }
378  if (dist < jetX->NN_dist) {
379  if (jetI != jetX) {
380  jetX->NN_dist = dist;
381  jetX->NN = jetI;}
382  }
383 }
384 
385 
386 inline void LazyTiling9SeparateGhosts::_set_NN(TiledJet3 * jetI,
387  vector<TiledJet3 *> & jets_for_minheap) {
388  assert(! jetI->is_ghost);
389  jetI->NN_dist = _R2;
390  jetI->NN = NULL;
391  // label jetI as needing heap action...
392  if (!jetI->minheap_update_needed()) {
393  jetI->label_minheap_update_needed();
394  jets_for_minheap.push_back(jetI);}
395  // now go over tiles that are neighbours of I (include own tile)
396  Tile3 * tile_ptr = &_tiles[jetI->tile_index];
397  //if (tile_ptr->is_near_zero_phi(_tile_size_phi)) {
398  for (Tile3 ** near_tile = tile_ptr->begin_tiles;
399  near_tile != tile_ptr->end_tiles; near_tile++) {
400  // for own tile, this will be zero automatically: should we be clever
401  // and skip the test? (With some doubling of code?)
402  if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
403  // and then over the contents of that tile
404  for (TiledJet3 * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
405  double dist = _bj_dist(jetI,jetJ);
406  if (dist < jetI->NN_dist && jetJ != jetI) {
407  jetI->NN_dist = dist; jetI->NN = jetJ;
408  }
409  }
410  // deal with the ghosts
411  for (TiledJet3 * jetJ = (*near_tile)->ghost_head; jetJ != NULL; jetJ = jetJ->next) {
412  double dist = _bj_dist(jetI,jetJ);
413  if (dist < jetI->NN_dist) {
414  jetI->NN_dist = dist; jetI->NN = jetJ;
415  }
416  }
417  }
418  // } else {
419  // // second copy that exploits the fact that for this tile we needn't worry
420  // // about periodicity
421  // for (Tile3 ** near_tile = tile_ptr->begin_tiles;
422  // near_tile != tile_ptr->end_tiles; near_tile++) {
423  // // for own tile, this will be zero automatically: should we be clever
424  // // and skip the test? (With some doubling of code?)
425  // if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
426  // // and then over the contents of that tile
427  // for (TiledJet3 * jetJ = (*near_tile)->head;
428  // jetJ != NULL; jetJ = jetJ->next) {
429  // double dist = _bj_dist_not_periodic(jetI,jetJ);
430  // if (dist < jetI->NN_dist && jetJ != jetI) {
431  // jetI->NN_dist = dist; jetI->NN = jetJ;
432  // }
433  // }
434  // }
435  // }
436 }
437 
438 
439 void LazyTiling9SeparateGhosts::run() {
440 
441  //_initialise_tiles();
442 
443  int ntot = _jets.size();
444  if (ntot == 0) return;
445 
446  TiledJet3 * briefjets = new TiledJet3[ntot];
447  TiledJet3 * jetA = briefjets, * jetB;
448  // avoid warning about uninitialised oldB below;
449  // only valid for ntot>=1 (hence the test ntot==0 test above)
450  TiledJet3 oldB = briefjets[0];
451 
452 
453  // will be used quite deep inside loops, but declare it here so that
454  // memory (de)allocation gets done only once
455  vector<int> tile_union(3*n_tile_neighbours);
456 
457  TiledJet3 * head = briefjets; // a nicer way of naming start
458 
459  // initialise the basic jet info
460  //
461  // Note that the threshold is a static member of the class
462  // first get the particles we'll keep as "real"
463  for (int i = 0; i< ntot; i++) {
464  bool is_ghost = _jets[i].perp2() < ghost_pt2_threshold;
465  if (!is_ghost) {
466  _tj_set_jetinfo(jetA, i, is_ghost);
467  jetA++; // move on to next entry of briefjets
468  }
469  }
470  int nreal = jetA - briefjets;
471  // then the ones we will label as ghosts
472  for (int i = 0; i< ntot; i++) {
473  bool is_ghost = _jets[i].perp2() < ghost_pt2_threshold;
474  if (is_ghost) {
475  _tj_set_jetinfo(jetA, i, is_ghost);
476  jetA++; // move on to next entry of briefjets
477  }
478  }
479 
480 
481  // set up the initial nearest neighbour information
482  vector<Tile3>::iterator tile;
483  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
484  // first do it on this tile
485  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
486  // real particles cluster with real particles
487  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
488  double dist = _bj_dist_not_periodic(jetA,jetB);
489  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
490  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
491  }
492  // they can also cluster with ghosts
493  for (jetB = tile->ghost_head; jetB != NULL; jetB = jetB->next) {
494  double dist = _bj_dist_not_periodic(jetA,jetB);
495  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
496  }
497  }
498  // only look out for NN dists of real particles, because ghosts never
499  // have NN in our structure
500  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
501  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
502  }
503  }
504  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
505  //if (tile->is_near_zero_phi(_tile_size_phi)) {
506  // then do it for RH tiles;
507  for (Tile3 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
508  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
509  double dist_to_tile = _distance_to_tile(jetA, *RTile);
510  // it only makes sense to do a tile if jetA is close enough to the Rtile
511  // either for a jet in the Rtile to be closer to jetA than it's current NN
512  // or if jetA could be closer to something in the Rtile than the largest
513  // NN distance within the RTile.
514  //
515  // GPS note: also tried approach where we perform only the
516  // first test and run over all surrounding tiles
517  // (not just RH ones). The test is passed less
518  // frequently, but one is running over more tiles
519  // and on balance, for the trial event we used, it's
520  // a bit slower.
521  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
522  bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
523  if (relevant_for_jetA || relevant_for_RTile) {
524  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
525  double dist = _bj_dist(jetA,jetB);
526  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
527  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
528  }
529  }
530  // now do the check over ghosts
531  if (relevant_for_jetA) {
532  for (jetB = (*RTile)->ghost_head; jetB != NULL; jetB = jetB->next) {
533  double dist = _bj_dist(jetA,jetB);
534  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
535  }
536  }
537  }
538  }
539  // and do a special loop to catch ghosts that are among the LH tiles
540  for (Tile3 ** LTile = tile->surrounding_tiles; LTile != tile->RH_tiles; LTile++) {
541  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
542  double dist_to_tile = _distance_to_tile(jetA, *LTile);
543  // it only makes sense to do a tile if jetA is close enough to the Rtile
544  // for a (ghost) jet in the Ltile to be closer to jetA than it's current NN.
545  bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
546  if (relevant_for_jetA) {
547  for (jetB = (*LTile)->ghost_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  }
551  }
552  }
553  }
554  // } else {
555  // // this second version of the code uses the faster
556  // // "not_periodic" version because it knows that the tile is
557  // // sufficiently far from the edge.
558  // for (Tile3 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
559  // for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
560  // double dist_to_tile = _distance_to_tile(jetA, *RTile);
561  // bool relevant_for_jetA = dist_to_tile <= jetA->NN_dist;
562  // bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
563  // if (relevant_for_jetA || relevant_for_RTile) {
564  // for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
565  // double dist = _bj_dist_not_periodic(jetA,jetB);
566  // if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
567  // if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
568  // }
569  // }
570  // }
571  // }
572  // }
573  // no need to do it for LH tiles, since they are implicitly done
574  // when we set NN for both jetA and jetB on the RH tiles.
575  }
576  // Now update the max_NN_dist within each tile. Not strictly
577  // necessary, because existing max_NN_dist is an upper bound. but
578  // costs little and may give some efficiency gain later.
579  // (Do it only for real particles -- ghosts don't come into the game).
580  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
581  tile->max_NN_dist = 0;
582  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
583  if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
584  }
585  }
586 
587 
588  vector<double> diJs(nreal);
589  for (int i = 0; i < nreal; i++) {
590  diJs[i] = _bj_diJ(&briefjets[i]);
591  briefjets[i].label_minheap_update_done();
592  }
593  MinHeap minheap(diJs);
594  // have a stack telling us which jets we'll have to update on the heap
595  vector<TiledJet3 *> jets_for_minheap;
596  jets_for_minheap.reserve(ntot); // GPS TMP: return here
597 
598  // now run the recombination loop
599  while (nreal > 0) {
600 
601  double diJ_min = minheap.minval() *_invR2;
602  jetA = head + minheap.minloc();
603 
604  // do the recombination between A and B
605  jetB = jetA->NN;
606 
607  if (jetB != NULL) {
608  // jet-jet recombination
609  // If necessary relabel A & B to ensure jetB < jetA, that way if
610  // the larger of them == newtail then that ends up being jetA and
611  // the new jet that is added as jetB is inserted in a position that
612  // has a future!
613  if (jetA < jetB) {std::swap(jetA,jetB);}
614 
615  int nn; // new jet index
616  _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
617 
618  // what was jetB will now become the new jet
619  _bj_remove_from_tiles(jetA);
620  oldB = * jetB; // take a copy because we will need it...
621  _bj_remove_from_tiles(jetB);
622  bool is_ghost = false;
623  _tj_set_jetinfo(jetB, nn, is_ghost); // cause jetB to become _jets[nn]
624  // (also registers the jet in the tiling); it will never be a ghost
625  } else {
626  // jet-beam recombination
627  // get the hist_index
628  _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
629  _bj_remove_from_tiles(jetA);
630  }
631 
632  // remove the minheap entry for jetA if jetA is a real particle;
633  // only in this case
634  if (!jetA->is_ghost) {
635  minheap.remove(jetA-head);
636  // jetB cannot be a ghost; so when jetA is not a ghost we decrease
637  // our count of remaining real particles
638  nreal--;
639  }
640 
641 
642 
643 
644  // first establish the set of tiles over which we are going to
645  // have to run searches for updated and new nearest-neighbours --
646  // basically a combination of vicinity of the tiles of the two old
647  // and one new jet.
648  int n_near_tiles = 0;
649  // add tiles neighbouring A even if A is a ghost because a ghost
650  // could have been the NN of some other particle
651  _add_untagged_neighbours_to_tile_union_using_max_info(jetA,
652  tile_union, n_near_tiles);
653  if (jetB != NULL) {
654  _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
655  tile_union,n_near_tiles);
656  jetB->label_minheap_update_needed();
657  jets_for_minheap.push_back(jetB);
658  }
659 
660 
661  // Initialise jetB's NN distance as well as updating it for
662  // other particles.
663  // Run over all tiles in our union
664  if (jetB != NULL) {
665  Tile3 & jetB_tile = _tiles[jetB->tile_index];
666  for (Tile3 ** near_tile = jetB_tile.begin_tiles;
667  near_tile != jetB_tile.end_tiles; near_tile++) {
668 
669  double dist_to_tile = _distance_to_tile(jetB, *near_tile);
670  // use <= in next line so that on first tile, relevant_for_jetB is
671  // set to true
672  bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
673  bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
674  bool relevant = relevant_for_jetB || relevant_for_near_tile;
675  // this first option decides exactly what loop to do based on whether
676  // the near tile was tagged. You'd think it's more efficient, but
677  // not necessarily...
678  if (relevant) {
679  if ((*near_tile)->tagged) {
680  for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
681  if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
682  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
683  }
684  (*near_tile)->tagged = false;
685  } else {
686  for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
687  _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
688  }
689  }
690  }
691  // Now HANDLE GHOSTS
692  if (relevant_for_jetB) {
693  for (TiledJet3 * jetI = (*near_tile)->ghost_head; jetI != NULL; jetI = jetI->next) {
694  double dist = _bj_dist(jetB,jetI);
695  if (dist < jetB->NN_dist) {
696  jetB->NN_dist = dist; jetB->NN = jetI;
697  }
698  }
699  }
700 
701 
702  // this second option does everything independently of whether the near tile
703  // was tagged -- somehow you'd expect it to be slower, but it may actually be
704  // marginally faster.
705  // if (relevant_for_jetB || relevant_for_near_tile) {
706  // for (TiledJet3 * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
707  //
708  // if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
709  // _set_NN(jetI, jets_for_minheap);
710  // }
711  //
712  // _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
713  // // -- Keep this old inline code for later speed tests
714  // // double dist = _bj_dist(jetI,jetB);
715  // // if (dist < jetI->NN_dist) {
716  // // if (jetI != jetB) {
717  // // jetI->NN_dist = dist;
718  // // jetI->NN = jetB;
719  // // // label jetI as needing heap action...
720  // // if (!jetI->minheap_update_needed()) {
721  // // jetI->label_minheap_update_needed();
722  // // jets_for_minheap.push_back(jetI);
723  // // }
724  // // }
725  // // }
726  // // if (dist < jetB->NN_dist) {
727  // // if (jetI != jetB) {
728  // // jetB->NN_dist = dist;
729  // // jetB->NN = jetI;}
730  // // }
731  // }
732  // (*near_tile)->tagged = false;
733  // }
734  }
735  }
736 
737  // now run over the tiles that were tagged earlier and that we haven't yet
738  // had a change to visit.
739  for (int itile = 0; itile < n_near_tiles; itile++) {
740  Tile3 * tile_ptr = &_tiles[tile_union[itile]];
741  if (!tile_ptr->tagged) continue; // because earlier loop may have undone the tag
742  tile_ptr->tagged = false;
743  // run over all jets in the current tile
744  for (TiledJet3 * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
745  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
746  if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
747  _set_NN(jetI, jets_for_minheap);
748  }
749  }
750  }
751 
752  // deal with jets whose minheap entry needs updating
753  //if (verbose) cout << " jets whose NN was modified: " << endl;
754  while (jets_for_minheap.size() > 0) {
755  TiledJet3 * jetI = jets_for_minheap.back();
756  jets_for_minheap.pop_back();
757  minheap.update(jetI-head, _bj_diJ(jetI));
758  jetI->label_minheap_update_done();
759  // handle max_NN_dist update for all jets that might have
760  // seen a change (increase) of distance
761  Tile3 & tile_I = _tiles[jetI->tile_index];
762  if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
763  }
764  }
765 
766  // final cleaning up;
767  delete[] briefjets;
768 }
769 
770 
771 FASTJET_END_NAMESPACE