Fork me on GitHub

source: git/external/fastjet/LazyTiling25.cc@ bdf9e2e

Last change on this file since bdf9e2e was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4

  • Property mode set to 100644
File size: 30.4 KB
Line 
1//FJSTARTHEADER
2// $Id: LazyTiling25.cc 3807 2015-02-20 11:16:55Z 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"
34using namespace std;
35
36
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39LazyTiling25::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///
75void 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
238int 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
260inline 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//----------------------------------------------------------------------
280void 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
300void 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)
325void 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.
340inline 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.
360inline 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 // << " -> diff = " << dist-(*near_tile)->max_NN_dist << 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
396inline 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...
440inline 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
461inline 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
507void 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 // GS comment:
668 //
669 // At this stage, we have perforned the clustering in
670 // ClusterSequence and we need to update the NNs. The objects we
671 // deal with are jetA and oldB (the once that have been clustered)
672 // as well as jetB (the result of the clustering)
673 //
674 // There are two types of objects we need to update:
675 // - find jetB NN
676 // - update the NN of points which had jetA or jetB as their NN
677 //
678 // Wile we find jetB's NN, browsing relevant tiles near jetB, we
679 // also search for points which had jetA or jetB as their
680 // NN. These are tagged. Then we list the relevant tiles where we
681 // can potentially have points to update (i.e. points which had
682 // jetA and oldB as NN) in the yet untagged relevant tiles near
683 // jetA and oldB.
684 //
685 // DEBUG:
686 // if (jetB != NULL) {
687 // cout << "jetA = " << jetA->_jets_index << " (tile " << jetA->tile_index << "), "
688 // << "oldB = " << oldB._jets_index << " (tile " << oldB. tile_index << "), "
689 // << "jetB = " << jetB->_jets_index << " (tile " << jetB->tile_index << ")" << endl;
690 // } else {
691 // cout << "jetA = " << jetA->_jets_index << " (tile " << jetA->tile_index << ")" << endl;
692 // }
693
694 // Initialise jetB's NN distance as well as updating it for other
695 // particles. While doing so, examine whether jetA or old jetB was
696 // some other particle's NN.
697 if (jetB != NULL) {
698 Tile25 & jetB_tile = _tiles[jetB->tile_index];
699 for (Tile25 ** near_tile = jetB_tile.begin_tiles;
700 near_tile != jetB_tile.end_tiles; near_tile++) {
701
702 double dist_to_tile = _distance_to_tile(jetB, *near_tile);
703 // use <= in next line so that on first tile, relevant_for_jetB is
704 // set to true
705 bool relevant_for_jetB = dist_to_tile <= jetB->NN_dist;
706 bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
707 bool relevant = relevant_for_jetB || relevant_for_near_tile;
708
709 // cout << " Relevance of tile " << *near_tile - & _tiles[0]
710 // << " wrt jetB is " << relevant << endl;
711
712 if (! relevant) continue;
713 // now label this tile as having been considered (so that we
714 // don't go over it again later)
715 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
716 (*near_tile)->tagged = true;
717 n_near_tiles++;
718
719 // if going over the neighbouring tile's jets, check anyway
720 // whether A or B were nearest neighbours, since it comes at a
721 // modest cost relative to the distance computation (and we would
722 // in most cases have to do it again later anyway).
723 for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
724 if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
725 _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
726 }
727 }
728 }
729
730 // first establish the set of tiles over which we are going to
731 // have to run searches for updated and new nearest-neighbours --
732 // basically a combination of vicinity of the tiles of the two old
733 // and one new jet.
734 int n_done_tiles = n_near_tiles;
735 //cout << "Looking at relevant tiles to update for jetA" << endl;
736 _add_untagged_neighbours_to_tile_union_using_max_info(jetA,
737 tile_union, n_near_tiles);
738 if (jetB != NULL) {
739 // cout << "Looking at relevant tiles to update for oldB" << endl;
740 _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
741 tile_union,n_near_tiles);
742 jetB->label_minheap_update_needed();
743 jets_for_minheap.push_back(jetB);
744 }
745
746
747 // first untag the tiles we have already dealt with
748 for (int itile = 0; itile < n_done_tiles; itile++) {
749 _tiles[tile_union[itile]].tagged = false;
750 }
751 // now run over the tiles that were tagged earlier and that we haven't yet
752 // had a change to visit.
753 for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
754 Tile25 * tile_ptr = &_tiles[tile_union[itile]];
755 tile_ptr->tagged = false;
756 // run over all jets in the current tile
757 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
758 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
759 if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
760 _set_NN(jetI, jets_for_minheap);
761 }
762 }
763 }
764
765 // deal with jets whose minheap entry needs updating
766 //if (verbose) cout << " jets whose NN was modified: " << endl;
767 while (jets_for_minheap.size() > 0) {
768 TiledJet * jetI = jets_for_minheap.back();
769 jets_for_minheap.pop_back();
770 minheap.update(jetI-head, _bj_diJ(jetI));
771 jetI->label_minheap_update_done();
772 // handle max_NN_dist update for all jets that might have
773 // seen a change (increase) of distance
774 Tile25 & tile_I = _tiles[jetI->tile_index];
775 if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
776 }
777 n--;
778 }
779
780 // final cleaning up;
781 delete[] briefjets;
782#ifdef INSTRUMENT2
783 cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
784#endif // INSTRUMENT2
785
786}
787
788FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.