FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
Dnn3piCylinder.cc
1 //FJSTARTHEADER
2 // $Id: Dnn3piCylinder.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 
32 #ifndef DROP_CGAL // in case we do not have the code for CGAL
33 #include <set>
34 #include "fastjet/internal/Dnn3piCylinder.hh"
35 using namespace std;
36 
37 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38 
39 //----------------------------------------------------------------------
40 /// initialiser...
41 Dnn3piCylinder::Dnn3piCylinder(
42  const vector<EtaPhi> & input_points,
43  const bool & ignore_nearest_is_mirror,
44  const bool & verbose) {
45 
46  _verbose = verbose;
47  _ignore_nearest_is_mirror = ignore_nearest_is_mirror;
48  vector<EtaPhi> plane_points;
49  //plane_points.reserve(2*input_points.size());
50 
51  for (unsigned int i=0; i < input_points.size(); i++) {
52  _RegisterCylinderPoint(input_points[i], plane_points);
53  }
54 
55  if (_verbose) cout << "============== Preparing _DNN" << endl;
56  _DNN = new DnnPlane(plane_points, verbose);
57 }
58 
59 
60 //----------------------------------------------------------------------
61 /// What on earth does this do?
62 ///
63 /// Example: last true "cylinder" index was 15
64 /// last plane index was 23
65 ///
66 /// Then: _cylinder_index_of_plane_vertex.size() = 24 and
67 /// _mirror_info.size() = 16
68 ///
69 /// IF cylinder_point's phi < pi then
70 /// create: _mirror_info[16] = (main_index = 24, mirror_index=25)
71 /// _cylinder_index_of_plane_vertex[24] = 16
72 /// _cylinder_index_of_plane_vertex[25] = 16
73 /// ELSE
74 /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..)
75 /// _cylinder_index_of_plane_vertex[24] = 16
76 ///
77 /// ADDITIONALLY push the cylinder_point (and if it exists the mirror
78 /// copy) onto the vector plane_points.
79 void Dnn3piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point,
80  vector<EtaPhi> & plane_points) {
81  double phi = cylinder_point.second;
82  assert(phi >= 0.0 && phi < 2*pi);
83 
84  // do main point
85  MirrorVertexInfo mvi;
86  mvi.main_index = _cylinder_index_of_plane_vertex.size();
87  _cylinder_index_of_plane_vertex.push_back(_mirror_info.size());
88  plane_points.push_back(cylinder_point);
89 
90  // do mirror point if need be
91  if (phi < pi) {
92  mvi.mirror_index = _cylinder_index_of_plane_vertex.size();
93  _cylinder_index_of_plane_vertex.push_back(_mirror_info.size());
94  plane_points.push_back(_remap_phi(cylinder_point));
95  } else {
96  mvi.mirror_index = INEXISTENT_VERTEX;
97  }
98 
99  //
100  _mirror_info.push_back(mvi);
101 }
102 
103 
104 //----------------------------------------------------------------------
105 /// insertion and removal of points
106 void Dnn3piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove,
107  const vector<EtaPhi> & points_to_add,
108  vector<int> & indices_added,
109  vector<int> & indices_of_updated_neighbours) {
110 
111  // translate from "cylinder" indices of points to remove to the
112  // plane indices of points to remove, bearing in mind that sometimes
113  // there are multple plane points to remove.
114  vector<int> plane_indices_to_remove;
115  for (unsigned int i=0; i < indices_to_remove.size(); i++) {
116  MirrorVertexInfo * mvi;
117  mvi = & _mirror_info[indices_to_remove[i]];
118  plane_indices_to_remove.push_back(mvi->main_index);
119  if (mvi->mirror_index != INEXISTENT_VERTEX) {
120  plane_indices_to_remove.push_back(mvi->mirror_index);
121  }
122  }
123 
124  // given "cylinder" points to add get hold of the list of
125  // plane-points to add.
126  vector<EtaPhi> plane_points_to_add;
127  indices_added.clear();
128  for (unsigned int i=0; i < points_to_add.size(); i++) {
129  indices_added.push_back(_mirror_info.size());
130  _RegisterCylinderPoint(points_to_add[i], plane_points_to_add);
131  }
132 
133  // now get the hard work done (note that we need to supply the
134  // plane_indices_added vector but that we will not actually check
135  // its contents in any way -- the indices_added that is actually
136  // returned has been calculated above).
137  vector<int> updated_plane_neighbours, plane_indices_added;
138  _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add,
139  plane_indices_added, updated_plane_neighbours);
140 
141  // extract, from the updated_plane_neighbours, the set of cylinder
142  // neighbours that have changed
143  set<int> index_set;
144  unsigned int i;
145  for (i=0; i < updated_plane_neighbours.size(); i++) {
146  index_set.insert(
147  _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);}
148 
149  // decant the set into the vector that needs to be returned
150  indices_of_updated_neighbours.clear();
151  for (set<int>::iterator iter = index_set.begin();
152  iter != index_set.end(); iter++) {
153  indices_of_updated_neighbours.push_back(*iter);
154  }
155 }
156 
157 
158 FASTJET_END_NAMESPACE
159 
160 #endif // DROP_CGAL
161