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