Fork me on GitHub

source: git/external/fastjet/Dnn2piCylinder.cc@ ded7435

ImprovedOutputFile Timing dual_readout llp
Last change on this file since ded7435 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 9.4 KB
Line 
1//STARTHEADER
2// $Id$
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"
33using namespace std;
34
35FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
36
37//----------------------------------------------------------------------
38/// initialiser...
39Dnn2piCylinder::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.
71void 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)
105void 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
181void 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
240FASTJET_END_NAMESPACE
241
242#endif // DROP_CGAL
Note: See TracBrowser for help on using the repository browser.