Fork me on GitHub

source: svn/trunk/external/fastjet/internal/Dnn3piCylinder.hh@ 1204

Last change on this file since 1204 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3

File size: 9.8 KB
Line 
1//STARTHEADER
2// $Id: Dnn3piCylinder.hh 2577 2011-09-13 15:11:38Z salam $
3//
4// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
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, see <http://www.gnu.org/licenses/>.
26//----------------------------------------------------------------------
27//ENDHEADER
28
29
30#ifndef DROP_CGAL // in case we do not have the code for CGAL
31#ifndef __FASTJET_DNN3PICYLINDER_HH__
32#define __FASTJET_DNN3PICYLINDER_HH__
33
34#include "fastjet/internal/DynamicNearestNeighbours.hh"
35#include "fastjet/internal/DnnPlane.hh"
36#include "fastjet/internal/numconsts.hh"
37
38FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39
40/// \if internal_doc
41/// @ingroup internal
42/// \class Dnn3piCylinder
43/// class derived from DynamicNearestNeighbours that provides an
44/// implementation for the surface of cylinder (using one
45/// DnnPlane object spanning 0--3pi).
46/// \endif
47class Dnn3piCylinder : public DynamicNearestNeighbours {
48 public:
49 /// empty initaliser
50 Dnn3piCylinder() {}
51
52 /// Initialiser from a set of points on an Eta-Phi plane, where
53 /// eta can have an arbitrary ranges and phi must be in range
54 /// 0 <= phi < 2pi;
55 ///
56 /// NB: this class is more efficient than the plain Dnn4piCylinder
57 /// class, but can give wrong answers when the nearest neighbour is
58 /// further away than 2pi (in this case a point's nearest neighbour
59 /// becomes itself, because it is considered to be a distance 2pi
60 /// away). For the kt-algorithm (e.g.) this is actually not a
61 /// problem (the distance need only be accurate when it is less than
62 /// R), so we can tell the routine to ignore this problem --
63 /// alternatively the routine will crash if it detects it occurring
64 /// (only when finding the nearest neighbour index, not its
65 /// distance).
66 Dnn3piCylinder(const std::vector<EtaPhi> &,
67 const bool & ignore_nearest_is_mirror = false,
68 const bool & verbose = false );
69
70 /// Returns the index of the nearest neighbour of point labelled
71 /// by ii (assumes ii is valid)
72 int NearestNeighbourIndex(const int & ii) const ;
73
74 /// Returns the distance to the nearest neighbour of point labelled
75 /// by index ii (assumes ii is valid)
76 double NearestNeighbourDistance(const int & ii) const ;
77
78 /// Returns true iff the given index corresponds to a point that
79 /// exists in the DNN structure (meaning that it has been added, and
80 /// not removed in the meantime)
81 bool Valid(const int & index) const;
82
83 void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
84 const std::vector<EtaPhi> & points_to_add,
85 std::vector<int> & indices_added,
86 std::vector<int> & indices_of_updated_neighbours);
87
88 ~Dnn3piCylinder();
89
90 private:
91
92 // our extras to help us navigate, find distance, etc.
93 const static int INEXISTENT_VERTEX=-3;
94
95 bool _verbose;
96
97 bool _ignore_nearest_is_mirror;
98
99 /// Picture of how things will work... Copy 0--pi part of the 0--2pi
100 /// cylinder into a region 2pi--3pi of a Euclidean plane. Below we
101 /// show points labelled by + that have a mirror image in this
102 /// manner, while points labelled by * do not have a mirror image.
103 ///
104 /// | . |
105 /// | . |
106 /// | . |
107 /// | . |
108 /// | 2 . |
109 /// | * . |
110 /// | + . + |
111 /// | 0 . 1 |
112 /// | . |
113 /// 0 2pi 3pi
114 ///
115 /// Each "true" point has its true "cylinder" index (the index that
116 /// is known externally to this class) as well as euclidean plane
117 /// indices (main_index and mirror index in the MirrorVertexInfo
118 /// structure), which are private concepts of this class.
119 ///
120 /// In above picture our structures would hold the following info
121 /// (the picture shows the euclidean-plane numbering)
122 ///
123 /// _mirror_info[cylinder_index = 0] = (0, 1)
124 /// _mirror_info[cylinder_index = 1] = (2, INEXISTENT_VERTEX)
125 ///
126 /// We also need to be able to go from the euclidean plane indices
127 /// back to the "true" cylinder index, and for this purpose we use
128 /// the vector _cylinder_index_of_plane_vertex[...], which in the above example has
129 /// the following contents
130 ///
131 /// _cylinder_index_of_plane_vertex[0] = 0
132 /// _cylinder_index_of_plane_vertex[1] = 0
133 /// _cylinder_index_of_plane_vertex[2] = 1
134 ///
135
136 ///
137 struct MirrorVertexInfo {
138 /// index of the given point (appearing in the range 0--2pi) in the
139 /// 0--3pi euclidean plane structure (position will coincide with
140 /// that on the 0--2pi cylinder, but index labelling it will be
141 /// different)
142 int main_index;
143 /// index of the mirror point (appearing in the range 2pi--3pi) in the
144 /// 0--3pi euclidean plane structure
145 int mirror_index;
146 };
147
148 // for each "true" vertex we have reference to indices in the euclidean
149 // plane structure
150 std::vector<MirrorVertexInfo> _mirror_info;
151 // for each index in the euclidean 0--3pi plane structure we want to
152 // be able to get back to the "true" vertex index on the overall
153 // 0--2pi cylinder structure
154 std::vector<int> _cylinder_index_of_plane_vertex;
155
156 // NB: we define POINTERS here because the initialisation gave
157 // us problems (things crashed!), perhaps because in practice
158 // we were making a copy without being careful and defining
159 // a proper copy constructor.
160 DnnPlane * _DNN;
161
162 /// given a phi value in the 0--2pi range return one
163 /// in the pi--3pi range.
164 inline EtaPhi _remap_phi(const EtaPhi & point) {
165 double phi = point.second;
166 if (phi < pi) { phi += twopi ;}
167 return EtaPhi(point.first, phi);}
168
169
170 //----------------------------------------------------------------------
171 /// What on earth does this do?
172 ///
173 /// Example: last true "cylinder" index was 15
174 /// last plane index was 23
175 ///
176 /// Then: _cylinder_index_of_plane_vertex.size() = 24 and
177 /// _mirror_info.size() = 16
178 ///
179 /// IF cylinder_point's phi < pi then
180 /// create: _mirror_info[16] = (main_index = 24, mirror_index=25)
181 /// _cylinder_index_of_plane_vertex[24] = 16
182 /// _cylinder_index_of_plane_vertex[25] = 16
183 /// ELSE
184 /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..)
185 /// _cylinder_index_of_plane_vertex[24] = 16
186 ///
187 /// ADDITIONALLY push the cylinder_point (and if it exists the mirror
188 /// copy) onto the vector plane_points.
189 void _RegisterCylinderPoint (const EtaPhi & cylinder_point,
190 std::vector<EtaPhi> & plane_points);
191};
192
193
194// here follow some inline implementations of the simpler of the
195// functions defined above
196
197//----------------------------------------------------------------------
198/// Note: one of the difficulties of the 0--3pi mapping is that
199/// a point may have its mirror copy as its own nearest neighbour
200/// (if no other point is within a distance of 2pi). This does
201/// not matter for the kt_algorithm with
202/// reasonable values of radius, but might matter for a general use
203/// of this algorithm -- depending on whether or not the user has
204/// initialised the class with instructions to ignore this problem the
205/// program will detect and ignore it, or crash.
206inline int Dnn3piCylinder::NearestNeighbourIndex(const int & current) const {
207 int main_index = _mirror_info[current].main_index;
208 int mirror_index = _mirror_info[current].mirror_index;
209 int plane_index;
210 if (mirror_index == INEXISTENT_VERTEX ) {
211 plane_index = _DNN->NearestNeighbourIndex(main_index);
212 } else {
213 plane_index = (
214 _DNN->NearestNeighbourDistance(main_index) <
215 _DNN->NearestNeighbourDistance(mirror_index)) ?
216 _DNN->NearestNeighbourIndex(main_index) :
217 _DNN->NearestNeighbourIndex(mirror_index) ;
218 }
219 int this_cylinder_index = _cylinder_index_of_plane_vertex[plane_index];
220 // either the user has acknowledged the fact that they may get the
221 // mirror copy as the closest point, or crash if it should occur
222 // that mirror copy is the closest point.
223 assert(_ignore_nearest_is_mirror || this_cylinder_index != current);
224 //if (this_cylinder_index == current) {
225 // std::cerr << "WARNING point "<<current<<
226 // " has its mirror copy as its own nearest neighbour"<<endl;
227 //}
228 return this_cylinder_index;
229}
230
231inline double Dnn3piCylinder::NearestNeighbourDistance(const int & current) const {
232 int main_index = _mirror_info[current].main_index;
233 int mirror_index = _mirror_info[current].mirror_index;
234 if (mirror_index == INEXISTENT_VERTEX ) {
235 return _DNN->NearestNeighbourDistance(main_index);
236 } else {
237 return (
238 _DNN->NearestNeighbourDistance(main_index) <
239 _DNN->NearestNeighbourDistance(mirror_index)) ?
240 _DNN->NearestNeighbourDistance(main_index) :
241 _DNN->NearestNeighbourDistance(mirror_index) ;
242 }
243
244}
245
246inline bool Dnn3piCylinder::Valid(const int & index) const {
247 return (_DNN->Valid(_mirror_info[index].main_index));
248}
249
250
251inline Dnn3piCylinder::~Dnn3piCylinder() {
252 delete _DNN;
253}
254
255
256FASTJET_END_NAMESPACE
257
258#endif // __FASTJET_DNN3PICYLINDER_HH__
259#endif // DROP_CGAL
Note: See TracBrowser for help on using the repository browser.