1 | //STARTHEADER |
---|
2 | // $Id: Dnn2piCylinder.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
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 | #include <set> |
---|
32 | #include "fastjet/internal/Dnn2piCylinder.hh" |
---|
33 | using namespace std; |
---|
34 | |
---|
35 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
36 | |
---|
37 | //---------------------------------------------------------------------- |
---|
38 | /// initialiser... |
---|
39 | Dnn2piCylinder::Dnn2piCylinder( |
---|
40 | const vector<EtaPhi> & input_points, |
---|
41 | const bool & ignore_nearest_is_mirror, |
---|
42 | const bool & verbose) { |
---|
43 | |
---|
44 | _verbose = verbose; |
---|
45 | _ignore_nearest_is_mirror = ignore_nearest_is_mirror; |
---|
46 | vector<EtaPhi> plane_points; |
---|
47 | vector<int> plane_point_indices(input_points.size()); |
---|
48 | //plane_points.reserve(2*input_points.size()); |
---|
49 | |
---|
50 | for (unsigned int i=0; i < input_points.size(); i++) { |
---|
51 | _RegisterCylinderPoint(input_points[i], plane_points); |
---|
52 | plane_point_indices[i] = i; |
---|
53 | } |
---|
54 | |
---|
55 | if (_verbose) cout << "============== Preparing _DNN" << endl; |
---|
56 | _DNN = new DnnPlane(plane_points, verbose); |
---|
57 | |
---|
58 | |
---|
59 | vector<int> updated_point_indices; // we'll not use information from this |
---|
60 | _CreateNecessaryMirrorPoints(plane_point_indices,updated_point_indices); |
---|
61 | } |
---|
62 | |
---|
63 | |
---|
64 | //---------------------------------------------------------------------- |
---|
65 | /// Actions here are similar to those in the |
---|
66 | /// Dnn3piCylinder::_RegisterCylinderPoint case, however here we do |
---|
67 | /// NOT create the mirror point -- instead we initialise the structure |
---|
68 | /// as if there were no need for the mirror point. |
---|
69 | /// |
---|
70 | /// ADDITIONALLY push the cylinder_point onto the vector plane_points. |
---|
71 | void Dnn2piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point, |
---|
72 | vector<EtaPhi> & plane_points) { |
---|
73 | double phi = cylinder_point.second; |
---|
74 | assert(phi >= 0.0 && phi < 2*pi); |
---|
75 | |
---|
76 | // do main point |
---|
77 | MirrorVertexInfo mvi; |
---|
78 | mvi.main_index = _cylinder_index_of_plane_vertex.size(); |
---|
79 | _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); |
---|
80 | plane_points.push_back(cylinder_point); |
---|
81 | mvi.mirror_index = INEXISTENT_VERTEX; |
---|
82 | |
---|
83 | // |
---|
84 | _mirror_info.push_back(mvi); |
---|
85 | } |
---|
86 | |
---|
87 | |
---|
88 | |
---|
89 | //---------------------------------------------------------------------- |
---|
90 | /// For each plane point specified in the vector plane_indices, |
---|
91 | /// establish whether there is a need to create a mirror point |
---|
92 | /// according to the following criteria: |
---|
93 | /// |
---|
94 | /// . phi < pi |
---|
95 | /// . mirror does not already exist |
---|
96 | /// . phi < NearestNeighbourDistance |
---|
97 | /// (if this is not true then there is no way that its mirror point |
---|
98 | /// could have a nearer neighbour). |
---|
99 | /// |
---|
100 | /// If conditions all hold, then create the mirror point, insert it |
---|
101 | /// into the _DNN structure, adjusting any nearest neighbours, and |
---|
102 | /// return the list of plane points whose nearest neighbours have |
---|
103 | /// changed (this will include the new neighbours that have just been |
---|
104 | /// added) |
---|
105 | void Dnn2piCylinder::_CreateNecessaryMirrorPoints( |
---|
106 | const vector<int> & plane_indices, |
---|
107 | vector<int> & updated_plane_points) { |
---|
108 | |
---|
109 | vector<EtaPhi> new_plane_points; |
---|
110 | |
---|
111 | for (size_t i = 0; i < plane_indices.size(); i++) { |
---|
112 | |
---|
113 | int ip = plane_indices[i]; // plane index |
---|
114 | EtaPhi position = _DNN->etaphi(ip); |
---|
115 | double phi = position.second; |
---|
116 | |
---|
117 | //BAD // require phi < pi |
---|
118 | //BAD if (phi >= pi) {continue;} |
---|
119 | |
---|
120 | // require absence of mirror |
---|
121 | int ic = _cylinder_index_of_plane_vertex[ip]; |
---|
122 | if (_mirror_info[ic].mirror_index != INEXISTENT_VERTEX) {continue;} |
---|
123 | |
---|
124 | //printf("%3d %3d %10.5f %10.5f %3d\n",ic, ip, phi, |
---|
125 | // _DNN->NearestNeighbourDistance(ip),_DNN->NearestNeighbourIndex(ip)); |
---|
126 | |
---|
127 | |
---|
128 | // check that we are sufficiently close to the border -- |
---|
129 | // i.e. closer than nearest neighbour distance. But RECALL: |
---|
130 | // nearest neighbourDistance actually returns a squared distance |
---|
131 | // (this was really stupid on our part -- caused considerable loss |
---|
132 | // of time ... ) |
---|
133 | double nndist = _DNN->NearestNeighbourDistance(ip); |
---|
134 | if (phi*phi >= nndist && (twopi-phi)*(twopi-phi) >= nndist) {continue;} |
---|
135 | |
---|
136 | // now proceed to prepare the point for addition |
---|
137 | new_plane_points.push_back(_remap_phi(position)); |
---|
138 | _mirror_info[ic].mirror_index = _cylinder_index_of_plane_vertex.size(); |
---|
139 | _cylinder_index_of_plane_vertex.push_back(ic); |
---|
140 | } |
---|
141 | |
---|
142 | vector<int> indices_to_remove; // empty |
---|
143 | vector<int> indices_added; // will be filled as result of call |
---|
144 | _DNN->RemoveAndAddPoints(indices_to_remove,new_plane_points,indices_added, |
---|
145 | updated_plane_points); |
---|
146 | |
---|
147 | // occasionally, addition of points might cause a change in the |
---|
148 | // nearest neighbour of a point in the 0--pi range? (But should be |
---|
149 | // impossible -- we add points beyond 2pi, so they can only be |
---|
150 | // nearest neighbours of points in the range pi--2pi; there is one |
---|
151 | // exception -- the nearest neighbour of one's self -- but in that |
---|
152 | // case this will already have been discovered, so there should be |
---|
153 | // nothing left to do). |
---|
154 | |
---|
155 | // To be on the safe side, check to see if we have updated points |
---|
156 | // with phi<pi and no current mirror point. BUT: this check, as |
---|
157 | // written, only makes sense when the mirror image was created only |
---|
158 | // beyond 2pi, which is no longer the case. Only apparent |
---|
159 | // alternative is to run separate insertions for beyond 2pi and |
---|
160 | // below phi=0, with separate checks in the two cases. But, given |
---|
161 | // that across a large number of recombinations and events in the |
---|
162 | // old method (single mirror), we never ran into this problem, it's |
---|
163 | // probably safe to assume that the arguments given above are OK! So |
---|
164 | // comment out the check... |
---|
165 | //NOTNEEDED for (size_t i = 0; i < updated_plane_points.size(); i++) { |
---|
166 | //NOTNEEDED int ip = updated_plane_points[i]; |
---|
167 | //NOTNEEDED double phi = _DNN->phi(ip); |
---|
168 | //NOTNEEDED int ic = _cylinder_index_of_plane_vertex[ip]; |
---|
169 | //NOTNEEDED assert (!(phi < pi && _mirror_info[ic].mirror_index == INEXISTENT_VERTEX)); |
---|
170 | //NOTNEEDED } |
---|
171 | // alternative recursive code |
---|
172 | //vector<int> extra_updated_plane_points; |
---|
173 | //_CreateNecessaryMirrorPoints(updated_plane_points,extra_updated_plane_points) |
---|
174 | //updated_plane_points.push_back(extra_updated_plane_points); |
---|
175 | } |
---|
176 | |
---|
177 | |
---|
178 | |
---|
179 | //---------------------------------------------------------------------- |
---|
180 | /// insertion and removal of points |
---|
181 | void Dnn2piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove, |
---|
182 | const vector<EtaPhi> & points_to_add, |
---|
183 | vector<int> & indices_added, |
---|
184 | vector<int> & indices_of_updated_neighbours) { |
---|
185 | |
---|
186 | // translate from "cylinder" indices of points to remove to the |
---|
187 | // plane indices of points to remove, bearing in mind that sometimes |
---|
188 | // there are multple plane points to remove. |
---|
189 | vector<int> plane_indices_to_remove; |
---|
190 | for (unsigned int i=0; i < indices_to_remove.size(); i++) { |
---|
191 | MirrorVertexInfo * mvi; |
---|
192 | mvi = & _mirror_info[indices_to_remove[i]]; |
---|
193 | plane_indices_to_remove.push_back(mvi->main_index); |
---|
194 | if (mvi->mirror_index != INEXISTENT_VERTEX) { |
---|
195 | plane_indices_to_remove.push_back(mvi->mirror_index); |
---|
196 | } |
---|
197 | } |
---|
198 | |
---|
199 | // given "cylinder" points to add get hold of the list of |
---|
200 | // plane-points to add. |
---|
201 | vector<EtaPhi> plane_points_to_add; |
---|
202 | indices_added.clear(); |
---|
203 | for (unsigned int i=0; i < points_to_add.size(); i++) { |
---|
204 | indices_added.push_back(_mirror_info.size()); |
---|
205 | _RegisterCylinderPoint(points_to_add[i], plane_points_to_add); |
---|
206 | } |
---|
207 | |
---|
208 | // now get the hard work done (note that we need to supply the |
---|
209 | // plane_indices_added vector but that we will not actually check |
---|
210 | // its contents in any way -- the indices_added that is actually |
---|
211 | // returned has been calculated above). |
---|
212 | vector<int> updated_plane_neighbours, plane_indices_added; |
---|
213 | _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add, |
---|
214 | plane_indices_added, updated_plane_neighbours); |
---|
215 | |
---|
216 | vector<int> extra_updated_plane_neighbours; |
---|
217 | _CreateNecessaryMirrorPoints(updated_plane_neighbours, |
---|
218 | extra_updated_plane_neighbours); |
---|
219 | |
---|
220 | // extract, from the updated_plane_neighbours, and |
---|
221 | // extra_updated_plane_neighbours, the set of cylinder neighbours |
---|
222 | // that have changed |
---|
223 | set<int> index_set; |
---|
224 | unsigned int i; |
---|
225 | for (i=0; i < updated_plane_neighbours.size(); i++) { |
---|
226 | index_set.insert( |
---|
227 | _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);} |
---|
228 | for (i=0; i < extra_updated_plane_neighbours.size(); i++) { |
---|
229 | index_set.insert( |
---|
230 | _cylinder_index_of_plane_vertex[extra_updated_plane_neighbours[i]]);} |
---|
231 | |
---|
232 | // decant the set into the vector that needs to be returned |
---|
233 | indices_of_updated_neighbours.clear(); |
---|
234 | for (set<int>::iterator iter = index_set.begin(); |
---|
235 | iter != index_set.end(); iter++) { |
---|
236 | indices_of_updated_neighbours.push_back(*iter); |
---|
237 | } |
---|
238 | } |
---|
239 | |
---|
240 | FASTJET_END_NAMESPACE |
---|
241 | |
---|
242 | #endif // DROP_CGAL |
---|