Fork me on GitHub

source: git/external/fastjet/NNFJN2Tiled.hh@ 5eaa70a

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 5eaa70a 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: 26.8 KB
Line 
1#ifndef __FASTJET_NNFJN2TILED_HH__
2#define __FASTJET_NNFJN2TILED_HH__
3
4//FJSTARTHEADER
5// $Id: NNFJN2Tiled.hh 4056 2016-03-03 15:27:35Z soyez $
6//
7// Copyright (c) 2016, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8//
9//----------------------------------------------------------------------
10// This file is part of FastJet.
11//
12// FastJet is free software; you can redistribute it and/or modify
13// it under the terms of the GNU General Public License as published by
14// the Free Software Foundation; either version 2 of the License, or
15// (at your option) any later version.
16//
17// The algorithms that underlie FastJet have required considerable
18// development. They are described in the original FastJet paper,
19// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
20// FastJet as part of work towards a scientific publication, please
21// quote the version you use and include a citation to the manual and
22// optionally also to hep-ph/0512210.
23//
24// FastJet is distributed in the hope that it will be useful,
25// but WITHOUT ANY WARRANTY; without even the implied warranty of
26// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27// GNU General Public License for more details.
28//
29// You should have received a copy of the GNU General Public License
30// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
31//----------------------------------------------------------------------
32//FJENDHEADER
33
34#include <fastjet/NNBase.hh>
35#include <fastjet/internal/TilingExtent.hh>
36
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39//----------------------------------------------------------------------
40/// @ingroup advanced_usage
41/// \class NNFJN2Tiled
42///
43/// Helps solve closest pair problems with factorised interparticle
44/// and beam distances (ie satisfying the FastJet lemma) that are on
45/// a cylindrical geometry and allow tiling.
46///
47/// (see NNBase.hh for an introductory description)
48///
49/// This variant provides an implementation based on the N2Tiled
50/// clustering strategy in FastJet. As for the NNFJN2Plain case, the
51/// interparticle and beam distances should be of the form
52///
53/// \code
54/// dij = min(mom_factor(i), mom_factor(j)) * geometrical_distance(i,j)
55/// diB = mom_factor(i) * geometrical_beam_distance(i)
56/// \endcode
57///
58/// Additionally, the NNFJN2Tiled class takes a tile_size parameter
59/// that controls the size of the tiles. It must be such that, for any
60/// two points in non-neighbouring (and non-identical) tiles, the
61/// geometrical distance between the 2 points is larger than the
62/// geometrical beam distance of each of the 2 points.
63///
64/// It is templated with a BJ (brief jet) class and can be used with or
65/// without an extra "Information" template, i.e. NNFJN2Tiled<BJ> or
66/// NNFJN2Tiled<BJ,I>
67///
68/// For the NNFJN2Tiled<BJ> version of the class to function, BJ must provide
69/// three member functions
70///
71/// \code
72/// void BJ::init(const PseudoJet & jet); // initialise with a PseudoJet
73/// double BJ::geometrical_distance(const BJ * other_bj_jet); // distance between this and other_bj_jet (geometrical part)
74/// double BJ::geometrical_beam_distance(); // distance to the beam (geometrical part)
75/// double BJ::momentum_factor(); // extra momentum factor
76/// \endcode
77///
78/// For the NNFJN2Tiled<BJ,I> version to function, the BJ::init(...) member
79/// must accept an extra argument
80///
81/// \code
82/// void BJ::init(const PseudoJet & jet, I * info); // initialise with a PseudoJet + info
83/// \endcode
84///
85/// NOTE: THE DISTANCE MUST BE SYMMETRIC I.E. SATISFY
86/// \code
87/// a.geometrical_distance(b) == b.geometrical_distance(a)
88/// \endcode
89///
90/// Finally, the BJ class needs to provide access to the variables used
91/// for the rectangular tiling:
92///
93/// \code
94/// double BJ::rap(); // rapidity-like variable
95/// double BJ::phi(); // azimutal-angle-like variable (should be > -2pi)
96/// \endcode
97///
98/// Note that you are strongly advised to add the following lines
99/// to your BJ class to allow it to be used also with NNH:
100///
101/// \code
102/// /// make this BJ class compatible with the use of NNH
103/// double BJ::distance(const BJ * other_bj_jet){
104/// double mom1 = momentum_factor();
105/// double mom2 = other_bj_jet->momentum_factor();
106/// return (mom1<mom2 ? mom1 : mom2) * geometrical_distance(other_bj_jet);
107/// }
108/// double BJ::beam_distance(){
109/// return momentum_factor() * geometrical_beam_distance();
110/// }
111/// \endcode
112///
113template<class BJ, class I = _NoInfo> class NNFJN2Tiled : public NNBase<I> {
114public:
115
116 /// constructor with an initial set of jets (which will be assigned indices
117 /// `0...jets.size()-1`)
118 NNFJN2Tiled(const std::vector<PseudoJet> & jets, double requested_tile_size)
119 : NNBase<I>(), _requested_tile_size(requested_tile_size) {start(jets);}
120 NNFJN2Tiled(const std::vector<PseudoJet> & jets, double requested_tile_size, I * info)
121 : NNBase<I>(info), _requested_tile_size(requested_tile_size) {start(jets);}
122
123 void start(const std::vector<PseudoJet> & jets);
124
125 /// return the dij_min and indices iA, iB, for the corresponding jets.
126 /// If iB < 0 then iA recombines with the beam
127 double dij_min(int & iA, int & iB);
128
129 /// remove the jet pointed to by index iA
130 void remove_jet(int iA);
131
132 /// merge the jets pointed to by indices A and B and replace them with
133 /// jet, assigning it an index jet_index.
134 void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
135
136 /// a destructor
137 ~NNFJN2Tiled() {
138 delete[] briefjets;
139 delete[] diJ;
140 }
141
142private:
143 class TiledJet; // forward declaration
144 class Tile;
145 class diJ_plus_link;
146
147 // Set up the tiles:
148 void _initialise_tiles(const std::vector<PseudoJet> & particles);
149
150 // return the full distance of a particle to its NN
151 inline double _compute_diJ(const TiledJet * const jet) const {
152 double mom_fact = jet->momentum_factor();
153 if (jet->NN != NULL) {
154 double other_mom_fact = jet->NN->momentum_factor();
155 if (other_mom_fact < mom_fact) {mom_fact = other_mom_fact;}
156 }
157 return jet->NN_dist * mom_fact;
158 }
159
160 // reasonably robust return of tile index given irap and iphi, in particular
161 // it works even if iphi is negative
162 inline int _tile_index (int irap, int iphi) const {
163 // note that (-1)%n = -1 so that we have to add _n_tiles_phi
164 // before performing modulo operation
165 return (irap-_tiles_irap_min)*_n_tiles_phi
166 + (iphi+_n_tiles_phi) % _n_tiles_phi;
167 }
168
169 int _tile_index(const double rap, const double phi) const;
170 void _tiledjet_set_jetinfo ( TiledJet * const tiled_jet, const PseudoJet &jet, int index);
171 void _bj_remove_from_tiles(TiledJet * const jet);
172 void _initialise_tiles();
173 void _print_tiles(TiledJet * briefjets ) const;
174 void _add_neighbours_to_tile_union(const int tile_index, int & n_near_tiles) const;
175 void _add_untagged_neighbours_to_tile_union(const int tile_index, int & n_near_tiles);
176
177
178 /// contains the briefjets
179 TiledJet * briefjets;
180
181 /// semaphores for the current extent of our structure
182 TiledJet * head;
183
184 /// currently active number of jets
185 int n;
186
187 /// where_is[i] contains a pointer to the jet with index i
188 std::vector<TiledJet *> where_is;
189
190 /// helper to keep tracks of tiles to be checked for updates
191 std::vector<int> tile_union;
192
193 /// a table containing all the (full) distances to each object's NN
194 diJ_plus_link * diJ;
195
196 /// tiling information
197 std::vector<Tile> _tiles;
198 double _requested_tile_size;
199 double _tiles_rap_min, _tiles_rap_max;
200 double _tile_size_rap, _tile_size_phi;
201 int _n_tiles_phi,_tiles_irap_min,_tiles_irap_max;
202
203 /// a class that wraps around the BJ, supplementing it with extra information
204 /// such as pointers to neighbours, etc.
205 class TiledJet : public BJ {
206 public:
207 void init(const PseudoJet & jet, int index_in) {
208 BJ::init(jet);
209 other_init(index_in);
210 }
211 void init(const PseudoJet & jet, int index_in, I * info) {
212 BJ::init(jet, info);
213 other_init(index_in);
214 }
215 void other_init(int index_in) {
216 _index = index_in;
217 NN_dist = BJ::geometrical_beam_distance();
218 NN = NULL;
219 }
220 int jet_index() const {return _index;}
221
222 double NN_dist;
223 TiledJet * NN, *previous, * next;
224 int tile_index, diJ_posn;
225 // routines that are useful in the minheap version of tiled
226 // clustering ("misuse" the otherwise unused diJ_posn, so as
227 // to indicate whether jets need to have their minheap entries
228 // updated).
229 inline void label_minheap_update_needed() {diJ_posn = 1;}
230 inline void label_minheap_update_done() {diJ_posn = 0;}
231 inline bool minheap_update_needed() const {return diJ_posn==1;}
232
233 private:
234 int _index;
235 };
236
237 /// number of neighbours that a tile will have (rectangular geometry
238 /// gives 9 neighbours).
239 static const int n_tile_neighbours = 9;
240 //----------------------------------------------------------------------
241 /// The fundamental structures to be used for the tiled N^2 algorithm
242 /// (see CCN27-44 for some discussion of pattern of tiling)
243 class Tile {
244 public:
245 /// pointers to neighbouring tiles, including self
246 Tile * begin_tiles[n_tile_neighbours];
247 /// neighbouring tiles, excluding self
248 Tile ** surrounding_tiles;
249 /// half of neighbouring tiles, no self
250 Tile ** RH_tiles;
251 /// just beyond end of tiles
252 Tile ** end_tiles;
253 /// start of list of BriefJets contained in this tile
254 TiledJet * head;
255 /// sometimes useful to be able to tag a tile
256 bool tagged;
257 };
258
259 // structure that holds the real, full, distance (as well as a pointer to the corresponding TiledJet)
260 struct diJ_plus_link {
261 double diJ; // the distance
262 TiledJet * jet; // the jet (i) for which we've found this distance
263 // (whose NN will the J).
264 };
265
266};
267
268
269
270//----------------------------------------------------------------------
271template<class BJ, class I> void NNFJN2Tiled<BJ,I>::start(const std::vector<PseudoJet> & jets) {
272
273 _initialise_tiles(jets);
274
275 n = jets.size();
276
277 briefjets = new TiledJet[n];
278 where_is.resize(2*n);
279
280 TiledJet * jetA = briefjets, * jetB;
281
282 // will be used quite deep inside loops, but declare it here so that
283 // memory (de)allocation gets done only once
284 tile_union.resize(3*n_tile_neighbours);
285
286 // initialise the basic jet info
287 for (int i = 0; i< n; i++) {
288 _tiledjet_set_jetinfo(jetA, jets[i], i);
289 where_is[i] = jetA;
290 jetA++; // move on to next entry of briefjets
291 }
292
293 head = briefjets; // a nicer way of naming start
294
295 // set up the initial nearest neighbour information
296 typename std::vector<Tile>::const_iterator tile;
297 for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
298 // first do it on this tile
299 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
300 for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
301 double dist = jetA->geometrical_distance(jetB);
302 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
303 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
304 }
305 }
306 // then do it for RH tiles
307 for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
308 for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
309 for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
310 double dist = jetA->geometrical_distance(jetB);
311 if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
312 if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
313 }
314 }
315 }
316 // no need to do it for LH tiles, since they are implicitly done
317 // when we set NN for both jetA and jetB on the RH tiles.
318 }
319
320 diJ = new diJ_plus_link[n];
321 jetA = head;
322 for (int i = 0; i < n; i++) {
323 diJ[i].diJ = _compute_diJ(jetA); // kt distance * R^2
324 diJ[i].jet = jetA; // our compact diJ table will not be in
325 jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
326 // so set up bi-directional correspondence here.
327 jetA++; // have jetA follow i
328 }
329}
330
331
332//----------------------------------------------------------------------
333template<class BJ, class I> double NNFJN2Tiled<BJ,I>::dij_min(int & iA, int & iB) {
334 // find the minimum of the diJ on this round
335 diJ_plus_link * best, *stop; // pointers a bit faster than indices
336 // could use best to keep track of diJ
337 // min, but it turns out to be
338 // marginally faster to have a separate
339 // variable (avoids n dereferences at
340 // the expense of n/2 assignments).
341 double diJ_min = diJ[0].diJ; // initialise the best one here.
342 best = diJ; // and here
343 stop = diJ+n;
344 for (diJ_plus_link * here = diJ+1; here != stop; here++) {
345 if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
346 }
347
348 // return information to user about recombination
349 TiledJet * jetA = best->jet;
350 iA = jetA->jet_index();
351 iB = jetA->NN ? jetA->NN->jet_index() : -1;
352 return diJ_min;
353}
354
355
356//----------------------------------------------------------------------
357// remove jetA from the list
358template<class BJ, class I> void NNFJN2Tiled<BJ,I>::remove_jet(int iA) {
359 TiledJet * jetA = where_is[iA];
360
361 _bj_remove_from_tiles(jetA);
362
363 // first establish the set of tiles over which we are going to
364 // have to run searches for updated and new nearest-neighbours --
365 // basically a combination of vicinity of the tiles of the two old
366 // and one new jet.
367 int n_near_tiles = 0;
368 _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
369
370 // now update our nearest neighbour info and diJ table
371 // first reduce size of table
372 n--;
373 // then compactify the diJ by taking the last of the diJ and copying
374 // it to the position occupied by the diJ for jetA
375 diJ[n].jet->diJ_posn = jetA->diJ_posn;
376 diJ[jetA->diJ_posn] = diJ[n];
377
378 // updating other particles' NN.
379 // Run over all tiles in our union
380 for (int itile = 0; itile < n_near_tiles; itile++) {
381 Tile * tile_ptr = &_tiles[tile_union[itile]];
382 tile_ptr->tagged = false; // reset tag, since we're done with unions
383 // run over all jets in the current tile
384 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
385 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
386 if (jetI->NN == jetA) {
387 jetI->NN_dist = jetI->geometrical_beam_distance();
388 jetI->NN = NULL;
389 // now go over tiles that are neighbours of I (include own tile)
390 for (Tile ** near_tile = tile_ptr->begin_tiles;
391 near_tile != tile_ptr->end_tiles; near_tile++) {
392 // and then over the contents of that tile
393 for (TiledJet * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
394 double dist = jetI->geometrical_distance(jetJ);
395 if (dist < jetI->NN_dist && jetJ != jetI) {
396 jetI->NN_dist = dist; jetI->NN = jetJ;
397 }
398 }
399 }
400 diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
401 }
402 }
403 }
404
405}
406
407
408//----------------------------------------------------------------------
409template<class BJ, class I> void NNFJN2Tiled<BJ,I>::merge_jets(int iA, int iB,
410 const PseudoJet & jet, int index) {
411
412 TiledJet * jetA = where_is[iA];
413 TiledJet * jetB = where_is[iB];
414
415 // jet-jet recombination
416 // If necessary relabel A & B to ensure jetB < jetA, that way if
417 // the larger of them == newtail then that ends up being jetA and
418 // the new jet that is added as jetB is inserted in a position that
419 // has a future!
420 if (jetA < jetB) {std::swap(jetA,jetB);}
421
422 // what was jetB will now become the new jet
423 _bj_remove_from_tiles(jetA);
424 TiledJet oldB = * jetB; // take a copy because we will need it...
425 _bj_remove_from_tiles(jetB);
426 _tiledjet_set_jetinfo(jetB, jet, index); // cause jetB to become _jets[nn]
427 // (also registers the jet in the tiling)
428 where_is[index] = jetB;
429
430 // first establish the set of tiles over which we are going to
431 // have to run searches for updated and new nearest-neighbours --
432 // basically a combination of vicinity of the tiles of the two old
433 // and one new jet.
434 int n_near_tiles = 0;
435 _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
436 if (jetB->tile_index != jetA->tile_index) {
437 _add_untagged_neighbours_to_tile_union(jetB->tile_index, n_near_tiles);
438 }
439 if (oldB.tile_index != jetA->tile_index &&
440 oldB.tile_index != jetB->tile_index) {
441 _add_untagged_neighbours_to_tile_union(oldB.tile_index, n_near_tiles);
442 }
443
444 // now update our nearest neighbour info and diJ table
445 // first reduce size of table
446 n--;
447 // then compactify the diJ by taking the last of the diJ and copying
448 // it to the position occupied by the diJ for jetA
449 diJ[n].jet->diJ_posn = jetA->diJ_posn;
450 diJ[jetA->diJ_posn] = diJ[n];
451
452 // Initialise jetB's NN distance as well as updating it for
453 // other particles.
454 // Run over all tiles in our union
455 for (int itile = 0; itile < n_near_tiles; itile++) {
456 Tile * tile_ptr = &_tiles[tile_union[itile]];
457 tile_ptr->tagged = false; // reset tag, since we're done with unions
458 // run over all jets in the current tile
459 for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
460 // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
461 if ((jetI->NN == jetA) || (jetI->NN == jetB)) {
462 jetI->NN_dist = jetI->geometrical_beam_distance();
463 jetI->NN = NULL;
464 // now go over tiles that are neighbours of I (include own tile)
465 for (Tile ** near_tile = tile_ptr->begin_tiles; near_tile != tile_ptr->end_tiles; near_tile++) {
466 // and then over the contents of that tile
467 for (TiledJet * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
468 double dist = jetI->geometrical_distance(jetJ);
469 if (dist < jetI->NN_dist && jetJ != jetI) {
470 jetI->NN_dist = dist; jetI->NN = jetJ;
471 }
472 }
473 }
474 diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
475 }
476 // check whether new jetB is closer than jetI's current NN and
477 // if jetI is closer than jetB's current (evolving) nearest
478 // neighbour. Where relevant update things
479 double dist = jetI->geometrical_distance(jetB);
480 if (dist < jetI->NN_dist) {
481 if (jetI != jetB) {
482 jetI->NN_dist = dist;
483 jetI->NN = jetB;
484 diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ...
485 }
486 }
487 if (dist < jetB->NN_dist) {
488 if (jetI != jetB) {
489 jetB->NN_dist = dist;
490 jetB->NN = jetI;}
491 }
492 }
493 }
494
495 // finally, register the updated kt distance for B
496 diJ[jetB->diJ_posn].diJ = _compute_diJ(jetB);
497}
498
499
500//----------------------------------------------------------------------
501/// Set up the tiles:
502/// - decide the range in eta
503/// - allocate the tiles
504/// - set up the cross-referencing info between tiles
505///
506/// The neighbourhood of a tile is set up as follows
507///
508/// LRR
509/// LXR
510/// LLR
511///
512/// such that tiles is an array containing XLLLLRRRR with pointers
513/// | \ RH_tiles
514/// \ surrounding_tiles
515///
516/// with appropriate precautions when close to the edge of the tiled
517/// region.
518///
519template <class BJ, class I>
520void NNFJN2Tiled<BJ,I>::_initialise_tiles(const std::vector<PseudoJet> &particles) {
521
522 // first decide tile sizes (with a lower bound to avoid huge memory use with
523 // very small R)
524 double default_size = _requested_tile_size>0.1 ? _requested_tile_size : 0.1;
525 _tile_size_rap = default_size;
526 // it makes no sense to go below 3 tiles in phi -- 3 tiles is
527 // sufficient to make sure all pair-wise combinations up to pi in
528 // phi are possible
529 _n_tiles_phi = int(floor(twopi/default_size));
530 if (_n_tiles_phi<3) _n_tiles_phi = 3;
531 _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
532
533 TilingExtent tiling_analysis(particles);
534 _tiles_rap_min = tiling_analysis.minrap();
535 _tiles_rap_max = tiling_analysis.maxrap();
536
537 // now adjust the values
538 _tiles_irap_min = int(floor(_tiles_rap_min/_tile_size_rap));
539 _tiles_irap_max = int(floor( _tiles_rap_max/_tile_size_rap));
540 _tiles_rap_min = _tiles_irap_min * _tile_size_rap;
541 _tiles_rap_max = _tiles_irap_max * _tile_size_rap;
542
543 // allocate the tiles
544 _tiles.resize((_tiles_irap_max-_tiles_irap_min+1)*_n_tiles_phi);
545
546 // now set up the cross-referencing between tiles
547 for (int irap = _tiles_irap_min; irap <= _tiles_irap_max; irap++) {
548 for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
549 Tile * tile = & _tiles[_tile_index(irap,iphi)];
550 // no jets in this tile yet
551 tile->head = NULL; // first element of tiles points to itself
552 tile->begin_tiles[0] = tile;
553 Tile ** pptile = & (tile->begin_tiles[0]);
554 pptile++;
555 //
556 // set up L's in column to the left of X
557 tile->surrounding_tiles = pptile;
558 if (irap > _tiles_irap_min) {
559 // with the itile subroutine, we can safely run tiles from
560 // idphi=-1 to idphi=+1, because it takes care of
561 // negative and positive boundaries
562 for (int idphi = -1; idphi <=+1; idphi++) {
563 *pptile = & _tiles[_tile_index(irap-1,iphi+idphi)];
564 pptile++;
565 }
566 }
567 // now set up last L (below X)
568 *pptile = & _tiles[_tile_index(irap,iphi-1)];
569 pptile++;
570 // set up first R (above X)
571 tile->RH_tiles = pptile;
572 *pptile = & _tiles[_tile_index(irap,iphi+1)];
573 pptile++;
574 // set up remaining R's, to the right of X
575 if (irap < _tiles_irap_max) {
576 for (int idphi = -1; idphi <= +1; idphi++) {
577 *pptile = & _tiles[_tile_index(irap+1,iphi+idphi)];
578 pptile++;
579 }
580 }
581 // now put semaphore for end tile
582 tile->end_tiles = pptile;
583 // finally make sure tiles are untagged
584 tile->tagged = false;
585 }
586 }
587
588}
589
590//----------------------------------------------------------------------
591/// return the tile index corresponding to the given rap,phi point
592template <class BJ, class I>
593int NNFJN2Tiled<BJ,I>::_tile_index(const double rap, const double phi) const {
594 int irap, iphi;
595 if (rap <= _tiles_rap_min) {irap = 0;}
596 else if (rap >= _tiles_rap_max) {irap = _tiles_irap_max-_tiles_irap_min;}
597 else {
598 //irap = int(floor((rap - _tiles_rap_min) / _tile_size_rap));
599 irap = int(((rap - _tiles_rap_min) / _tile_size_rap));
600 // following needed in case of rare but nasty rounding errors
601 if (irap > _tiles_irap_max-_tiles_irap_min) {
602 irap = _tiles_irap_max-_tiles_irap_min;}
603 }
604 // allow for some extent of being beyond range in calculation of phi
605 // as well
606 //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
607 // with just int and no floor, things run faster but beware
608 iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
609 return (iphi + irap * _n_tiles_phi);
610}
611
612//----------------------------------------------------------------------
613template <class BJ, class I>
614void NNFJN2Tiled<BJ,I>::_bj_remove_from_tiles(TiledJet * const jet) {
615 Tile * tile = & _tiles[jet->tile_index];
616
617 if (jet->previous == NULL) {
618 // we are at head of the tile, so reset it.
619 // If this was the only jet on the tile then tile->head will now be NULL
620 tile->head = jet->next;
621 } else {
622 // adjust link from previous jet in this tile
623 jet->previous->next = jet->next;
624 }
625 if (jet->next != NULL) {
626 // adjust backwards-link from next jet in this tile
627 jet->next->previous = jet->previous;
628 }
629}
630
631
632//----------------------------------------------------------------------
633// overloaded version which additionally sets up information regarding the
634// tiling
635template <class BJ, class I>
636inline void NNFJN2Tiled<BJ,I>::_tiledjet_set_jetinfo(TiledJet * const tile_jet,
637 const PseudoJet &jet,
638 int index) {
639 // the this-> in the next line is required by standard compiler
640 // see e.g. http://stackoverflow.com/questions/10639053/name-lookups-in-c-templates
641 this->init_jet(tile_jet, jet, index);
642
643 // Then do the setup specific to the tiled case.
644
645 // Find out which tile it belonds to
646 tile_jet->tile_index = _tile_index(tile_jet->rap(), tile_jet->phi());
647
648 // Insert it into the tile's linked list of jets
649 Tile * tile = &_tiles[tile_jet->tile_index];
650 tile_jet->previous = NULL;
651 tile_jet->next = tile->head;
652 if (tile_jet->next != NULL) {tile_jet->next->previous = tile_jet;}
653 tile->head = tile_jet;
654}
655
656//----------------------------------------------------------------------
657/// Add to the vector tile_union the tiles that are in the neighbourhood
658/// of the specified tile_index, including itself -- start adding
659/// from position n_near_tiles-1, and increase n_near_tiles as
660/// you go along (could have done it more C++ like with vector with reserved
661/// space, but fear is that it would have been slower, e.g. checking
662/// for end of vector at each stage to decide whether to resize it)
663template <class BJ, class I>
664void NNFJN2Tiled<BJ,I>::_add_neighbours_to_tile_union(const int tile_index,
665 int & n_near_tiles) const {
666 for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
667 near_tile != _tiles[tile_index].end_tiles; near_tile++){
668 // get the tile number
669 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
670 n_near_tiles++;
671 }
672}
673
674//----------------------------------------------------------------------
675/// Like _add_neighbours_to_tile_union, but only adds neighbours if
676/// their "tagged" status is false; when a neighbour is added its
677/// tagged status is set to true.
678///
679/// Note that with a high level of warnings (-pedantic -Wextra -ansi,
680/// gcc complains about tile_index maybe being used uninitialised for
681/// oldB in ClusterSequence::_minheap_faster_tiled_N2_cluster(). We
682/// have explicitly checked that it was harmless so we could disable
683/// the gcc warning by hand using the construct below
684///
685/// #pragma GCC diagnostic push
686/// #pragma GCC diagnostic ignored "-Wpragmas"
687/// #pragma GCC diagnostic ignored "-Wuninitialized"
688/// #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
689/// ...
690/// #pragma GCC diagnostic pop
691///
692/// the @GCC diagnostic push/pop directive was only introduced in
693/// gcc-4.6, so for broader usage, we'd need to insert #pragma GCC
694/// diagnostic ignored "-Wpragmas" at the top of this file
695template <class BJ, class I>
696inline void NNFJN2Tiled<BJ,I>::_add_untagged_neighbours_to_tile_union(
697 const int tile_index,
698 int & n_near_tiles) {
699 for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
700 near_tile != _tiles[tile_index].end_tiles; near_tile++){
701 if (! (*near_tile)->tagged) {
702 (*near_tile)->tagged = true;
703 // get the tile number
704 tile_union[n_near_tiles] = *near_tile - & _tiles[0];
705 n_near_tiles++;
706 }
707 }
708}
709
710
711
712FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
713
714
715#endif // __FASTJET_NNFJN2TILED_HH__
Note: See TracBrowser for help on using the repository browser.