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