FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: Dnn3piCylinder.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_DNN3PICYLINDER_HH__ 00032 #define __FASTJET_DNN3PICYLINDER_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 /// \if internal_doc 00041 /// @ingroup internal 00042 /// \class Dnn3piCylinder 00043 /// class derived from DynamicNearestNeighbours that provides an 00044 /// implementation for the surface of cylinder (using one 00045 /// DnnPlane object spanning 0--3pi). 00046 /// \endif 00047 class Dnn3piCylinder : public DynamicNearestNeighbours { 00048 public: 00049 /// empty initaliser 00050 Dnn3piCylinder() {} 00051 00052 /// Initialiser from a set of points on an Eta-Phi plane, where 00053 /// eta can have an arbitrary ranges and phi must be in range 00054 /// 0 <= phi < 2pi; 00055 /// 00056 /// NB: this class is more efficient than the plain Dnn4piCylinder 00057 /// class, but can give wrong answers when the nearest neighbour is 00058 /// further away than 2pi (in this case a point's nearest neighbour 00059 /// becomes itself, because it is considered to be a distance 2pi 00060 /// away). For the kt-algorithm (e.g.) this is actually not a 00061 /// problem (the distance need only be accurate when it is less than 00062 /// R), so we can tell the routine to ignore this problem -- 00063 /// alternatively the routine will crash if it detects it occurring 00064 /// (only when finding the nearest neighbour index, not its 00065 /// distance). 00066 Dnn3piCylinder(const std::vector<EtaPhi> &, 00067 const bool & ignore_nearest_is_mirror = false, 00068 const bool & verbose = false ); 00069 00070 /// Returns the index of the nearest neighbour of point labelled 00071 /// by ii (assumes ii is valid) 00072 int NearestNeighbourIndex(const int & ii) const ; 00073 00074 /// Returns the distance to the nearest neighbour of point labelled 00075 /// by index ii (assumes ii is valid) 00076 double NearestNeighbourDistance(const int & ii) const ; 00077 00078 /// Returns true iff the given index corresponds to a point that 00079 /// exists in the DNN structure (meaning that it has been added, and 00080 /// not removed in the meantime) 00081 bool Valid(const int & index) const; 00082 00083 void RemoveAndAddPoints(const std::vector<int> & indices_to_remove, 00084 const std::vector<EtaPhi> & points_to_add, 00085 std::vector<int> & indices_added, 00086 std::vector<int> & indices_of_updated_neighbours); 00087 00088 ~Dnn3piCylinder(); 00089 00090 private: 00091 00092 // our extras to help us navigate, find distance, etc. 00093 const static int INEXISTENT_VERTEX=-3; 00094 00095 bool _verbose; 00096 00097 bool _ignore_nearest_is_mirror; 00098 00099 /// Picture of how things will work... Copy 0--pi part of the 0--2pi 00100 /// cylinder into a region 2pi--3pi of a Euclidean plane. Below we 00101 /// show points labelled by + that have a mirror image in this 00102 /// manner, while points labelled by * do not have a mirror image. 00103 /// 00104 /// | . | 00105 /// | . | 00106 /// | . | 00107 /// | . | 00108 /// | 2 . | 00109 /// | * . | 00110 /// | + . + | 00111 /// | 0 . 1 | 00112 /// | . | 00113 /// 0 2pi 3pi 00114 /// 00115 /// Each "true" point has its true "cylinder" index (the index that 00116 /// is known externally to this class) as well as euclidean plane 00117 /// indices (main_index and mirror index in the MirrorVertexInfo 00118 /// structure), which are private concepts of this class. 00119 /// 00120 /// In above picture our structures would hold the following info 00121 /// (the picture shows the euclidean-plane numbering) 00122 /// 00123 /// _mirror_info[cylinder_index = 0] = (0, 1) 00124 /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX) 00125 /// 00126 /// We also need to be able to go from the euclidean plane indices 00127 /// back to the "true" cylinder index, and for this purpose we use 00128 /// the vector _cylinder_index_of_plane_vertex[...], which in the above example has 00129 /// the following contents 00130 /// 00131 /// _cylinder_index_of_plane_vertex[0] = 0 00132 /// _cylinder_index_of_plane_vertex[1] = 0 00133 /// _cylinder_index_of_plane_vertex[2] = 1 00134 /// 00135 00136 /// 00137 struct MirrorVertexInfo { 00138 /// index of the given point (appearing in the range 0--2pi) in the 00139 /// 0--3pi euclidean plane structure (position will coincide with 00140 /// that on the 0--2pi cylinder, but index labelling it will be 00141 /// different) 00142 int main_index; 00143 /// index of the mirror point (appearing in the range 2pi--3pi) in the 00144 /// 0--3pi euclidean plane structure 00145 int mirror_index; 00146 }; 00147 00148 // for each "true" vertex we have reference to indices in the euclidean 00149 // plane structure 00150 std::vector<MirrorVertexInfo> _mirror_info; 00151 // for each index in the euclidean 0--3pi plane structure we want to 00152 // be able to get back to the "true" vertex index on the overall 00153 // 0--2pi cylinder structure 00154 std::vector<int> _cylinder_index_of_plane_vertex; 00155 00156 // NB: we define POINTERS here because the initialisation gave 00157 // us problems (things crashed!), perhaps because in practice 00158 // we were making a copy without being careful and defining 00159 // a proper copy constructor. 00160 DnnPlane * _DNN; 00161 00162 /// given a phi value in the 0--2pi range return one 00163 /// in the pi--3pi range. 00164 inline EtaPhi _remap_phi(const EtaPhi & point) { 00165 double phi = point.second; 00166 if (phi < pi) { phi += twopi ;} 00167 return EtaPhi(point.first, phi);} 00168 00169 00170 //---------------------------------------------------------------------- 00171 /// What on earth does this do? 00172 /// 00173 /// Example: last true "cylinder" index was 15 00174 /// last plane index was 23 00175 /// 00176 /// Then: _cylinder_index_of_plane_vertex.size() = 24 and 00177 /// _mirror_info.size() = 16 00178 /// 00179 /// IF cylinder_point's phi < pi then 00180 /// create: _mirror_info[16] = (main_index = 24, mirror_index=25) 00181 /// _cylinder_index_of_plane_vertex[24] = 16 00182 /// _cylinder_index_of_plane_vertex[25] = 16 00183 /// ELSE 00184 /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..) 00185 /// _cylinder_index_of_plane_vertex[24] = 16 00186 /// 00187 /// ADDITIONALLY push the cylinder_point (and if it exists the mirror 00188 /// copy) onto the vector plane_points. 00189 void _RegisterCylinderPoint (const EtaPhi & cylinder_point, 00190 std::vector<EtaPhi> & plane_points); 00191 }; 00192 00193 00194 // here follow some inline implementations of the simpler of the 00195 // functions defined above 00196 00197 //---------------------------------------------------------------------- 00198 /// Note: one of the difficulties of the 0--3pi mapping is that 00199 /// a point may have its mirror copy as its own nearest neighbour 00200 /// (if no other point is within a distance of 2pi). This does 00201 /// not matter for the kt_algorithm with 00202 /// reasonable values of radius, but might matter for a general use 00203 /// of this algorithm -- depending on whether or not the user has 00204 /// initialised the class with instructions to ignore this problem the 00205 /// program will detect and ignore it, or crash. 00206 inline int Dnn3piCylinder::NearestNeighbourIndex(const int & current) const { 00207 int main_index = _mirror_info[current].main_index; 00208 int mirror_index = _mirror_info[current].mirror_index; 00209 int plane_index; 00210 if (mirror_index == INEXISTENT_VERTEX ) { 00211 plane_index = _DNN->NearestNeighbourIndex(main_index); 00212 } else { 00213 plane_index = ( 00214 _DNN->NearestNeighbourDistance(main_index) < 00215 _DNN->NearestNeighbourDistance(mirror_index)) ? 00216 _DNN->NearestNeighbourIndex(main_index) : 00217 _DNN->NearestNeighbourIndex(mirror_index) ; 00218 } 00219 int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index]; 00220 // either the user has acknowledged the fact that they may get the 00221 // mirror copy as the closest point, or crash if it should occur 00222 // that mirror copy is the closest point. 00223 assert(_ignore_nearest_is_mirror || this_cylinder_index != current); 00224 //if (this_cylinder_index == current) { 00225 // std::cerr << "WARNING point "<<current<< 00226 // " has its mirror copy as its own nearest neighbour"<<endl; 00227 //} 00228 return this_cylinder_index; 00229 } 00230 00231 inline double Dnn3piCylinder::NearestNeighbourDistance(const int & current) const { 00232 int main_index = _mirror_info[current].main_index; 00233 int mirror_index = _mirror_info[current].mirror_index; 00234 if (mirror_index == INEXISTENT_VERTEX ) { 00235 return _DNN->NearestNeighbourDistance(main_index); 00236 } else { 00237 return ( 00238 _DNN->NearestNeighbourDistance(main_index) < 00239 _DNN->NearestNeighbourDistance(mirror_index)) ? 00240 _DNN->NearestNeighbourDistance(main_index) : 00241 _DNN->NearestNeighbourDistance(mirror_index) ; 00242 } 00243 00244 } 00245 00246 inline bool Dnn3piCylinder::Valid(const int & index) const { 00247 return (_DNN->Valid(_mirror_info[index].main_index)); 00248 } 00249 00250 00251 inline Dnn3piCylinder::~Dnn3piCylinder() { 00252 delete _DNN; 00253 } 00254 00255 00256 FASTJET_END_NAMESPACE 00257 00258 #endif // __FASTJET_DNN3PICYLINDER_HH__ 00259 #endif // DROP_CGAL