FastJet 3.0.4
Dnn2piCylinder.hh
00001 //STARTHEADER
00002 // $Id: Dnn2piCylinder.hh 2577 2011-09-13 15:11:38Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00031 #ifndef __FASTJET_DNN2PICYLINDER_HH__
00032 #define __FASTJET_DNN2PICYLINDER_HH__
00033 
00034 #include "fastjet/internal/DynamicNearestNeighbours.hh"
00035 #include "fastjet/internal/DnnPlane.hh"
00036 #include "fastjet/internal/numconsts.hh"
00037 
00038 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00039 
00040 
00041 /// \if internal_doc
00042 /// @ingroup internal
00043 /// \class Dnn2piCylinder
00044 /// class derived from DynamicNearestNeighbours that provides an
00045 /// implementation for the surface of cylinder (using one 
00046 /// DnnPlane object spanning 0--2pi).
00047 /// \endif
00048 class Dnn2piCylinder : public DynamicNearestNeighbours {
00049  public:
00050   /// empty initaliser
00051   Dnn2piCylinder() {}
00052 
00053   /// Initialiser from a set of points on an Eta-Phi plane, where
00054   /// eta can have an arbitrary ranges and phi must be in range
00055   /// 0 <= phi < 2pi;
00056   /// 
00057   /// NB: this class is more efficient than the plain Dnn4piCylinder
00058   /// class, but can give wrong answers when the nearest neighbour is
00059   /// further away than 2pi (in this case a point's nearest neighbour
00060   /// becomes itself, because it is considered to be a distance 2pi
00061   /// away). For the kt-algorithm (e.g.) this is actually not a
00062   /// problem (the distance need only be accurate when it is less than
00063   /// R, assuming R<2pi [not necessarily always the case as of
00064   /// 2010-11-19, when we've removed the requirement R<pi/2 in the
00065   /// JetDefinition constructor]), so we can tell the routine to
00066   /// ignore this problem -- alternatively the routine will crash if
00067   /// it detects it occurring (only when finding the nearest neighbour
00068   /// index, not its distance).
00069   Dnn2piCylinder(const std::vector<EtaPhi> &,
00070                  const bool & ignore_nearest_is_mirror = false,
00071                  const bool & verbose = false );
00072 
00073   /// Returns the index of  the nearest neighbour of point labelled
00074   /// by ii (assumes ii is valid)
00075   int NearestNeighbourIndex(const int & ii) const ;
00076 
00077   /// Returns the distance to the nearest neighbour of point labelled
00078   /// by index ii (assumes ii is valid)
00079   double NearestNeighbourDistance(const int & ii) const ;
00080 
00081   /// Returns true iff the given index corresponds to a point that
00082   /// exists in the DNN structure (meaning that it has been added, and
00083   /// not removed in the meantime)
00084   bool Valid(const int & index) const;
00085 
00086   void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
00087                           const std::vector<EtaPhi> & points_to_add,
00088                           std::vector<int> & indices_added,
00089                           std::vector<int> & indices_of_updated_neighbours);
00090 
00091   ~Dnn2piCylinder();
00092 
00093  private:
00094 
00095   // our extras to help us navigate, find distance, etc.
00096   const static int INEXISTENT_VERTEX=-3;
00097 
00098   bool _verbose;
00099 
00100   bool _ignore_nearest_is_mirror;
00101 
00102   /// Picture of how things will work... Copy 0--pi part of the 0--2pi
00103   /// cylinder into a region 2pi--2pi+ a bit of a Euclidean plane. Below we
00104   /// show points labelled by + that have a mirror image in this
00105   /// manner, while points labelled by * do not have a mirror image.
00106   ///                                     
00107   ///      |           .     |            
00108   ///      |           .     |            
00109   ///      |           .     |            
00110   ///      |           .     |            
00111   ///      |        2  .     |            
00112   ///      |        *  .     |            
00113   ///      | +         . +   |            
00114   ///      | 0         . 1   |
00115   ///      |           .     |
00116   ///      0          2pi   2pi + a bit
00117   ///                
00118   /// Each "true" point has its true "cylinder" index (the index that
00119   /// is known externally to this class) as well as euclidean plane
00120   /// indices (main_index and mirror index in the MirrorVertexInfo
00121   /// structure), which are private concepts of this class.
00122   /// 
00123   /// In above picture our structures would hold the following info
00124   /// (the picture shows the euclidean-plane numbering)
00125   ///
00126   /// _mirror_info[cylinder_index = 0] = (0, 1)
00127   /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX)
00128   ///
00129   /// We also need to be able to go from the euclidean plane indices
00130   /// back to the "true" cylinder index, and for this purpose we use
00131   /// the std::vector _cylinder_index_of_plane_vertex[...], which in the above example has
00132   /// the following contents
00133   ///
00134   /// _cylinder_index_of_plane_vertex[0] = 0
00135   /// _cylinder_index_of_plane_vertex[1] = 0
00136   /// _cylinder_index_of_plane_vertex[2] = 1
00137   ///
00138 
00139   /// 
00140   struct MirrorVertexInfo {
00141     /// index of the given point (appearing in the range 0--2pi) in the 
00142     /// 0--2pi euclidean plane structure (position will coincide with
00143     /// that on the 0--2pi cylinder, but index labelling it will be
00144     /// different)
00145     int main_index; 
00146     /// index of the mirror point (appearing in the range 2pi--3pi) in the
00147     /// 0--3pi euclidean plane structure
00148     int mirror_index; 
00149   };
00150 
00151   // for each "true" vertex we have reference to indices in the euclidean
00152   // plane structure
00153   std::vector<MirrorVertexInfo> _mirror_info;
00154   // for each index in the euclidean 0--2pi plane structure we want to
00155   // be able to get back to the "true" vertex index on the overall
00156   // 0--2pi cylinder structure
00157   std::vector<int> _cylinder_index_of_plane_vertex;
00158 
00159   // NB: we define POINTERS here because the initialisation gave
00160   //     us problems (things crashed!), perhaps because in practice
00161   //     we were making a copy without being careful and defining
00162   //     a proper copy constructor.
00163   DnnPlane * _DNN;
00164 
00165   /// given a phi value in the 0--pi range return one 
00166   /// in the 2pi--3pi range; whereas if it is in the pi-2pi range then
00167   /// remap it to be inthe range (-pi)--0.
00168   inline EtaPhi _remap_phi(const EtaPhi & point) {
00169     double phi = point.second;
00170     if (phi < pi) { phi += twopi ;} else {phi -= twopi;}
00171     return EtaPhi(point.first, phi);}
00172 
00173 
00174   //----------------------------------------------------------------------
00175   /// Actions here are similar to those in the
00176   /// Dnn3piCylinder::_RegisterCylinderPoint case, however here we do
00177   /// NOT create the mirror point -- instead we initialise the structure
00178   /// as if there were no need for the mirror point.
00179   ///
00180   /// ADDITIONALLY push the cylinder_point onto the vector plane_points.
00181   void _RegisterCylinderPoint (const EtaPhi & cylinder_point,
00182                                std::vector<EtaPhi> & plane_points);
00183 
00184   /// For each plane point specified in the vector plane_indices,
00185   /// establish whether there is a need to create a mirror point
00186   /// according to the following criteria:
00187   ///
00188   /// . phi < pi
00189   /// . mirror does not already exist
00190   /// . phi < NearestNeighbourDistance 
00191   ///   (if this is not true then there is no way that its mirror point
00192   ///   could have a nearer neighbour).
00193   ///
00194   /// If conditions all hold, then create the mirror point, insert it
00195   /// into the _DNN structure, adjusting any nearest neighbours, and
00196   /// return the list of plane points whose nearest neighbours have
00197   /// changed (this will include the new neighbours that have just been
00198   /// added)
00199   void _CreateNecessaryMirrorPoints(
00200                           const std::vector<int> & plane_indices,
00201                           std::vector<int> & updated_plane_points);
00202 
00203 };
00204 
00205 
00206 // here follow some inline implementations of the simpler of the
00207 // functions defined above
00208 
00209 //----------------------------------------------------------------------
00210 /// Note: one of the difficulties of the 0--2pi mapping is that
00211 /// a point may have its mirror copy as its own nearest neighbour
00212 /// (if no other point is within a distance of 2pi). This does
00213 /// not matter for the kt_algorithm with
00214 /// reasonable values of radius, but might matter for a general use
00215 /// of this algorithm -- depending on whether or not the user has
00216 /// initialised the class with instructions to ignore this problem the
00217 /// program will detect and ignore it, or crash.
00218 inline int Dnn2piCylinder::NearestNeighbourIndex(const int & current) const {
00219   int main_index = _mirror_info[current].main_index;
00220   int mirror_index = _mirror_info[current].mirror_index;
00221   int plane_index;
00222   if (mirror_index == INEXISTENT_VERTEX ) {
00223     plane_index = _DNN->NearestNeighbourIndex(main_index);
00224   } else {
00225     plane_index = (
00226         _DNN->NearestNeighbourDistance(main_index) < 
00227         _DNN->NearestNeighbourDistance(mirror_index)) ? 
00228       _DNN->NearestNeighbourIndex(main_index) : 
00229       _DNN->NearestNeighbourIndex(mirror_index) ; 
00230   }
00231   int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index];
00232   // either the user has acknowledged the fact that they may get the
00233   // mirror copy as the closest point, or crash if it should occur
00234   // that mirror copy is the closest point.
00235   assert(_ignore_nearest_is_mirror || this_cylinder_index != current);
00236   //if (this_cylinder_index == current) {
00237   //  cerr << "WARNING point "<<current<<
00238   //    " has its mirror copy as its own nearest neighbour"<<endl;
00239   //}
00240   return this_cylinder_index;
00241 }
00242 
00243 inline double Dnn2piCylinder::NearestNeighbourDistance(const int & current) const {
00244   int main_index = _mirror_info[current].main_index;
00245   int mirror_index = _mirror_info[current].mirror_index;
00246   if (mirror_index == INEXISTENT_VERTEX ) {
00247     return _DNN->NearestNeighbourDistance(main_index);
00248   } else {
00249     return (
00250         _DNN->NearestNeighbourDistance(main_index) < 
00251         _DNN->NearestNeighbourDistance(mirror_index)) ? 
00252       _DNN->NearestNeighbourDistance(main_index) : 
00253       _DNN->NearestNeighbourDistance(mirror_index) ; 
00254   }
00255  
00256 }
00257 
00258 inline bool Dnn2piCylinder::Valid(const int & index) const {
00259   return (_DNN->Valid(_mirror_info[index].main_index));
00260 }
00261 
00262 
00263 inline Dnn2piCylinder::~Dnn2piCylinder() {
00264   delete _DNN; 
00265 }
00266 
00267 
00268 FASTJET_END_NAMESPACE
00269 
00270 #endif //  __FASTJET_DNN2PICYLINDER_HH__
00271 #endif //DROP_CGAL 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends