1 | //STARTHEADER |
---|
2 | // $Id: Dnn3piCylinder.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/Dnn3piCylinder.hh" |
---|
33 | using namespace std; |
---|
34 | |
---|
35 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
36 | |
---|
37 | //---------------------------------------------------------------------- |
---|
38 | /// initialiser... |
---|
39 | Dnn3piCylinder::Dnn3piCylinder( |
---|
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 | //plane_points.reserve(2*input_points.size()); |
---|
48 | |
---|
49 | for (unsigned int i=0; i < input_points.size(); i++) { |
---|
50 | _RegisterCylinderPoint(input_points[i], plane_points); |
---|
51 | } |
---|
52 | |
---|
53 | if (_verbose) cout << "============== Preparing _DNN" << endl; |
---|
54 | _DNN = new DnnPlane(plane_points, verbose); |
---|
55 | } |
---|
56 | |
---|
57 | |
---|
58 | //---------------------------------------------------------------------- |
---|
59 | /// What on earth does this do? |
---|
60 | /// |
---|
61 | /// Example: last true "cylinder" index was 15 |
---|
62 | /// last plane index was 23 |
---|
63 | /// |
---|
64 | /// Then: _cylinder_index_of_plane_vertex.size() = 24 and |
---|
65 | /// _mirror_info.size() = 16 |
---|
66 | /// |
---|
67 | /// IF cylinder_point's phi < pi then |
---|
68 | /// create: _mirror_info[16] = (main_index = 24, mirror_index=25) |
---|
69 | /// _cylinder_index_of_plane_vertex[24] = 16 |
---|
70 | /// _cylinder_index_of_plane_vertex[25] = 16 |
---|
71 | /// ELSE |
---|
72 | /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..) |
---|
73 | /// _cylinder_index_of_plane_vertex[24] = 16 |
---|
74 | /// |
---|
75 | /// ADDITIONALLY push the cylinder_point (and if it exists the mirror |
---|
76 | /// copy) onto the vector plane_points. |
---|
77 | void Dnn3piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point, |
---|
78 | vector<EtaPhi> & plane_points) { |
---|
79 | double phi = cylinder_point.second; |
---|
80 | assert(phi >= 0.0 && phi < 2*pi); |
---|
81 | |
---|
82 | // do main point |
---|
83 | MirrorVertexInfo mvi; |
---|
84 | mvi.main_index = _cylinder_index_of_plane_vertex.size(); |
---|
85 | _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); |
---|
86 | plane_points.push_back(cylinder_point); |
---|
87 | |
---|
88 | // do mirror point if need be |
---|
89 | if (phi < pi) { |
---|
90 | mvi.mirror_index = _cylinder_index_of_plane_vertex.size(); |
---|
91 | _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); |
---|
92 | plane_points.push_back(_remap_phi(cylinder_point)); |
---|
93 | } else { |
---|
94 | mvi.mirror_index = INEXISTENT_VERTEX; |
---|
95 | } |
---|
96 | |
---|
97 | // |
---|
98 | _mirror_info.push_back(mvi); |
---|
99 | } |
---|
100 | |
---|
101 | |
---|
102 | //---------------------------------------------------------------------- |
---|
103 | /// insertion and removal of points |
---|
104 | void Dnn3piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove, |
---|
105 | const vector<EtaPhi> & points_to_add, |
---|
106 | vector<int> & indices_added, |
---|
107 | vector<int> & indices_of_updated_neighbours) { |
---|
108 | |
---|
109 | // translate from "cylinder" indices of points to remove to the |
---|
110 | // plane indices of points to remove, bearing in mind that sometimes |
---|
111 | // there are multple plane points to remove. |
---|
112 | vector<int> plane_indices_to_remove; |
---|
113 | for (unsigned int i=0; i < indices_to_remove.size(); i++) { |
---|
114 | MirrorVertexInfo * mvi; |
---|
115 | mvi = & _mirror_info[indices_to_remove[i]]; |
---|
116 | plane_indices_to_remove.push_back(mvi->main_index); |
---|
117 | if (mvi->mirror_index != INEXISTENT_VERTEX) { |
---|
118 | plane_indices_to_remove.push_back(mvi->mirror_index); |
---|
119 | } |
---|
120 | } |
---|
121 | |
---|
122 | // given "cylinder" points to add get hold of the list of |
---|
123 | // plane-points to add. |
---|
124 | vector<EtaPhi> plane_points_to_add; |
---|
125 | indices_added.clear(); |
---|
126 | for (unsigned int i=0; i < points_to_add.size(); i++) { |
---|
127 | indices_added.push_back(_mirror_info.size()); |
---|
128 | _RegisterCylinderPoint(points_to_add[i], plane_points_to_add); |
---|
129 | } |
---|
130 | |
---|
131 | // now get the hard work done (note that we need to supply the |
---|
132 | // plane_indices_added vector but that we will not actually check |
---|
133 | // its contents in any way -- the indices_added that is actually |
---|
134 | // returned has been calculated above). |
---|
135 | vector<int> updated_plane_neighbours, plane_indices_added; |
---|
136 | _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add, |
---|
137 | plane_indices_added, updated_plane_neighbours); |
---|
138 | |
---|
139 | // extract, from the updated_plane_neighbours, the set of cylinder |
---|
140 | // neighbours that have changed |
---|
141 | set<int> index_set; |
---|
142 | unsigned int i; |
---|
143 | for (i=0; i < updated_plane_neighbours.size(); i++) { |
---|
144 | index_set.insert( |
---|
145 | _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);} |
---|
146 | |
---|
147 | // decant the set into the vector that needs to be returned |
---|
148 | indices_of_updated_neighbours.clear(); |
---|
149 | for (set<int>::iterator iter = index_set.begin(); |
---|
150 | iter != index_set.end(); iter++) { |
---|
151 | indices_of_updated_neighbours.push_back(*iter); |
---|
152 | } |
---|
153 | } |
---|
154 | |
---|
155 | |
---|
156 | FASTJET_END_NAMESPACE |
---|
157 | |
---|
158 | #endif // DROP_CGAL |
---|
159 | |
---|