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