FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
DnnPlane.hh
1 //FJSTARTHEADER
2 // $Id: DnnPlane.hh 3918 2015-07-03 14:19:13Z 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 
34 #ifndef __FASTJET_DNNPLANE_HH__
35 #define __FASTJET_DNNPLANE_HH__
36 
37 #include "fastjet/internal/Triangulation.hh"
38 #include "fastjet/internal/DynamicNearestNeighbours.hh"
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 
43 /// \if internal_doc
44 /// @ingroup internal
45 /// \class DnnPlane
46 /// class derived from DynamicNearestNeighbours that provides an
47 /// implementation for the Euclidean plane
48 ///
49 /// This class that uses CGAL Delaunay triangulation for most of the
50 /// work (it allows for easy and efficient removal and addition of
51 /// points and circulation over a point's neighbours). The treatment
52 /// of coincident points is not supported by CGAL and is implemented
53 /// according to the method specified in
54 /// issue-tracker/2012-02-CGAL-coincident/METHOD
55 /// \endif
56 class DnnPlane : public DynamicNearestNeighbours {
57  public:
58  /// empty initaliser
59  DnnPlane() {}
60 
61  /// Initialiser from a set of points on an Eta-Phi plane, where both
62  /// eta and phi can have arbitrary ranges
63  DnnPlane(const std::vector<EtaPhi> &, const bool & verbose = false );
64 
65 
66  /// Returns the index of the nearest neighbour of point labelled
67  /// by ii (assumes ii is valid)
68  int NearestNeighbourIndex(const int ii) const ;
69 
70  /// Returns the distance to the nearest neighbour of point labelled
71  /// by index ii (assumes ii is valid)
72  double NearestNeighbourDistance(const int ii) const ;
73 
74  /// Returns true iff the given index corresponds to a point that
75  /// exists in the DNN structure (meaning that it has been added, and
76  /// not removed in the meantime)
77  bool Valid(const int index) const;
78 
79  void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
80  const std::vector<EtaPhi> & points_to_add,
81  std::vector<int> & indices_added,
82  std::vector<int> & indices_of_updated_neighbours);
83 
84  /// returns the EtaPhi of point with index i.
85  EtaPhi etaphi(const int i) const;
86  /// returns the eta point with index i.
87  double eta(const int i) const;
88  /// returns the phi point with index i.
89  double phi(const int i) const;
90 
91 private:
92 
93  /// Structure containing a vertex_handle and cached information on
94  /// the nearest neighbour.
95  struct SuperVertex {
96  Vertex_handle vertex; // NULL indicates inexistence...
97  double NNdistance;
98  int NNindex;
99  int coincidence; // ==vertex->info.val() if no coincidence
100  // points to the coinciding SV in case of coincidence
101  // later on for cylinder put a second vertex?
102  };
103 
104  std::vector<SuperVertex> _supervertex;
105  //set<Vertex_handle> _vertex_set;
106  bool _verbose;
107 
108  //static const bool _crash_on_coincidence = true;
109  static const bool _crash_on_coincidence = false;
110 
111  Triangulation _TR; /// CGAL object for dealing with triangulations
112 
113  /// calculates and returns the euclidean distance between points p1
114  /// and p2
115  inline double _euclid_distance(const Point& p1, const Point& p2) const {
116  double distx= p1.x()-p2.x();
117  double disty= p1.y()-p2.y();
118  return distx*distx+disty*disty;
119  }
120 
121  //----------------------------------------------------------------------
122  /// Determines the index and distance of the nearest neighbour to
123  /// point j and puts the information into the _supervertex entry for j
124  void _SetNearest(const int j);
125 
126  //----------------------------------------------------------------------
127  /// Determines and stores the nearest neighbour of j.
128  ///
129  /// For each voronoi neighbour D of j if the distance between j and D
130  /// is less than D's own nearest neighbour, then update the
131  /// nearest-neighbour info in D; push D's index onto
132  /// indices_of_updated_neighbours
133  ///
134  /// Note that j is NOT pushed onto indices_of_updated_neighbours --
135  /// if you want it there, put it there yourself.
136  void _SetAndUpdateNearest(const int j,
137  std::vector<int> & indices_of_updated_neighbours);
138 
139  /// given a vertex_handle returned by CGAL on insertion of a new
140  /// points, returns the coinciding vertex's value if it turns out
141  /// that it corresponds to a vertex that we already knew about
142  /// (usually because two points coincide)
143  int _CheckIfVertexPresent(const Vertex_handle & vertex,
144  const int its_index);
145 
146  //----------------------------------------------------------------------
147  /// if the distance between 'pref' and 'candidate' is smaller (or
148  /// equal) than the one between 'pref' and 'near', return true and
149  /// set 'mindist' to that distance. Note that it is assumed that
150  /// 'mindist' is the euclidian distance between 'pref' and 'near'
151  ///
152  /// Note that the 'near' point is passed through its vertex rather
153  /// than as a point. This allows us to handle cases where we have no min
154  /// yet (near is the infinite vertex)
155  inline bool _is_closer_to(const Point &pref,
156  const Point &candidate,
157  const Vertex_handle &near,
158  double & dist,
159  double & mindist){
160  dist = _euclid_distance(pref, candidate);
161  return _is_closer_to_with_hint(pref, candidate, near, dist, mindist);
162  }
163 
164  /// same as '_is_closer_to' except that 'dist' already contains the
165  /// distance between 'pref' and 'candidate'
166  inline bool _is_closer_to_with_hint(const Point &pref,
167  const Point &candidate,
168  const Vertex_handle &near,
169  const double & dist,
170  double & mindist){
171 
172  // check if 'dist', the pre-computed distance between 'candidate'
173  // and 'pref' is smaller than the distance between 'pref' and its
174  // currently registered nearest neighbour 'near' (and update
175  // things if it is)
176  //
177  // Interestingly enough, it has to be pointed out that the use of
178  // 'abs' instead of 'std::abs' returns wrong results (apparently
179  // ints without any compiler warning)
180  //
181  // The (near != NULL) test is there for one single reason: when
182  // checking that a newly inserted point is not closer than a
183  // previous NN, if that distance comparison involves a "nearly
184  // degenerate" distance we need to access near->point. But
185  // sometimes, in the course of RemoveAndAddPoints, its previous NN
186  // has been deleted and its vertex (corresponding to 'near') set
187  // to NULL. This is not a problem as all points having a deleted
188  // point as NN will have their NN explicitly recomputed at the end
189  // of RemoveAndAddPoints so here we should just make sure there is
190  // no crash... that's done by checking (near != NULL)
191  if ((std::abs(dist-mindist)<DISTANCE_FOR_CGAL_CHECKS) &&
192  _is_not_null(near) && // GPS cf. note below about != NULL with clang/CGAL
193  (_euclid_distance(candidate, near->point())<DISTANCE_FOR_CGAL_CHECKS)){
194  // we're in a situation where there might be a rounding issue,
195  // use CGAL's distance computation to get it right
196  //
197  // Note that in the test right above,
198  // (abs(dist-mindist)<1e-12) guarantees that the current
199  // nearest point is not the infinite vertex and thus
200  // nearest->point() is not ill-defined
201  if (_verbose) std::cout << "using CGAL's distance ordering" << std::endl;
202  if (CGAL::compare_distance_to_point(pref, candidate, near->point())!=CGAL::LARGER){
203  mindist = dist;
204  return true;
205  }
206  } else if (dist <= mindist) {
207  // Note that the use of a <= in the above expression (instead of
208  // a strict ordering <) is important in one case: when checking
209  // if a new point is the new NN of one of the points in its
210  // neighbourhood, in case of distances being ==, we are sure
211  // that 'candidate' is in a cell adjacent to 'pref' while it may
212  // no longer be the case for 'near'
213  mindist = dist;
214  return true;
215  }
216 
217  return false;
218  }
219 
220  /// if a distance between a point and 2 others is smaller than this
221  /// and the distance between the two points is also smaller than this
222  /// then use CGAL to compare the distances.
223  static const double DISTANCE_FOR_CGAL_CHECKS;
224 
225 
226  /// small routine to check if if a vertex handle is not null.
227  ///
228  /// This is centralised here because of the need for a not-very
229  /// readable workaround for certain CGAL/clang++ compiler
230  /// combinations.
231  inline static bool _is_not_null(const Vertex_handle & vh) {
232  // as found in issue 2015-07-clang61+CGAL, the form
233  // vh != NULL
234  // appears to be broken with CGAL-3.6 and clang-3.6.0/.1
235  //
236  // The not very "readable"
237  // vh.operator->()
238  // does the job instead
239  return vh.operator->();
240  }
241 
242 };
243 
244 
245 // here follow some inline implementations of the simpler of the
246 // functions defined above
247 
248 inline int DnnPlane::NearestNeighbourIndex(const int ii) const {
249  return _supervertex[ii].NNindex;}
250 
251 inline double DnnPlane::NearestNeighbourDistance(const int ii) const {
252  return _supervertex[ii].NNdistance;}
253 
254 inline bool DnnPlane::Valid(const int index) const {
255  if (index >= 0 && index < static_cast<int>(_supervertex.size())) {
256  return _is_not_null(_supervertex[index].vertex);
257  }
258  else {
259  return false;
260  }
261 }
262 
263 
264 inline EtaPhi DnnPlane::etaphi(const int i) const {
265  Point * p = & (_supervertex[i].vertex->point());
266  return EtaPhi(p->x(),p->y()); }
267 
268 inline double DnnPlane::eta(const int i) const {
269  return _supervertex[i].vertex->point().x(); }
270 
271 inline double DnnPlane::phi(const int i) const {
272  return _supervertex[i].vertex->point().y(); }
273 
274 
275 FASTJET_END_NAMESPACE
276 
277 #endif // __FASTJET_DNNPLANE_HH__
278 
279 #endif // DROP_CGAL