Fork me on GitHub

source: svn/trunk/external/fastjet/internal/Dnn2piCylinder.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: 10.5 KB
Line 
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
38FASTJET_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
48class 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.
218inline 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
243inline 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
258inline bool Dnn2piCylinder::Valid(const int & index) const {
259 return (_DNN->Valid(_mirror_info[index].main_index));
260}
261
262
263inline Dnn2piCylinder::~Dnn2piCylinder() {
264 delete _DNN;
265}
266
267
268FASTJET_END_NAMESPACE
269
270#endif // __FASTJET_DNN2PICYLINDER_HH__
271#endif //DROP_CGAL
Note: See TracBrowser for help on using the repository browser.