Fork me on GitHub

source: git/external/fastjet/LazyTiling9Alt.cc@ 81bedba

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 81bedba 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: 28.4 KB
RevLine 
[35cdc46]1//FJSTARTHEADER
[1d208a2]2// $Id: LazyTiling9Alt.cc 3807 2015-02-20 11:16:55Z 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 <iomanip>
32#include "fastjet/internal/LazyTiling9Alt.hh"
33#include "fastjet/internal/TilingExtent.hh"
34using namespace std;
35
36FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37
38LazyTiling9Alt::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///
68void 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
194int 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
216inline 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//----------------------------------------------------------------------
236void 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
256void 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)
279void 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.
294inline 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.
314inline 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;
[10e33bc]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;
[35cdc46]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--------------------
338ostream & 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...
391inline 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
412inline 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
458void 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
726FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.