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