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