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