FastJet  3.3.2
NNFJN2Tiled.hh
1 #ifndef __FASTJET_NNFJN2TILED_HH__
2 #define __FASTJET_NNFJN2TILED_HH__
3 
4 //FJSTARTHEADER
5 // $Id: NNFJN2Tiled.hh 4355 2018-04-22 15:38:54Z salam $
6 //
7 // Copyright (c) 2016-2018, 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 
37 FASTJET_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 ///
113 template<class BJ, class I = _NoInfo> class NNFJN2Tiled : public NNBase<I> {
114 public:
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
138  delete[] briefjets;
139  delete[] diJ;
140  }
141 
142 private:
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  class diJ_plus_link {
261  public:
262  double diJ; // the distance
263  TiledJet * jet; // the jet (i) for which we've found this distance
264  // (whose NN will the J).
265  };
266 
267 };
268 
269 
270 
271 //----------------------------------------------------------------------
272 template<class BJ, class I> void NNFJN2Tiled<BJ,I>::start(const std::vector<PseudoJet> & jets) {
273 
274  _initialise_tiles(jets);
275 
276  n = jets.size();
277 
278  briefjets = new TiledJet[n];
279  where_is.resize(2*n);
280 
281  TiledJet * jetA = briefjets, * jetB;
282 
283  // will be used quite deep inside loops, but declare it here so that
284  // memory (de)allocation gets done only once
285  tile_union.resize(3*n_tile_neighbours);
286 
287  // initialise the basic jet info
288  for (int i = 0; i< n; i++) {
289  _tiledjet_set_jetinfo(jetA, jets[i], i);
290  where_is[i] = jetA;
291  jetA++; // move on to next entry of briefjets
292  }
293 
294  head = briefjets; // a nicer way of naming start
295 
296  // set up the initial nearest neighbour information
297  typename std::vector<Tile>::const_iterator tile;
298  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
299  // first do it on this tile
300  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
301  for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
302  double dist = jetA->geometrical_distance(jetB);
303  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
304  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
305  }
306  }
307  // then do it for RH tiles
308  for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
309  for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
310  for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
311  double dist = jetA->geometrical_distance(jetB);
312  if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
313  if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
314  }
315  }
316  }
317  // no need to do it for LH tiles, since they are implicitly done
318  // when we set NN for both jetA and jetB on the RH tiles.
319  }
320 
321  diJ = new diJ_plus_link[n];
322  jetA = head;
323  for (int i = 0; i < n; i++) {
324  diJ[i].diJ = _compute_diJ(jetA); // kt distance * R^2
325  diJ[i].jet = jetA; // our compact diJ table will not be in
326  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
327  // so set up bi-directional correspondence here.
328  jetA++; // have jetA follow i
329  }
330 }
331 
332 
333 //----------------------------------------------------------------------
334 template<class BJ, class I> double NNFJN2Tiled<BJ,I>::dij_min(int & iA, int & iB) {
335  // find the minimum of the diJ on this round
336  diJ_plus_link * best, *stop; // pointers a bit faster than indices
337  // could use best to keep track of diJ
338  // min, but it turns out to be
339  // marginally faster to have a separate
340  // variable (avoids n dereferences at
341  // the expense of n/2 assignments).
342  double diJ_min = diJ[0].diJ; // initialise the best one here.
343  best = diJ; // and here
344  stop = diJ+n;
345  for (diJ_plus_link * here = diJ+1; here != stop; here++) {
346  if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
347  }
348 
349  // return information to user about recombination
350  TiledJet * jetA = best->jet;
351  iA = jetA->jet_index();
352  iB = jetA->NN ? jetA->NN->jet_index() : -1;
353  return diJ_min;
354 }
355 
356 
357 //----------------------------------------------------------------------
358 // remove jetA from the list
359 template<class BJ, class I> void NNFJN2Tiled<BJ,I>::remove_jet(int iA) {
360  TiledJet * jetA = where_is[iA];
361 
362  _bj_remove_from_tiles(jetA);
363 
364  // first establish the set of tiles over which we are going to
365  // have to run searches for updated and new nearest-neighbours --
366  // basically a combination of vicinity of the tiles of the two old
367  // and one new jet.
368  int n_near_tiles = 0;
369  _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
370 
371  // now update our nearest neighbour info and diJ table
372  // first reduce size of table
373  n--;
374  // then compactify the diJ by taking the last of the diJ and copying
375  // it to the position occupied by the diJ for jetA
376  diJ[n].jet->diJ_posn = jetA->diJ_posn;
377  diJ[jetA->diJ_posn] = diJ[n];
378 
379  // updating other particles' NN.
380  // Run over all tiles in our union
381  for (int itile = 0; itile < n_near_tiles; itile++) {
382  Tile * tile_ptr = &_tiles[tile_union[itile]];
383  tile_ptr->tagged = false; // reset tag, since we're done with unions
384  // run over all jets in the current tile
385  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
386  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
387  if (jetI->NN == jetA) {
388  jetI->NN_dist = jetI->geometrical_beam_distance();
389  jetI->NN = NULL;
390  // now go over tiles that are neighbours of I (include own tile)
391  for (Tile ** near_tile = tile_ptr->begin_tiles;
392  near_tile != tile_ptr->end_tiles; near_tile++) {
393  // and then over the contents of that tile
394  for (TiledJet * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
395  double dist = jetI->geometrical_distance(jetJ);
396  if (dist < jetI->NN_dist && jetJ != jetI) {
397  jetI->NN_dist = dist; jetI->NN = jetJ;
398  }
399  }
400  }
401  diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
402  }
403  }
404  }
405 
406 }
407 
408 
409 //----------------------------------------------------------------------
410 template<class BJ, class I> void NNFJN2Tiled<BJ,I>::merge_jets(int iA, int iB,
411  const PseudoJet & jet, int index) {
412 
413  TiledJet * jetA = where_is[iA];
414  TiledJet * jetB = where_is[iB];
415 
416  // jet-jet recombination
417  // If necessary relabel A & B to ensure jetB < jetA, that way if
418  // the larger of them == newtail then that ends up being jetA and
419  // the new jet that is added as jetB is inserted in a position that
420  // has a future!
421  if (jetA < jetB) {std::swap(jetA,jetB);}
422 
423  // what was jetB will now become the new jet
424  _bj_remove_from_tiles(jetA);
425  TiledJet oldB = * jetB; // take a copy because we will need it...
426  _bj_remove_from_tiles(jetB);
427  _tiledjet_set_jetinfo(jetB, jet, index); // cause jetB to become _jets[nn]
428  // (also registers the jet in the tiling)
429  where_is[index] = jetB;
430 
431  // first establish the set of tiles over which we are going to
432  // have to run searches for updated and new nearest-neighbours --
433  // basically a combination of vicinity of the tiles of the two old
434  // and one new jet.
435  int n_near_tiles = 0;
436  _add_untagged_neighbours_to_tile_union(jetA->tile_index, n_near_tiles);
437  if (jetB->tile_index != jetA->tile_index) {
438  _add_untagged_neighbours_to_tile_union(jetB->tile_index, n_near_tiles);
439  }
440  if (oldB.tile_index != jetA->tile_index &&
441  oldB.tile_index != jetB->tile_index) {
442  _add_untagged_neighbours_to_tile_union(oldB.tile_index, n_near_tiles);
443  }
444 
445  // now update our nearest neighbour info and diJ table
446  // first reduce size of table
447  n--;
448  // then compactify the diJ by taking the last of the diJ and copying
449  // it to the position occupied by the diJ for jetA
450  diJ[n].jet->diJ_posn = jetA->diJ_posn;
451  diJ[jetA->diJ_posn] = diJ[n];
452 
453  // Initialise jetB's NN distance as well as updating it for
454  // other particles.
455  // Run over all tiles in our union
456  for (int itile = 0; itile < n_near_tiles; itile++) {
457  Tile * tile_ptr = &_tiles[tile_union[itile]];
458  tile_ptr->tagged = false; // reset tag, since we're done with unions
459  // run over all jets in the current tile
460  for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
461  // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
462  if ((jetI->NN == jetA) || (jetI->NN == jetB)) {
463  jetI->NN_dist = jetI->geometrical_beam_distance();
464  jetI->NN = NULL;
465  // now go over tiles that are neighbours of I (include own tile)
466  for (Tile ** near_tile = tile_ptr->begin_tiles; near_tile != tile_ptr->end_tiles; near_tile++) {
467  // and then over the contents of that tile
468  for (TiledJet * jetJ = (*near_tile)->head; jetJ != NULL; jetJ = jetJ->next) {
469  double dist = jetI->geometrical_distance(jetJ);
470  if (dist < jetI->NN_dist && jetJ != jetI) {
471  jetI->NN_dist = dist; jetI->NN = jetJ;
472  }
473  }
474  }
475  diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ kt-dist
476  }
477  // check whether new jetB is closer than jetI's current NN and
478  // if jetI is closer than jetB's current (evolving) nearest
479  // neighbour. Where relevant update things
480  double dist = jetI->geometrical_distance(jetB);
481  if (dist < jetI->NN_dist) {
482  if (jetI != jetB) {
483  jetI->NN_dist = dist;
484  jetI->NN = jetB;
485  diJ[jetI->diJ_posn].diJ = _compute_diJ(jetI); // update diJ...
486  }
487  }
488  if (dist < jetB->NN_dist) {
489  if (jetI != jetB) {
490  jetB->NN_dist = dist;
491  jetB->NN = jetI;}
492  }
493  }
494  }
495 
496  // finally, register the updated kt distance for B
497  diJ[jetB->diJ_posn].diJ = _compute_diJ(jetB);
498 }
499 
500 
501 //----------------------------------------------------------------------
502 /// Set up the tiles:
503 /// - decide the range in eta
504 /// - allocate the tiles
505 /// - set up the cross-referencing info between tiles
506 ///
507 /// The neighbourhood of a tile is set up as follows
508 ///
509 /// LRR
510 /// LXR
511 /// LLR
512 ///
513 /// such that tiles is an array containing XLLLLRRRR with pointers
514 /// | \ RH_tiles
515 /// \ surrounding_tiles
516 ///
517 /// with appropriate precautions when close to the edge of the tiled
518 /// region.
519 ///
520 template <class BJ, class I>
521 void NNFJN2Tiled<BJ,I>::_initialise_tiles(const std::vector<PseudoJet> &particles) {
522 
523  // first decide tile sizes (with a lower bound to avoid huge memory use with
524  // very small R)
525  double default_size = _requested_tile_size>0.1 ? _requested_tile_size : 0.1;
526  _tile_size_rap = default_size;
527  // it makes no sense to go below 3 tiles in phi -- 3 tiles is
528  // sufficient to make sure all pair-wise combinations up to pi in
529  // phi are possible
530  _n_tiles_phi = int(floor(twopi/default_size));
531  if (_n_tiles_phi<3) _n_tiles_phi = 3;
532  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
533 
534  TilingExtent tiling_analysis(particles);
535  _tiles_rap_min = tiling_analysis.minrap();
536  _tiles_rap_max = tiling_analysis.maxrap();
537 
538  // now adjust the values
539  _tiles_irap_min = int(floor(_tiles_rap_min/_tile_size_rap));
540  _tiles_irap_max = int(floor( _tiles_rap_max/_tile_size_rap));
541  _tiles_rap_min = _tiles_irap_min * _tile_size_rap;
542  _tiles_rap_max = _tiles_irap_max * _tile_size_rap;
543 
544  // allocate the tiles
545  _tiles.resize((_tiles_irap_max-_tiles_irap_min+1)*_n_tiles_phi);
546 
547  // now set up the cross-referencing between tiles
548  for (int irap = _tiles_irap_min; irap <= _tiles_irap_max; irap++) {
549  for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
550  Tile * tile = & _tiles[_tile_index(irap,iphi)];
551  // no jets in this tile yet
552  tile->head = NULL; // first element of tiles points to itself
553  tile->begin_tiles[0] = tile;
554  Tile ** pptile = & (tile->begin_tiles[0]);
555  pptile++;
556  //
557  // set up L's in column to the left of X
558  tile->surrounding_tiles = pptile;
559  if (irap > _tiles_irap_min) {
560  // with the itile subroutine, we can safely run tiles from
561  // idphi=-1 to idphi=+1, because it takes care of
562  // negative and positive boundaries
563  for (int idphi = -1; idphi <=+1; idphi++) {
564  *pptile = & _tiles[_tile_index(irap-1,iphi+idphi)];
565  pptile++;
566  }
567  }
568  // now set up last L (below X)
569  *pptile = & _tiles[_tile_index(irap,iphi-1)];
570  pptile++;
571  // set up first R (above X)
572  tile->RH_tiles = pptile;
573  *pptile = & _tiles[_tile_index(irap,iphi+1)];
574  pptile++;
575  // set up remaining R's, to the right of X
576  if (irap < _tiles_irap_max) {
577  for (int idphi = -1; idphi <= +1; idphi++) {
578  *pptile = & _tiles[_tile_index(irap+1,iphi+idphi)];
579  pptile++;
580  }
581  }
582  // now put semaphore for end tile
583  tile->end_tiles = pptile;
584  // finally make sure tiles are untagged
585  tile->tagged = false;
586  }
587  }
588 
589 }
590 
591 //----------------------------------------------------------------------
592 /// return the tile index corresponding to the given rap,phi point
593 template <class BJ, class I>
594 int NNFJN2Tiled<BJ,I>::_tile_index(const double rap, const double phi) const {
595  int irap, iphi;
596  if (rap <= _tiles_rap_min) {irap = 0;}
597  else if (rap >= _tiles_rap_max) {irap = _tiles_irap_max-_tiles_irap_min;}
598  else {
599  //irap = int(floor((rap - _tiles_rap_min) / _tile_size_rap));
600  irap = int(((rap - _tiles_rap_min) / _tile_size_rap));
601  // following needed in case of rare but nasty rounding errors
602  if (irap > _tiles_irap_max-_tiles_irap_min) {
603  irap = _tiles_irap_max-_tiles_irap_min;}
604  }
605  // allow for some extent of being beyond range in calculation of phi
606  // as well
607  //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
608  // with just int and no floor, things run faster but beware
609  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
610  return (iphi + irap * _n_tiles_phi);
611 }
612 
613 //----------------------------------------------------------------------
614 template <class BJ, class I>
615 void NNFJN2Tiled<BJ,I>::_bj_remove_from_tiles(TiledJet * const jet) {
616  Tile * tile = & _tiles[jet->tile_index];
617 
618  if (jet->previous == NULL) {
619  // we are at head of the tile, so reset it.
620  // If this was the only jet on the tile then tile->head will now be NULL
621  tile->head = jet->next;
622  } else {
623  // adjust link from previous jet in this tile
624  jet->previous->next = jet->next;
625  }
626  if (jet->next != NULL) {
627  // adjust backwards-link from next jet in this tile
628  jet->next->previous = jet->previous;
629  }
630 }
631 
632 
633 //----------------------------------------------------------------------
634 // overloaded version which additionally sets up information regarding the
635 // tiling
636 template <class BJ, class I>
637 inline void NNFJN2Tiled<BJ,I>::_tiledjet_set_jetinfo(TiledJet * const tile_jet,
638  const PseudoJet &jet,
639  int index) {
640  // the this-> in the next line is required by standard compiler
641  // see e.g. http://stackoverflow.com/questions/10639053/name-lookups-in-c-templates
642  this->init_jet(tile_jet, jet, index);
643 
644  // Then do the setup specific to the tiled case.
645 
646  // Find out which tile it belonds to
647  tile_jet->tile_index = _tile_index(tile_jet->rap(), tile_jet->phi());
648 
649  // Insert it into the tile's linked list of jets
650  Tile * tile = &_tiles[tile_jet->tile_index];
651  tile_jet->previous = NULL;
652  tile_jet->next = tile->head;
653  if (tile_jet->next != NULL) {tile_jet->next->previous = tile_jet;}
654  tile->head = tile_jet;
655 }
656 
657 //----------------------------------------------------------------------
658 /// Add to the vector tile_union the tiles that are in the neighbourhood
659 /// of the specified tile_index, including itself -- start adding
660 /// from position n_near_tiles-1, and increase n_near_tiles as
661 /// you go along (could have done it more C++ like with vector with reserved
662 /// space, but fear is that it would have been slower, e.g. checking
663 /// for end of vector at each stage to decide whether to resize it)
664 template <class BJ, class I>
665 void NNFJN2Tiled<BJ,I>::_add_neighbours_to_tile_union(const int tile_index,
666  int & n_near_tiles) const {
667  for (Tile * const * near_tile = _tiles[tile_index].begin_tiles;
668  near_tile != _tiles[tile_index].end_tiles; near_tile++){
669  // get the tile number
670  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
671  n_near_tiles++;
672  }
673 }
674 
675 //----------------------------------------------------------------------
676 /// Like _add_neighbours_to_tile_union, but only adds neighbours if
677 /// their "tagged" status is false; when a neighbour is added its
678 /// tagged status is set to true.
679 ///
680 /// Note that with a high level of warnings (-pedantic -Wextra -ansi,
681 /// gcc complains about tile_index maybe being used uninitialised for
682 /// oldB in ClusterSequence::_minheap_faster_tiled_N2_cluster(). We
683 /// have explicitly checked that it was harmless so we could disable
684 /// the gcc warning by hand using the construct below
685 ///
686 /// #pragma GCC diagnostic push
687 /// #pragma GCC diagnostic ignored "-Wpragmas"
688 /// #pragma GCC diagnostic ignored "-Wuninitialized"
689 /// #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
690 /// ...
691 /// #pragma GCC diagnostic pop
692 ///
693 /// the @GCC diagnostic push/pop directive was only introduced in
694 /// gcc-4.6, so for broader usage, we'd need to insert #pragma GCC
695 /// diagnostic ignored "-Wpragmas" at the top of this file
696 template <class BJ, class I>
698  const int tile_index,
699  int & n_near_tiles) {
700  for (Tile ** near_tile = _tiles[tile_index].begin_tiles;
701  near_tile != _tiles[tile_index].end_tiles; near_tile++){
702  if (! (*near_tile)->tagged) {
703  (*near_tile)->tagged = true;
704  // get the tile number
705  tile_union[n_near_tiles] = *near_tile - & _tiles[0];
706  n_near_tiles++;
707  }
708  }
709 }
710 
711 
712 
713 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
714 
715 
716 #endif // __FASTJET_NNFJN2TILED_HH__
Helps solve closest pair problems with factorised interparticle and beam distances (ie satisfying the...
Definition: NNFJN2Tiled.hh:113
class to perform a fast analysis of the appropriate rapidity range in which to perform tiling ...
Definition: TilingExtent.hh:41
NNFJN2Tiled(const std::vector< PseudoJet > &jets, double requested_tile_size)
constructor with an initial set of jets (which will be assigned indices 0...jets.size()-1) ...
Definition: NNFJN2Tiled.hh:118
~NNFJN2Tiled()
a destructor
Definition: NNFJN2Tiled.hh:137
double minrap() const
returns the suggested minimum rapidity for the tiling
Definition: TilingExtent.hh:52
Helps solve closest pair problems with generic interparticle and particle-beam distances.
Definition: NNBase.hh:164
double maxrap() const
returns the suggested maximum rapidity for the tiling
Definition: TilingExtent.hh:54
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67