Fork me on GitHub

source: git/external/fastjet/LazyTiling9SeparateGhosts.cc@ d77b51d

ImprovedOutputFile Timing dual_readout llp
Last change on this file since d77b51d was 10e33bc, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

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