FastJet  3.2.2
LazyTiling9Alt.hh
1 #ifndef __FASTJET_LAZYTILING9ALT_HH__
2 #define __FASTJET_LAZYTILING9ALT_HH__
3 
4 //FJSTARTHEADER
5 // $Id: LazyTiling9Alt.hh 3807 2015-02-20 11:16:55Z soyez $
6 //
7 // Copyright (c) 2005-2014, 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/PseudoJet.hh"
35 #include "fastjet/internal/MinHeap.hh"
36 #include "fastjet/ClusterSequence.hh"
37 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 /// Rounding errors in the Lazy strategies may cause the following
41 /// problem: when browsing tiles in the vicinity of the particles
42 /// being clustered in order to decide which of these tiles may
43 /// contain particles that need to be updated (because theit NN is one
44 /// of the particles that are currently clustered), we discard tiles
45 /// that are deemed "too far from the cell" by the "max_NN_dist"
46 /// criterion. Because of rounding error, this condition can sometimes
47 /// miss cases where an update is needed.
48 ///
49 /// An example of this happens if a particle '1' is, say, at the lower
50 /// edge of the rapidity of a given tile, with a particle '2' in the
51 /// tile directly on its left at the same rapidity. Assume also that
52 /// max_NN_dist in 2's tile corresponds to the distance between 2 and
53 /// teh tile of 1. If 2 is 1's NN then in case 2 gets clustered, 1's
54 /// NN needs to be updated. However, rounding errors in the
55 /// calculation of the distance between 1 and 2 may result is
56 /// something slightly larger than the max_NN_dist in 2's tile.
57 ///
58 /// This situation corresponds to the bug reported by Jochen Olt on
59 /// February 12 2015 [see issue-tracker/2015-02-infinite-loop],
60 /// causing an infinite loop.
61 ///
62 /// To prevent this, the simplest solution is, when looking at tiles
63 /// to browse for updateds, to add a margin of security close to the
64 /// edges of the cell, i.e. instead of updating only tiles for which
65 /// distance<=max_NN_dist, we will update tiles for which
66 /// distance<=max_NN_dist+tile_edge_security_margin.
67 ///
68 /// Note that this does not need to be done when computing nearest
69 /// neighbours [rounding errors are tolerated there] but it is
70 /// critical when tracking points that have to be updated.
71 const double tile_edge_security_margin=1.0e-7;
72 
73 /// structure analogous to BriefJet, but with the extra information
74 /// needed for dealing with tiles
75 class TiledJet {
76 public:
77  double eta, phi, kt2, NN_dist;
78  TiledJet * NN, *previous, * next;
79  int _jets_index, tile_index;
80  bool _minheap_update_needed;
81 
82  // indicate whether jets need to have their minheap entries
83  // updated).
84  inline void label_minheap_update_needed() {_minheap_update_needed = true;}
85  inline void label_minheap_update_done() {_minheap_update_needed = false;}
86  inline bool minheap_update_needed() const {return _minheap_update_needed;}
87 };
88 
89 const int n_tile_neighbours = 9;
90 
91 class Tile {
92 public:
93  typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
94  typedef std::pair<Tile *, DistToTileFn> TileFnPair;
95  /// pointers to neighbouring tiles, including self
96  TileFnPair begin_tiles[n_tile_neighbours];
97  /// neighbouring tiles, excluding self
98  TileFnPair * surrounding_tiles;
99  /// half of neighbouring tiles, no self
100  TileFnPair * RH_tiles;
101  /// just beyond end of tiles
102  TileFnPair * end_tiles;
103  /// start of list of BriefJets contained in this tile
104  TiledJet * head;
105  /// sometimes useful to be able to tag a tile
106  bool tagged;
107  /// true for tiles where the delta phi calculation needs
108  /// potentially to account for periodicity in phi
109  bool use_periodic_delta_phi;
110  /// for all particles in the tile, this stores the largest of the
111  /// (squared) nearest-neighbour distances.
112  double max_NN_dist;
113  double eta_min, eta_max, phi_min, phi_max;
114 
115  double distance_to_centre(const TiledJet *) const {return 0;}
116  double distance_to_left(const TiledJet * jet) const {
117  double deta = jet->eta - eta_min;
118  return deta*deta;
119  }
120  double distance_to_right(const TiledJet * jet) const {
121  double deta = jet->eta - eta_max;
122  return deta*deta;
123  }
124  double distance_to_bottom(const TiledJet * jet) const {
125  double dphi = jet->phi - phi_min;
126  return dphi*dphi;
127  }
128  double distance_to_top(const TiledJet * jet) const {
129  double dphi = jet->phi - phi_max;
130  return dphi*dphi;
131  }
132 
133  double distance_to_left_top(const TiledJet * jet) const {
134  double deta = jet->eta - eta_min;
135  double dphi = jet->phi - phi_max;
136  return deta*deta + dphi*dphi;
137  }
138  double distance_to_left_bottom(const TiledJet * jet) const {
139  double deta = jet->eta - eta_min;
140  double dphi = jet->phi - phi_min;
141  return deta*deta + dphi*dphi;
142  }
143  double distance_to_right_top(const TiledJet * jet) const {
144  double deta = jet->eta - eta_max;
145  double dphi = jet->phi - phi_max;
146  return deta*deta + dphi*dphi;
147  }
148  double distance_to_right_bottom(const TiledJet * jet) const {
149  double deta = jet->eta - eta_max;
150  double dphi = jet->phi - phi_min;
151  return deta*deta + dphi*dphi;
152  }
153 
154 
155 };
156 
157 //----------------------------------------------------------------------
158 class LazyTiling9Alt {
159 public:
160  LazyTiling9Alt(ClusterSequence & cs);
161 
162  void run();
163 
164  //void get_next_clustering(int & jetA_index, int & jetB_index, double & dij);
165 
166 
167 protected:
168  ClusterSequence & _cs;
169  const std::vector<PseudoJet> & _jets;
170  std::vector<Tile> _tiles;
171 
172 
173  double _Rparam, _R2, _invR2;
174  double _tiles_eta_min, _tiles_eta_max;
175  double _tile_size_eta, _tile_size_phi;
176  double _tile_half_size_eta, _tile_half_size_phi;
177  int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
178 
179  std::vector<TiledJet *> _jets_for_minheap;
180 
181  //MinHeap _minheap;
182 
183  void _initialise_tiles();
184 
185  // reasonably robust return of tile index given ieta and iphi, in particular
186  // it works even if iphi is negative
187  inline int _tile_index (int ieta, int iphi) const {
188  // note that (-1)%n = -1 so that we have to add _n_tiles_phi
189  // before performing modulo operation
190  return (ieta-_tiles_ieta_min)*_n_tiles_phi
191  + (iphi+_n_tiles_phi) % _n_tiles_phi;
192  }
193 
194  void _bj_remove_from_tiles(TiledJet * const jet);
195 
196  /// returns the tile index given the eta and phi values of a jet
197  int _tile_index(const double eta, const double phi) const;
198 
199  // sets up information regarding the tiling of the given jet
200  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
201 
202  void _print_tiles(TiledJet * briefjets ) const;
203  void _add_neighbours_to_tile_union(const int tile_index,
204  std::vector<int> & tile_union, int & n_near_tiles) const;
205  void _add_untagged_neighbours_to_tile_union(const int tile_index,
206  std::vector<int> & tile_union, int & n_near_tiles);
207  void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet,
208  std::vector<int> & tile_union, int & n_near_tiles);
209  //double _distance_to_tile(const TiledJet * bj, const Tile *) const;
210  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
211 
212  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
213 
214  // return the diJ (multiplied by _R2) for this jet assuming its NN
215  // info is correct
216  template <class J> double _bj_diJ(const J * const jet) const {
217  double kt2 = jet->kt2;
218  if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
219  return jet->NN_dist * kt2;
220  }
221 
222 
223  //----------------------------------------------------------------------
224  template <class J> inline void _bj_set_jetinfo(
225  J * const jetA, const int _jets_index) const {
226  jetA->eta = _jets[_jets_index].rap();
227  jetA->phi = _jets[_jets_index].phi_02pi();
228  jetA->kt2 = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
229  jetA->_jets_index = _jets_index;
230  // initialise NN info as well
231  jetA->NN_dist = _R2;
232  jetA->NN = NULL;
233  }
234 
235 
236  //----------------------------------------------------------------------
237  template <class J> inline double _bj_dist(
238  const J * const jetA, const J * const jetB) const {
239  double dphi = std::abs(jetA->phi - jetB->phi);
240  double deta = (jetA->eta - jetB->eta);
241  if (dphi > pi) {dphi = twopi - dphi;}
242  return dphi*dphi + deta*deta;
243  }
244 
245 
246  //----------------------------------------------------------------------
247  template <class J> inline double _bj_dist_not_periodic(
248  const J * const jetA, const J * const jetB) const {
249  double dphi = jetA->phi - jetB->phi;
250  double deta = (jetA->eta - jetB->eta);
251  return dphi*dphi + deta*deta;
252  }
253 
254 };
255 
256 
257 FASTJET_END_NAMESPACE
258 
259 #endif // __FASTJET_LAZYTILING9ALT_HH__
deals with clustering
const double tile_edge_security_margin
Rounding errors in the Lazy strategies may cause the following problem: when browsing tiles in the vi...
structure analogous to BriefJet, but with the extra information needed for dealing with tiles ...
double jet_scale_for_algorithm(const PseudoJet &jet) const
returns the scale associated with a jet as required for this clustering algorithm (kt^2 for the kt-al...