Fork me on GitHub

source: git/external/fastjet/LazyTiling9.cc@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 28.3 KB
Line 
1//FJSTARTHEADER
2// $Id: LazyTiling9.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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"
37using namespace std;
38
39// uncomment the line below to use TilingExtent in LazyTiling9
40#define _FASTJET_TILING2_USE_TILING_ANALYSIS_
41
42FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
43
44
45
46LazyTiling9::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///
80void 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
223int 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
245inline 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//----------------------------------------------------------------------
265void 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
285void 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)
308void 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.
323inline 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.
343inline 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
379inline 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...
423inline 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
444inline 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
490void 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
751FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.