|
||||
File indexing completed on 2025-01-18 09:57:15
0001 //FJSTARTHEADER 0002 // $Id$ 0003 // 0004 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 0005 // 0006 //---------------------------------------------------------------------- 0007 // This file is part of FastJet. 0008 // 0009 // FastJet is free software; you can redistribute it and/or modify 0010 // it under the terms of the GNU General Public License as published by 0011 // the Free Software Foundation; either version 2 of the License, or 0012 // (at your option) any later version. 0013 // 0014 // The algorithms that underlie FastJet have required considerable 0015 // development. They are described in the original FastJet paper, 0016 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use 0017 // FastJet as part of work towards a scientific publication, please 0018 // quote the version you use and include a citation to the manual and 0019 // optionally also to hep-ph/0512210. 0020 // 0021 // FastJet is distributed in the hope that it will be useful, 0022 // but WITHOUT ANY WARRANTY; without even the implied warranty of 0023 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0024 // GNU General Public License for more details. 0025 // 0026 // You should have received a copy of the GNU General Public License 0027 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 0028 //---------------------------------------------------------------------- 0029 //FJENDHEADER 0030 0031 0032 #ifndef DROP_CGAL // in case we do not have the code for CGAL 0033 0034 #ifndef __FASTJET_DNNPLANE_HH__ 0035 #define __FASTJET_DNNPLANE_HH__ 0036 0037 #include "fastjet/internal/Triangulation.hh" 0038 #include "fastjet/internal/DynamicNearestNeighbours.hh" 0039 0040 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 0041 0042 0043 /// \if internal_doc 0044 /// @ingroup internal 0045 /// \class DnnPlane 0046 /// class derived from DynamicNearestNeighbours that provides an 0047 /// implementation for the Euclidean plane 0048 /// 0049 /// This class that uses CGAL Delaunay triangulation for most of the 0050 /// work (it allows for easy and efficient removal and addition of 0051 /// points and circulation over a point's neighbours). The treatment 0052 /// of coincident points is not supported by CGAL and is implemented 0053 /// according to the method specified in 0054 /// issue-tracker/2012-02-CGAL-coincident/METHOD 0055 /// \endif 0056 class DnnPlane : public DynamicNearestNeighbours { 0057 public: 0058 /// empty initaliser 0059 DnnPlane() {} 0060 0061 /// Initialiser from a set of points on an Eta-Phi plane, where both 0062 /// eta and phi can have arbitrary ranges 0063 DnnPlane(const std::vector<EtaPhi> &, const bool & verbose = false ); 0064 0065 0066 /// Returns the index of the nearest neighbour of point labelled 0067 /// by ii (assumes ii is valid) 0068 int NearestNeighbourIndex(const int ii) const ; 0069 0070 /// Returns the distance to the nearest neighbour of point labelled 0071 /// by index ii (assumes ii is valid) 0072 double NearestNeighbourDistance(const int ii) const ; 0073 0074 /// Returns true iff the given index corresponds to a point that 0075 /// exists in the DNN structure (meaning that it has been added, and 0076 /// not removed in the meantime) 0077 bool Valid(const int index) const; 0078 0079 void RemoveAndAddPoints(const std::vector<int> & indices_to_remove, 0080 const std::vector<EtaPhi> & points_to_add, 0081 std::vector<int> & indices_added, 0082 std::vector<int> & indices_of_updated_neighbours); 0083 0084 /// returns the EtaPhi of point with index i. 0085 EtaPhi etaphi(const int i) const; 0086 /// returns the eta point with index i. 0087 double eta(const int i) const; 0088 /// returns the phi point with index i. 0089 double phi(const int i) const; 0090 0091 private: 0092 0093 /// Structure containing a vertex_handle and cached information on 0094 /// the nearest neighbour. 0095 struct SuperVertex { 0096 Vertex_handle vertex; // NULL indicates inexistence... 0097 double NNdistance; 0098 int NNindex; 0099 int coincidence; // ==vertex->info.val() if no coincidence 0100 // points to the coinciding SV in case of coincidence 0101 // later on for cylinder put a second vertex? 0102 }; 0103 0104 std::vector<SuperVertex> _supervertex; 0105 //set<Vertex_handle> _vertex_set; 0106 bool _verbose; 0107 0108 //static const bool _crash_on_coincidence = true; 0109 static const bool _crash_on_coincidence = false; 0110 0111 Triangulation _TR; /// CGAL object for dealing with triangulations 0112 0113 /// calculates and returns the euclidean distance between points p1 0114 /// and p2 0115 inline double _euclid_distance(const Point& p1, const Point& p2) const { 0116 double distx= p1.x()-p2.x(); 0117 double disty= p1.y()-p2.y(); 0118 return distx*distx+disty*disty; 0119 } 0120 0121 //---------------------------------------------------------------------- 0122 /// Determines the index and distance of the nearest neighbour to 0123 /// point j and puts the information into the _supervertex entry for j 0124 void _SetNearest(const int j); 0125 0126 //---------------------------------------------------------------------- 0127 /// Determines and stores the nearest neighbour of j. 0128 /// 0129 /// For each voronoi neighbour D of j if the distance between j and D 0130 /// is less than D's own nearest neighbour, then update the 0131 /// nearest-neighbour info in D; push D's index onto 0132 /// indices_of_updated_neighbours 0133 /// 0134 /// Note that j is NOT pushed onto indices_of_updated_neighbours -- 0135 /// if you want it there, put it there yourself. 0136 void _SetAndUpdateNearest(const int j, 0137 std::vector<int> & indices_of_updated_neighbours); 0138 0139 /// given a vertex_handle returned by CGAL on insertion of a new 0140 /// points, returns the coinciding vertex's value if it turns out 0141 /// that it corresponds to a vertex that we already knew about 0142 /// (usually because two points coincide) 0143 int _CheckIfVertexPresent(const Vertex_handle & vertex, 0144 const int its_index); 0145 0146 //---------------------------------------------------------------------- 0147 /// if the distance between 'pref' and 'candidate' is smaller (or 0148 /// equal) than the one between 'pref' and 'near', return true and 0149 /// set 'mindist' to that distance. Note that it is assumed that 0150 /// 'mindist' is the euclidian distance between 'pref' and 'near' 0151 /// 0152 /// Note that the 'near' point is passed through its vertex rather 0153 /// than as a point. This allows us to handle cases where we have no min 0154 /// yet (near is the infinite vertex) 0155 inline bool _is_closer_to(const Point &pref, 0156 const Point &candidate, 0157 const Vertex_handle &near, 0158 double & dist, 0159 double & mindist){ 0160 dist = _euclid_distance(pref, candidate); 0161 return _is_closer_to_with_hint(pref, candidate, near, dist, mindist); 0162 } 0163 0164 /// same as '_is_closer_to' except that 'dist' already contains the 0165 /// distance between 'pref' and 'candidate' 0166 inline bool _is_closer_to_with_hint(const Point &pref, 0167 const Point &candidate, 0168 const Vertex_handle &near, 0169 const double & dist, 0170 double & mindist){ 0171 0172 // check if 'dist', the pre-computed distance between 'candidate' 0173 // and 'pref' is smaller than the distance between 'pref' and its 0174 // currently registered nearest neighbour 'near' (and update 0175 // things if it is) 0176 // 0177 // Interestingly enough, it has to be pointed out that the use of 0178 // 'abs' instead of 'std::abs' returns wrong results (apparently 0179 // ints without any compiler warning) 0180 // 0181 // The (near != NULL) test is there for one single reason: when 0182 // checking that a newly inserted point is not closer than a 0183 // previous NN, if that distance comparison involves a "nearly 0184 // degenerate" distance we need to access near->point. But 0185 // sometimes, in the course of RemoveAndAddPoints, its previous NN 0186 // has been deleted and its vertex (corresponding to 'near') set 0187 // to NULL. This is not a problem as all points having a deleted 0188 // point as NN will have their NN explicitly recomputed at the end 0189 // of RemoveAndAddPoints so here we should just make sure there is 0190 // no crash... that's done by checking (near != NULL) 0191 if ((std::abs(dist-mindist)<DISTANCE_FOR_CGAL_CHECKS) && 0192 _is_not_null(near) && // GPS cf. note below about != NULL with clang/CGAL 0193 (_euclid_distance(candidate, near->point())<DISTANCE_FOR_CGAL_CHECKS)){ 0194 // we're in a situation where there might be a rounding issue, 0195 // use CGAL's distance computation to get it right 0196 // 0197 // Note that in the test right above, 0198 // (abs(dist-mindist)<1e-12) guarantees that the current 0199 // nearest point is not the infinite vertex and thus 0200 // nearest->point() is not ill-defined 0201 if (_verbose) std::cout << "using CGAL's distance ordering" << std::endl; 0202 if (CGAL::compare_distance_to_point(pref, candidate, near->point())!=CGAL::LARGER){ 0203 mindist = dist; 0204 return true; 0205 } 0206 } else if (dist <= mindist) { 0207 // Note that the use of a <= in the above expression (instead of 0208 // a strict ordering <) is important in one case: when checking 0209 // if a new point is the new NN of one of the points in its 0210 // neighbourhood, in case of distances being ==, we are sure 0211 // that 'candidate' is in a cell adjacent to 'pref' while it may 0212 // no longer be the case for 'near' 0213 mindist = dist; 0214 return true; 0215 } 0216 0217 return false; 0218 } 0219 0220 /// if a distance between a point and 2 others is smaller than this 0221 /// and the distance between the two points is also smaller than this 0222 /// then use CGAL to compare the distances. 0223 static const double DISTANCE_FOR_CGAL_CHECKS; 0224 0225 0226 /// small routine to check if if a vertex handle is not null. 0227 /// 0228 /// This is centralised here because of the need for a not-very 0229 /// readable workaround for certain CGAL/clang++ compiler 0230 /// combinations. 0231 inline static bool _is_not_null(const Vertex_handle & vh) { 0232 // as found in issue 2015-07-clang61+CGAL, the form 0233 // vh != NULL 0234 // appears to be broken with CGAL-3.6 and clang-3.6.0/.1 0235 // 0236 // The not very "readable" 0237 // vh.operator->() 0238 // does the job instead 0239 return vh.operator->(); 0240 } 0241 0242 }; 0243 0244 0245 // here follow some inline implementations of the simpler of the 0246 // functions defined above 0247 0248 inline int DnnPlane::NearestNeighbourIndex(const int ii) const { 0249 return _supervertex[ii].NNindex;} 0250 0251 inline double DnnPlane::NearestNeighbourDistance(const int ii) const { 0252 return _supervertex[ii].NNdistance;} 0253 0254 inline bool DnnPlane::Valid(const int index) const { 0255 if (index >= 0 && index < static_cast<int>(_supervertex.size())) { 0256 return _is_not_null(_supervertex[index].vertex); 0257 } 0258 else { 0259 return false; 0260 } 0261 } 0262 0263 0264 inline EtaPhi DnnPlane::etaphi(const int i) const { 0265 Point * p = & (_supervertex[i].vertex->point()); 0266 return EtaPhi(p->x(),p->y()); } 0267 0268 inline double DnnPlane::eta(const int i) const { 0269 return _supervertex[i].vertex->point().x(); } 0270 0271 inline double DnnPlane::phi(const int i) const { 0272 return _supervertex[i].vertex->point().y(); } 0273 0274 0275 FASTJET_END_NAMESPACE 0276 0277 #endif // __FASTJET_DNNPLANE_HH__ 0278 0279 #endif // DROP_CGAL
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |