Back to home page

EIC code displayed by LXR

 
 

    


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