FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Dnn2piCylinder.cc
1 //STARTHEADER
2 // $Id: Dnn2piCylinder.cc 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 #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