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