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