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 | |
---|
38 | FASTJET_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 |
---|
47 | class 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. |
---|
206 | inline 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 | |
---|
231 | inline 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 | |
---|
246 | inline bool Dnn3piCylinder::Valid(const int & index) const { |
---|
247 | return (_DNN->Valid(_mirror_info[index].main_index)); |
---|
248 | } |
---|
249 | |
---|
250 | |
---|
251 | inline Dnn3piCylinder::~Dnn3piCylinder() { |
---|
252 | delete _DNN; |
---|
253 | } |
---|
254 | |
---|
255 | |
---|
256 | FASTJET_END_NAMESPACE |
---|
257 | |
---|
258 | #endif // __FASTJET_DNN3PICYLINDER_HH__ |
---|
259 | #endif // DROP_CGAL |
---|