FastJet  3.0.6
Dnn3piCylinder.hh
1 //STARTHEADER
2 // $Id: Dnn3piCylinder.hh 2577 2011-09-13 15:11:38Z salam $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 
30 #ifndef DROP_CGAL // in case we do not have the code for CGAL
31 #ifndef __FASTJET_DNN3PICYLINDER_HH__
32 #define __FASTJET_DNN3PICYLINDER_HH__
33 
34 #include "fastjet/internal/DynamicNearestNeighbours.hh"
35 #include "fastjet/internal/DnnPlane.hh"
36 #include "fastjet/internal/numconsts.hh"
37 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 /// \if internal_doc
41 /// @ingroup internal
42 /// \class Dnn3piCylinder
43 /// class derived from DynamicNearestNeighbours that provides an
44 /// implementation for the surface of cylinder (using one
45 /// DnnPlane object spanning 0--3pi).
46 /// \endif
47 class Dnn3piCylinder : public DynamicNearestNeighbours {
48  public:
49  /// empty initaliser
50  Dnn3piCylinder() {}
51 
52  /// Initialiser from a set of points on an Eta-Phi plane, where
53  /// eta can have an arbitrary ranges and phi must be in range
54  /// 0 <= phi < 2pi;
55  ///
56  /// NB: this class is more efficient than the plain Dnn4piCylinder
57  /// class, but can give wrong answers when the nearest neighbour is
58  /// further away than 2pi (in this case a point's nearest neighbour
59  /// becomes itself, because it is considered to be a distance 2pi
60  /// away). For the kt-algorithm (e.g.) this is actually not a
61  /// problem (the distance need only be accurate when it is less than
62  /// R), so we can tell the routine to ignore this problem --
63  /// alternatively the routine will crash if it detects it occurring
64  /// (only when finding the nearest neighbour index, not its
65  /// distance).
66  Dnn3piCylinder(const std::vector<EtaPhi> &,
67  const bool & ignore_nearest_is_mirror = false,
68  const bool & verbose = false );
69 
70  /// Returns the index of the nearest neighbour of point labelled
71  /// by ii (assumes ii is valid)
72  int NearestNeighbourIndex(const int & ii) const ;
73 
74  /// Returns the distance to the nearest neighbour of point labelled
75  /// by index ii (assumes ii is valid)
76  double NearestNeighbourDistance(const int & ii) const ;
77 
78  /// Returns true iff the given index corresponds to a point that
79  /// exists in the DNN structure (meaning that it has been added, and
80  /// not removed in the meantime)
81  bool Valid(const int & index) const;
82 
83  void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
84  const std::vector<EtaPhi> & points_to_add,
85  std::vector<int> & indices_added,
86  std::vector<int> & indices_of_updated_neighbours);
87 
88  ~Dnn3piCylinder();
89 
90  private:
91 
92  // our extras to help us navigate, find distance, etc.
93  const static int INEXISTENT_VERTEX=-3;
94 
95  bool _verbose;
96 
97  bool _ignore_nearest_is_mirror;
98 
99  /// Picture of how things will work... Copy 0--pi part of the 0--2pi
100  /// cylinder into a region 2pi--3pi of a Euclidean plane. Below we
101  /// show points labelled by + that have a mirror image in this
102  /// manner, while points labelled by * do not have a mirror image.
103  ///
104  /// | . |
105  /// | . |
106  /// | . |
107  /// | . |
108  /// | 2 . |
109  /// | * . |
110  /// | + . + |
111  /// | 0 . 1 |
112  /// | . |
113  /// 0 2pi 3pi
114  ///
115  /// Each "true" point has its true "cylinder" index (the index that
116  /// is known externally to this class) as well as euclidean plane
117  /// indices (main_index and mirror index in the MirrorVertexInfo
118  /// structure), which are private concepts of this class.
119  ///
120  /// In above picture our structures would hold the following info
121  /// (the picture shows the euclidean-plane numbering)
122  ///
123  /// _mirror_info[cylinder_index = 0] = (0, 1)
124  /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX)
125  ///
126  /// We also need to be able to go from the euclidean plane indices
127  /// back to the "true" cylinder index, and for this purpose we use
128  /// the vector _cylinder_index_of_plane_vertex[...], which in the above example has
129  /// the following contents
130  ///
131  /// _cylinder_index_of_plane_vertex[0] = 0
132  /// _cylinder_index_of_plane_vertex[1] = 0
133  /// _cylinder_index_of_plane_vertex[2] = 1
134  ///
135 
136  ///
137  struct MirrorVertexInfo {
138  /// index of the given point (appearing in the range 0--2pi) in the
139  /// 0--3pi euclidean plane structure (position will coincide with
140  /// that on the 0--2pi cylinder, but index labelling it will be
141  /// different)
142  int main_index;
143  /// index of the mirror point (appearing in the range 2pi--3pi) in the
144  /// 0--3pi euclidean plane structure
145  int mirror_index;
146  };
147 
148  // for each "true" vertex we have reference to indices in the euclidean
149  // plane structure
150  std::vector<MirrorVertexInfo> _mirror_info;
151  // for each index in the euclidean 0--3pi plane structure we want to
152  // be able to get back to the "true" vertex index on the overall
153  // 0--2pi cylinder structure
154  std::vector<int> _cylinder_index_of_plane_vertex;
155 
156  // NB: we define POINTERS here because the initialisation gave
157  // us problems (things crashed!), perhaps because in practice
158  // we were making a copy without being careful and defining
159  // a proper copy constructor.
160  DnnPlane * _DNN;
161 
162  /// given a phi value in the 0--2pi range return one
163  /// in the pi--3pi range.
164  inline EtaPhi _remap_phi(const EtaPhi & point) {
165  double phi = point.second;
166  if (phi < pi) { phi += twopi ;}
167  return EtaPhi(point.first, phi);}
168 
169 
170  //----------------------------------------------------------------------
171  /// What on earth does this do?
172  ///
173  /// Example: last true "cylinder" index was 15
174  /// last plane index was 23
175  ///
176  /// Then: _cylinder_index_of_plane_vertex.size() = 24 and
177  /// _mirror_info.size() = 16
178  ///
179  /// IF cylinder_point's phi < pi then
180  /// create: _mirror_info[16] = (main_index = 24, mirror_index=25)
181  /// _cylinder_index_of_plane_vertex[24] = 16
182  /// _cylinder_index_of_plane_vertex[25] = 16
183  /// ELSE
184  /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..)
185  /// _cylinder_index_of_plane_vertex[24] = 16
186  ///
187  /// ADDITIONALLY push the cylinder_point (and if it exists the mirror
188  /// copy) onto the vector plane_points.
189  void _RegisterCylinderPoint (const EtaPhi & cylinder_point,
190  std::vector<EtaPhi> & plane_points);
191 };
192 
193 
194 // here follow some inline implementations of the simpler of the
195 // functions defined above
196 
197 //----------------------------------------------------------------------
198 /// Note: one of the difficulties of the 0--3pi mapping is that
199 /// a point may have its mirror copy as its own nearest neighbour
200 /// (if no other point is within a distance of 2pi). This does
201 /// not matter for the kt_algorithm with
202 /// reasonable values of radius, but might matter for a general use
203 /// of this algorithm -- depending on whether or not the user has
204 /// initialised the class with instructions to ignore this problem the
205 /// program will detect and ignore it, or crash.
206 inline int Dnn3piCylinder::NearestNeighbourIndex(const int & current) const {
207  int main_index = _mirror_info[current].main_index;
208  int mirror_index = _mirror_info[current].mirror_index;
209  int plane_index;
210  if (mirror_index == INEXISTENT_VERTEX ) {
211  plane_index = _DNN->NearestNeighbourIndex(main_index);
212  } else {
213  plane_index = (
214  _DNN->NearestNeighbourDistance(main_index) <
215  _DNN->NearestNeighbourDistance(mirror_index)) ?
216  _DNN->NearestNeighbourIndex(main_index) :
217  _DNN->NearestNeighbourIndex(mirror_index) ;
218  }
219  int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index];
220  // either the user has acknowledged the fact that they may get the
221  // mirror copy as the closest point, or crash if it should occur
222  // that mirror copy is the closest point.
223  assert(_ignore_nearest_is_mirror || this_cylinder_index != current);
224  //if (this_cylinder_index == current) {
225  // std::cerr << "WARNING point "<<current<<
226  // " has its mirror copy as its own nearest neighbour"<<endl;
227  //}
228  return this_cylinder_index;
229 }
230 
231 inline double Dnn3piCylinder::NearestNeighbourDistance(const int & current) const {
232  int main_index = _mirror_info[current].main_index;
233  int mirror_index = _mirror_info[current].mirror_index;
234  if (mirror_index == INEXISTENT_VERTEX ) {
235  return _DNN->NearestNeighbourDistance(main_index);
236  } else {
237  return (
238  _DNN->NearestNeighbourDistance(main_index) <
239  _DNN->NearestNeighbourDistance(mirror_index)) ?
240  _DNN->NearestNeighbourDistance(main_index) :
241  _DNN->NearestNeighbourDistance(mirror_index) ;
242  }
243 
244 }
245 
246 inline bool Dnn3piCylinder::Valid(const int & index) const {
247  return (_DNN->Valid(_mirror_info[index].main_index));
248 }
249 
250 
251 inline Dnn3piCylinder::~Dnn3piCylinder() {
252  delete _DNN;
253 }
254 
255 
256 FASTJET_END_NAMESPACE
257 
258 #endif // __FASTJET_DNN3PICYLINDER_HH__
259 #endif // DROP_CGAL