Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-16 10:14:42

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 #include "fastjet/config.h"
0032 
0033 #ifndef DROP_CGAL // in case we do not have the code for CGAL
0034 #ifndef __FASTJET_DNN2PICYLINDER_HH__
0035 #define __FASTJET_DNN2PICYLINDER_HH__
0036 
0037 #include "fastjet/internal/DynamicNearestNeighbours.hh"
0038 #include "fastjet/internal/DnnPlane.hh"
0039 #include "fastjet/internal/numconsts.hh"
0040 
0041 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0042 
0043 
0044 /// \if internal_doc
0045 /// @ingroup internal
0046 /// \class Dnn2piCylinder
0047 /// class derived from DynamicNearestNeighbours that provides an
0048 /// implementation for the surface of cylinder (using one 
0049 /// DnnPlane object spanning 0--2pi).
0050 /// \endif
0051 class Dnn2piCylinder : public DynamicNearestNeighbours {
0052  public:
0053   /// empty initaliser
0054   Dnn2piCylinder() {}
0055 
0056   /// Initialiser from a set of points on an Eta-Phi plane, where
0057   /// eta can have an arbitrary ranges and phi must be in range
0058   /// 0 <= phi < 2pi;
0059   /// 
0060   /// NB: this class is more efficient than the plain Dnn4piCylinder
0061   /// class, but can give wrong answers when the nearest neighbour is
0062   /// further away than 2pi (in this case a point's nearest neighbour
0063   /// becomes itself, because it is considered to be a distance 2pi
0064   /// away). For the kt-algorithm (e.g.) this is actually not a
0065   /// problem (the distance need only be accurate when it is less than
0066   /// R, assuming R<2pi [not necessarily always the case as of
0067   /// 2010-11-19, when we've removed the requirement R<pi/2 in the
0068   /// JetDefinition constructor]), so we can tell the routine to
0069   /// ignore this problem -- alternatively the routine will crash if
0070   /// it detects it occurring (only when finding the nearest neighbour
0071   /// index, not its distance).
0072   Dnn2piCylinder(const std::vector<EtaPhi> &,
0073          const bool & ignore_nearest_is_mirror = false,
0074          const bool & verbose = false );
0075 
0076   /// Returns the index of  the nearest neighbour of point labelled
0077   /// by ii (assumes ii is valid)
0078   int NearestNeighbourIndex(const int ii) const ;
0079 
0080   /// Returns the distance to the nearest neighbour of point labelled
0081   /// by index ii (assumes ii is valid)
0082   double NearestNeighbourDistance(const int ii) const ;
0083 
0084   /// Returns true iff the given index corresponds to a point that
0085   /// exists in the DNN structure (meaning that it has been added, and
0086   /// not removed in the meantime)
0087   bool Valid(const int index) const;
0088 
0089   void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
0090               const std::vector<EtaPhi> & points_to_add,
0091               std::vector<int> & indices_added,
0092               std::vector<int> & indices_of_updated_neighbours);
0093 
0094   ~Dnn2piCylinder();
0095 
0096  private:
0097 
0098   // our extras to help us navigate, find distance, etc.
0099   const static int INEXISTENT_VERTEX=-3;
0100 
0101   bool _verbose;
0102 
0103   bool _ignore_nearest_is_mirror;
0104 
0105   /// Picture of how things will work... Copy 0--pi part of the 0--2pi
0106   /// cylinder into a region 2pi--2pi+ a bit of a Euclidean plane. Below we
0107   /// show points labelled by + that have a mirror image in this
0108   /// manner, while points labelled by * do not have a mirror image.
0109   ///                     
0110   ///      |           .     |        
0111   ///      |           .     |        
0112   ///      |           .     |        
0113   ///      |           .     |        
0114   ///      |        2  .     |        
0115   ///      |        *  .     |        
0116   ///      | +         . +   |        
0117   ///      | 0         . 1   |
0118   ///      |           .     |
0119   ///      0          2pi   2pi + a bit
0120   ///            
0121   /// Each "true" point has its true "cylinder" index (the index that
0122   /// is known externally to this class) as well as euclidean plane
0123   /// indices (main_index and mirror index in the MirrorVertexInfo
0124   /// structure), which are private concepts of this class.
0125   /// 
0126   /// In above picture our structures would hold the following info
0127   /// (the picture shows the euclidean-plane numbering)
0128   ///
0129   /// _mirror_info[cylinder_index = 0] = (0, 1)
0130   /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX)
0131   ///
0132   /// We also need to be able to go from the euclidean plane indices
0133   /// back to the "true" cylinder index, and for this purpose we use
0134   /// the std::vector _cylinder_index_of_plane_vertex[...], which in the above example has
0135   /// the following contents
0136   ///
0137   /// _cylinder_index_of_plane_vertex[0] = 0
0138   /// _cylinder_index_of_plane_vertex[1] = 0
0139   /// _cylinder_index_of_plane_vertex[2] = 1
0140   ///
0141 
0142   /// 
0143   struct MirrorVertexInfo {
0144     /// index of the given point (appearing in the range 0--2pi) in the 
0145     /// 0--2pi euclidean plane structure (position will coincide with
0146     /// that on the 0--2pi cylinder, but index labelling it will be
0147     /// different)
0148     int main_index; 
0149     /// index of the mirror point (appearing in the range 2pi--3pi) in the
0150     /// 0--3pi euclidean plane structure
0151     int mirror_index; 
0152   };
0153 
0154   // for each "true" vertex we have reference to indices in the euclidean
0155   // plane structure
0156   std::vector<MirrorVertexInfo> _mirror_info;
0157   // for each index in the euclidean 0--2pi plane structure we want to
0158   // be able to get back to the "true" vertex index on the overall
0159   // 0--2pi cylinder structure
0160   std::vector<int> _cylinder_index_of_plane_vertex;
0161 
0162   // NB: we define POINTERS here because the initialisation gave
0163   //     us problems (things crashed!), perhaps because in practice
0164   //     we were making a copy without being careful and defining
0165   //     a proper copy constructor.
0166   DnnPlane * _DNN;
0167 
0168   /// given a phi value in the 0--pi range return one 
0169   /// in the 2pi--3pi range; whereas if it is in the pi-2pi range then
0170   /// remap it to be inthe range (-pi)--0.
0171   inline EtaPhi _remap_phi(const EtaPhi & point) {
0172     double phi = point.second;
0173     if (phi < pi) { phi += twopi ;} else {phi -= twopi;}
0174     return EtaPhi(point.first, phi);}
0175 
0176 
0177   //----------------------------------------------------------------------
0178   /// Actions here are similar to those in the
0179   /// Dnn3piCylinder::_RegisterCylinderPoint case, however here we do
0180   /// NOT create the mirror point -- instead we initialise the structure
0181   /// as if there were no need for the mirror point.
0182   ///
0183   /// ADDITIONALLY push the cylinder_point onto the vector plane_points.
0184   void _RegisterCylinderPoint (const EtaPhi & cylinder_point,
0185                    std::vector<EtaPhi> & plane_points);
0186 
0187   /// For each plane point specified in the vector plane_indices,
0188   /// establish whether there is a need to create a mirror point
0189   /// according to the following criteria:
0190   ///
0191   /// . phi < pi
0192   /// . mirror does not already exist
0193   /// . phi < NearestNeighbourDistance 
0194   ///   (if this is not true then there is no way that its mirror point
0195   ///   could have a nearer neighbour).
0196   ///
0197   /// If conditions all hold, then create the mirror point, insert it
0198   /// into the _DNN structure, adjusting any nearest neighbours, and
0199   /// return the list of plane points whose nearest neighbours have
0200   /// changed (this will include the new neighbours that have just been
0201   /// added)
0202   void _CreateNecessaryMirrorPoints(
0203               const std::vector<int> & plane_indices,
0204               std::vector<int> & updated_plane_points);
0205 
0206 };
0207 
0208 
0209 // here follow some inline implementations of the simpler of the
0210 // functions defined above
0211 
0212 //----------------------------------------------------------------------
0213 /// Note: one of the difficulties of the 0--2pi mapping is that
0214 /// a point may have its mirror copy as its own nearest neighbour
0215 /// (if no other point is within a distance of 2pi). This does
0216 /// not matter for the kt_algorithm with
0217 /// reasonable values of radius, but might matter for a general use
0218 /// of this algorithm -- depending on whether or not the user has
0219 /// initialised the class with instructions to ignore this problem the
0220 /// program will detect and ignore it, or crash.
0221 inline int Dnn2piCylinder::NearestNeighbourIndex(const int current) const {
0222   int main_index = _mirror_info[current].main_index;
0223   int mirror_index = _mirror_info[current].mirror_index;
0224   int plane_index;
0225   if (mirror_index == INEXISTENT_VERTEX ) {
0226     plane_index = _DNN->NearestNeighbourIndex(main_index);
0227   } else {
0228     plane_index = (
0229     _DNN->NearestNeighbourDistance(main_index) < 
0230     _DNN->NearestNeighbourDistance(mirror_index)) ? 
0231       _DNN->NearestNeighbourIndex(main_index) : 
0232       _DNN->NearestNeighbourIndex(mirror_index) ; 
0233   }
0234   int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index];
0235   // either the user has acknowledged the fact that they may get the
0236   // mirror copy as the closest point, or crash if it should occur
0237   // that mirror copy is the closest point.
0238   assert(_ignore_nearest_is_mirror || this_cylinder_index != current);
0239   //if (this_cylinder_index == current) {
0240   //  cerr << "WARNING point "<<current<<
0241   //    " has its mirror copy as its own nearest neighbour"<<endl;
0242   //}
0243   return this_cylinder_index;
0244 }
0245 
0246 inline double Dnn2piCylinder::NearestNeighbourDistance(const int current) const {
0247   int main_index = _mirror_info[current].main_index;
0248   int mirror_index = _mirror_info[current].mirror_index;
0249   if (mirror_index == INEXISTENT_VERTEX ) {
0250     return _DNN->NearestNeighbourDistance(main_index);
0251   } else {
0252     return (
0253     _DNN->NearestNeighbourDistance(main_index) < 
0254     _DNN->NearestNeighbourDistance(mirror_index)) ? 
0255       _DNN->NearestNeighbourDistance(main_index) : 
0256       _DNN->NearestNeighbourDistance(mirror_index) ; 
0257   }
0258  
0259 }
0260 
0261 inline bool Dnn2piCylinder::Valid(const int index) const {
0262   return (_DNN->Valid(_mirror_info[index].main_index));
0263 }
0264 
0265 
0266 inline Dnn2piCylinder::~Dnn2piCylinder() {
0267   delete _DNN; 
0268 }
0269 
0270 
0271 FASTJET_END_NAMESPACE
0272 
0273 #endif //  __FASTJET_DNN2PICYLINDER_HH__
0274 #endif //DROP_CGAL