00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00033 #include <set>
00034 #include "fastjet/internal/Dnn2piCylinder.hh"
00035 using namespace std;
00036
00037 FASTJET_BEGIN_NAMESPACE
00038
00039
00041 Dnn2piCylinder::Dnn2piCylinder(
00042 const vector<EtaPhi> & input_points,
00043 const bool & ignore_nearest_is_mirror,
00044 const bool & verbose) {
00045
00046 _verbose = verbose;
00047 _ignore_nearest_is_mirror = ignore_nearest_is_mirror;
00048 vector<EtaPhi> plane_points;
00049 vector<int> plane_point_indices(input_points.size());
00050
00051
00052 for (unsigned int i=0; i < input_points.size(); i++) {
00053 _RegisterCylinderPoint(input_points[i], plane_points);
00054 plane_point_indices[i] = i;
00055 }
00056
00057 if (_verbose) cout << "============== Preparing _DNN" << endl;
00058 _DNN = new DnnPlane(plane_points, verbose);
00059
00060
00061 vector<int> updated_point_indices;
00062 _CreateNecessaryMirrorPoints(plane_point_indices,updated_point_indices);
00063 }
00064
00065
00066
00073 void Dnn2piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point,
00074 vector<EtaPhi> & plane_points) {
00075 double phi = cylinder_point.second;
00076 assert(phi >= 0.0 && phi < 2*pi);
00077
00078
00079 MirrorVertexInfo mvi;
00080 mvi.main_index = _cylinder_index_of_plane_vertex.size();
00081 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size());
00082 plane_points.push_back(cylinder_point);
00083 mvi.mirror_index = INEXISTENT_VERTEX;
00084
00085
00086 _mirror_info.push_back(mvi);
00087 }
00088
00089
00090
00091
00107 void Dnn2piCylinder::_CreateNecessaryMirrorPoints(
00108 const vector<int> & plane_indices,
00109 vector<int> & updated_plane_points) {
00110
00111 vector<EtaPhi> new_plane_points;
00112
00113 for (size_t i = 0; i < plane_indices.size(); i++) {
00114
00115 int ip = plane_indices[i];
00116 EtaPhi position = _DNN->etaphi(ip);
00117 double phi = position.second;
00118
00119
00120
00121
00122
00123 int ic = _cylinder_index_of_plane_vertex[ip];
00124 if (_mirror_info[ic].mirror_index != INEXISTENT_VERTEX) {continue;}
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135 double nndist = _DNN->NearestNeighbourDistance(ip);
00136 if (phi*phi >= nndist && (twopi-phi)*(twopi-phi) >= nndist) {continue;}
00137
00138
00139 new_plane_points.push_back(_remap_phi(position));
00140 _mirror_info[ic].mirror_index = _cylinder_index_of_plane_vertex.size();
00141 _cylinder_index_of_plane_vertex.push_back(ic);
00142 }
00143
00144 vector<int> indices_to_remove;
00145 vector<int> indices_added;
00146 _DNN->RemoveAndAddPoints(indices_to_remove,new_plane_points,indices_added,
00147 updated_plane_points);
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 }
00178
00179
00180
00181
00183 void Dnn2piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove,
00184 const vector<EtaPhi> & points_to_add,
00185 vector<int> & indices_added,
00186 vector<int> & indices_of_updated_neighbours) {
00187
00188
00189
00190
00191 vector<int> plane_indices_to_remove;
00192 for (unsigned int i=0; i < indices_to_remove.size(); i++) {
00193 MirrorVertexInfo * mvi;
00194 mvi = & _mirror_info[indices_to_remove[i]];
00195 plane_indices_to_remove.push_back(mvi->main_index);
00196 if (mvi->mirror_index != INEXISTENT_VERTEX) {
00197 plane_indices_to_remove.push_back(mvi->mirror_index);
00198 }
00199 }
00200
00201
00202
00203 vector<EtaPhi> plane_points_to_add;
00204 indices_added.clear();
00205 for (unsigned int i=0; i < points_to_add.size(); i++) {
00206 indices_added.push_back(_mirror_info.size());
00207 _RegisterCylinderPoint(points_to_add[i], plane_points_to_add);
00208 }
00209
00210
00211
00212
00213
00214 vector<int> updated_plane_neighbours, plane_indices_added;
00215 _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add,
00216 plane_indices_added, updated_plane_neighbours);
00217
00218 vector<int> extra_updated_plane_neighbours;
00219 _CreateNecessaryMirrorPoints(updated_plane_neighbours,
00220 extra_updated_plane_neighbours);
00221
00222
00223
00224
00225 set<int> index_set;
00226 unsigned int i;
00227 for (i=0; i < updated_plane_neighbours.size(); i++) {
00228 index_set.insert(
00229 _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);}
00230 for (i=0; i < extra_updated_plane_neighbours.size(); i++) {
00231 index_set.insert(
00232 _cylinder_index_of_plane_vertex[extra_updated_plane_neighbours[i]]);}
00233
00234
00235 indices_of_updated_neighbours.clear();
00236 for (set<int>::iterator iter = index_set.begin();
00237 iter != index_set.end(); iter++) {
00238 indices_of_updated_neighbours.push_back(*iter);
00239 }
00240 }
00241
00242 FASTJET_END_NAMESPACE
00243
00244 #endif // DROP_CGAL