JETAN/fastjet/fastjet/internal/Dnn2piCylinder.hh

   1 //STARTHEADER
   2 // $Id: Dnn2piCylinder.hh 431 2007-01-20 10:44:55Z salam $
   3 //
   4 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
   5 //
   6 //----------------------------------------------------------------------
   7 // This file is part of FastJet.
   8 //
   9 //  FastJet is free software; you can redistribute it and/or modify
  10 //  it under the terms of the GNU General Public License as published by
  11 //  the Free Software Foundation; either version 2 of the License, or
  12 //  (at your option) any later version.
  13 //
  14 //  The algorithms that underlie FastJet have required considerable
  15 //  development and are described in hep-ph/0512210. If you use
  16 //  FastJet as part of work towards a scientific publication, please
  17 //  include a citation to the FastJet paper.
  18 //
  19 //  FastJet is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 //  GNU General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU General Public License
  25 //  along with FastJet; if not, write to the Free Software
  26 //  Foundation, Inc.:
  27 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  28 //----------------------------------------------------------------------
  29 //ENDHEADER
  30
  31
  32 #ifndef DROP_CGAL // in case we do not have the code for CGAL
  33 #ifndef __FASTJET_DNN2PICYLINDER_HH__
  34 #define __FASTJET_DNN2PICYLINDER_HH__
  35
  36 #include "fastjet/internal/DynamicNearestNeighbours.hh"
  37 #include "fastjet/internal/DnnPlane.hh"
  38 #include "fastjet/internal/numconsts.hh"
  39
  40 namespace fastjet {      // defined in fastjet/internal/base.hh
  41
  42
  43 /// class derived from DynamicNearestNeighbours that provides an
  44 /// implementation for the surface of cylinder (using one
  45 /// DnnPlane object spanning 0--2pi).
  46 class Dnn2piCylinder : public DynamicNearestNeighbours {
  47  public:
  48   /// empty initaliser
  49   Dnn2piCylinder() {}
  50
  51   /// Initialiser from a set of points on an Eta-Phi plane, where
  52   /// eta can have an arbitrary ranges and phi must be in range
  53   /// 0 <= phi < 2pi;
  54   ///
  55   /// NB: this class is more efficient than the plain Dnn4piCylinder
  56   /// class, but can give wrong answers when the nearest neighbour is
  57   /// further away than 2pi (in this case a point's nearest neighbour
  58   /// becomes itself, because it is considered to be a distance 2pi
  59   /// away). For the kt-algorithm (e.g.) this is actually not a
  60   /// problem (the distance need only be accurate when it is less than
  61   /// R), so we can tell the routine to ignore this problem --
  62   /// alternatively the routine will crash if it detects it occurring
  63   /// (only when finding the nearest neighbour index, not its
  64   /// distance).
  65   Dnn2piCylinder(const std::vector<EtaPhi> &,
  66                  const bool & ignore_nearest_is_mirror = false,
  67                  const bool & verbose = false );
  68
  69   /// Returns the index of  the nearest neighbour of point labelled
  70   /// by ii (assumes ii is valid)
  71   int NearestNeighbourIndex(const int & ii) const ;
  72
  73   /// Returns the distance to the nearest neighbour of point labelled
  74   /// by index ii (assumes ii is valid)
  75   double NearestNeighbourDistance(const int & ii) const ;
  76
  77   /// Returns true iff the given index corresponds to a point that
  78   /// exists in the DNN structure (meaning that it has been added, and
  79   /// not removed in the meantime)
  80   bool Valid(const int & index) const;
  81
  82   void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
  83                           const std::vector<EtaPhi> & points_to_add,
  84                           std::vector<int> & indices_added,
  85                           std::vector<int> & indices_of_updated_neighbours);
  86
  87   ~Dnn2piCylinder();
  88
  89  private:
  90
  91   // our extras to help us navigate, find distance, etc.
  92   const static int INEXISTENT_VERTEX=-3;
  93
  94   bool _verbose;
  95
  96   bool _ignore_nearest_is_mirror;
  97
  98   /// Picture of how things will work... Copy 0--pi part of the 0--2pi
  99   /// cylinder into a region 2pi--2pi+ a bit of a Euclidean plane. Below we
 100   /// show points labelled by + that have a mirror image in this
 101   /// manner, while points labelled by * do not have a mirror image.
 102   ///
 103   ///      |           .     |
 104   ///      |           .     |
 105   ///      |           .     |
 106   ///      |           .     |
 107   ///      |        2  .     |
 108   ///      |        *  .     |
 109   ///      | +         . +   |
 110   ///      | 0         . 1   |
 111   ///      |           .     |
 112   ///      0          2pi   2pi + a bit
 113   ///
 114   /// Each "true" point has its true "cylinder" index (the index that
 115   /// is known externally to this class) as well as euclidean plane
 116   /// indices (main_index and mirror index in the MirrorVertexInfo
 117   /// structure), which are private concepts of this class.
 118   ///
 119   /// In above picture our structures would hold the following info
 120   /// (the picture shows the euclidean-plane numbering)
 121   ///
 122   /// _mirror_info[cylinder_index = 0] = (0, 1)
 123   /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX)
 124   ///
 125   /// We also need to be able to go from the euclidean plane indices
 126   /// back to the "true" cylinder index, and for this purpose we use
 127   /// the std::vector _cylinder_index_of_plane_vertex[...], which in the above example has
 128   /// the following contents
 129   ///
 130   /// _cylinder_index_of_plane_vertex[0] = 0
 131   /// _cylinder_index_of_plane_vertex[1] = 0
 132   /// _cylinder_index_of_plane_vertex[2] = 1
 133   ///
 134
 135   ///
 136   struct MirrorVertexInfo {
 137     /// index of the given point (appearing in the range 0--2pi) in the
 138     /// 0--2pi euclidean plane structure (position will coincide with
 139     /// that on the 0--2pi cylinder, but index labelling it will be
 140     /// different)
 141     int main_index;
 142     /// index of the mirror point (appearing in the range 2pi--3pi) in the
 143     /// 0--3pi euclidean plane structure
 144     int mirror_index;
 145   };
 146
 147   // for each "true" vertex we have reference to indices in the euclidean
 148   // plane structure
 149   std::vector<MirrorVertexInfo> _mirror_info;
 150   // for each index in the euclidean 0--2pi plane structure we want to
 151   // be able to get back to the "true" vertex index on the overall
 152   // 0--2pi cylinder structure
 153   std::vector<int> _cylinder_index_of_plane_vertex;
 154
 155   // NB: we define POINTERS here because the initialisation gave
 156   //     us problems (things crashed!), perhaps because in practice
 157   //     we were making a copy without being careful and defining
 158   //     a proper copy constructor.
 159   DnnPlane * _DNN;
 160
 161   /// given a phi value in the 0--pi range return one
 162   /// in the 2pi--3pi range; whereas if it is in the pi-2pi range then
 163   /// remap it to be inthe range (-pi)--0.
 164   inline EtaPhi _remap_phi(const EtaPhi & point) {
 165     double phi = point.second;
 166     if (phi < pi) { phi += twopi ;} else {phi -= twopi;}
 167     return EtaPhi(point.first, phi);}
 168
 169
 170   //----------------------------------------------------------------------
 171   /// Actions here are similar to those in the
 172   /// Dnn3piCylinder::_RegisterCylinderPoint case, however here we do
 173   /// NOT create the mirror point -- instead we initialise the structure
 174   /// as if there were no need for the mirror point.
 175   ///
 176   /// ADDITIONALLY push the cylinder_point onto the vector plane_points.
 177   void _RegisterCylinderPoint (const EtaPhi & cylinder_point,
 178                                std::vector<EtaPhi> & plane_points);
 179
 180   /// For each plane point specified in the vector plane_indices,
 181   /// establish whether there is a need to create a mirror point
 182   /// according to the following criteria:
 183   ///
 184   /// . phi < pi
 185   /// . mirror does not already exist
 186   /// . phi < NearestNeighbourDistance
 187   ///   (if this is not true then there is no way that its mirror point
 188   ///   could have a nearer neighbour).
 189   ///
 190   /// If conditions all hold, then create the mirror point, insert it
 191   /// into the _DNN structure, adjusting any nearest neighbours, and
 192   /// return the list of plane points whose nearest neighbours have
 193   /// changed (this will include the new neighbours that have just been
 194   /// added)
 195   void _CreateNecessaryMirrorPoints(
 196                           const std::vector<int> & plane_indices,
 197                           std::vector<int> & updated_plane_points);
 198
 199 };
 200
 201
 202 // here follow some inline implementations of the simpler of the
 203 // functions defined above
 204
 205 //----------------------------------------------------------------------
 206 /// Note: one of the difficulties of the 0--2pi mapping is that
 207 /// a point may have its mirror copy as its own nearest neighbour
 208 /// (if no other point is within a distance of 2pi). This does
 209 /// not matter for the kt_algorithm with
 210 /// reasonable values of radius, but might matter for a general use
 211 /// of this algorithm -- depending on whether or not the user has
 212 /// initialised the class with instructions to ignore this problem the
 213 /// program will detect and ignore it, or crash.
 214 inline int Dnn2piCylinder::NearestNeighbourIndex(const int & current) const {
 215   int main_index = _mirror_info[current].main_index;
 216   int mirror_index = _mirror_info[current].mirror_index;
 217   int plane_index;
 218   if (mirror_index == INEXISTENT_VERTEX ) {
 219     plane_index = _DNN->NearestNeighbourIndex(main_index);
 220   } else {
 221     plane_index = (
 222         _DNN->NearestNeighbourDistance(main_index) <
 223         _DNN->NearestNeighbourDistance(mirror_index)) ?
 224       _DNN->NearestNeighbourIndex(main_index) :
 225       _DNN->NearestNeighbourIndex(mirror_index) ;
 226   }
 227   int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index];
 228   // either the user has acknowledged the fact that they may get the
 229   // mirror copy as the closest point, or crash if it should occur
 230   // that mirror copy is the closest point.
 231   assert(_ignore_nearest_is_mirror || this_cylinder_index != current);
 232   //if (this_cylinder_index == current) {
 233   //  cerr << "WARNING point "<<current<<
 234   //    " has its mirror copy as its own nearest neighbour"<<endl;
 235   //}
 236   return this_cylinder_index;
 237 }
 238
 239 inline double Dnn2piCylinder::NearestNeighbourDistance(const int & current) const {
 240   int main_index = _mirror_info[current].main_index;
 241   int mirror_index = _mirror_info[current].mirror_index;
 242   if (mirror_index == INEXISTENT_VERTEX ) {
 243     return _DNN->NearestNeighbourDistance(main_index);
 244   } else {
 245     return (
 246         _DNN->NearestNeighbourDistance(main_index) <
 247         _DNN->NearestNeighbourDistance(mirror_index)) ?
 248       _DNN->NearestNeighbourDistance(main_index) :
 249       _DNN->NearestNeighbourDistance(mirror_index) ;
 250   }
 251
 252 }
 253
 254 inline bool Dnn2piCylinder::Valid(const int & index) const {
 255   return (_DNN->Valid(_mirror_info[index].main_index));
 256 }
 257
 258
 259 inline Dnn2piCylinder::~Dnn2piCylinder() {
 260   delete _DNN;
 261 }
 262
 263
 264 } // fastjet namespace
 265
 266 #endif //  __FASTJET_DNN2PICYLINDER_HH__
 267 #endif //DROP_CGAL