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