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