Back to home page

EIC code displayed by LXR

 
 

    


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