FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: Dnn3piCylinder.cc 2577 2011-09-13 15:11:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 00030 #ifndef DROP_CGAL // in case we do not have the code for CGAL 00031 #include <set> 00032 #include "fastjet/internal/Dnn3piCylinder.hh" 00033 using namespace std; 00034 00035 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00036 00037 //---------------------------------------------------------------------- 00038 /// initialiser... 00039 Dnn3piCylinder::Dnn3piCylinder( 00040 const vector<EtaPhi> & input_points, 00041 const bool & ignore_nearest_is_mirror, 00042 const bool & verbose) { 00043 00044 _verbose = verbose; 00045 _ignore_nearest_is_mirror = ignore_nearest_is_mirror; 00046 vector<EtaPhi> plane_points; 00047 //plane_points.reserve(2*input_points.size()); 00048 00049 for (unsigned int i=0; i < input_points.size(); i++) { 00050 _RegisterCylinderPoint(input_points[i], plane_points); 00051 } 00052 00053 if (_verbose) cout << "============== Preparing _DNN" << endl; 00054 _DNN = new DnnPlane(plane_points, verbose); 00055 } 00056 00057 00058 //---------------------------------------------------------------------- 00059 /// What on earth does this do? 00060 /// 00061 /// Example: last true "cylinder" index was 15 00062 /// last plane index was 23 00063 /// 00064 /// Then: _cylinder_index_of_plane_vertex.size() = 24 and 00065 /// _mirror_info.size() = 16 00066 /// 00067 /// IF cylinder_point's phi < pi then 00068 /// create: _mirror_info[16] = (main_index = 24, mirror_index=25) 00069 /// _cylinder_index_of_plane_vertex[24] = 16 00070 /// _cylinder_index_of_plane_vertex[25] = 16 00071 /// ELSE 00072 /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..) 00073 /// _cylinder_index_of_plane_vertex[24] = 16 00074 /// 00075 /// ADDITIONALLY push the cylinder_point (and if it exists the mirror 00076 /// copy) onto the vector plane_points. 00077 void Dnn3piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point, 00078 vector<EtaPhi> & plane_points) { 00079 double phi = cylinder_point.second; 00080 assert(phi >= 0.0 && phi < 2*pi); 00081 00082 // do main point 00083 MirrorVertexInfo mvi; 00084 mvi.main_index = _cylinder_index_of_plane_vertex.size(); 00085 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); 00086 plane_points.push_back(cylinder_point); 00087 00088 // do mirror point if need be 00089 if (phi < pi) { 00090 mvi.mirror_index = _cylinder_index_of_plane_vertex.size(); 00091 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); 00092 plane_points.push_back(_remap_phi(cylinder_point)); 00093 } else { 00094 mvi.mirror_index = INEXISTENT_VERTEX; 00095 } 00096 00097 // 00098 _mirror_info.push_back(mvi); 00099 } 00100 00101 00102 //---------------------------------------------------------------------- 00103 /// insertion and removal of points 00104 void Dnn3piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove, 00105 const vector<EtaPhi> & points_to_add, 00106 vector<int> & indices_added, 00107 vector<int> & indices_of_updated_neighbours) { 00108 00109 // translate from "cylinder" indices of points to remove to the 00110 // plane indices of points to remove, bearing in mind that sometimes 00111 // there are multple plane points to remove. 00112 vector<int> plane_indices_to_remove; 00113 for (unsigned int i=0; i < indices_to_remove.size(); i++) { 00114 MirrorVertexInfo * mvi; 00115 mvi = & _mirror_info[indices_to_remove[i]]; 00116 plane_indices_to_remove.push_back(mvi->main_index); 00117 if (mvi->mirror_index != INEXISTENT_VERTEX) { 00118 plane_indices_to_remove.push_back(mvi->mirror_index); 00119 } 00120 } 00121 00122 // given "cylinder" points to add get hold of the list of 00123 // plane-points to add. 00124 vector<EtaPhi> plane_points_to_add; 00125 indices_added.clear(); 00126 for (unsigned int i=0; i < points_to_add.size(); i++) { 00127 indices_added.push_back(_mirror_info.size()); 00128 _RegisterCylinderPoint(points_to_add[i], plane_points_to_add); 00129 } 00130 00131 // now get the hard work done (note that we need to supply the 00132 // plane_indices_added vector but that we will not actually check 00133 // its contents in any way -- the indices_added that is actually 00134 // returned has been calculated above). 00135 vector<int> updated_plane_neighbours, plane_indices_added; 00136 _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add, 00137 plane_indices_added, updated_plane_neighbours); 00138 00139 // extract, from the updated_plane_neighbours, the set of cylinder 00140 // neighbours that have changed 00141 set<int> index_set; 00142 unsigned int i; 00143 for (i=0; i < updated_plane_neighbours.size(); i++) { 00144 index_set.insert( 00145 _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);} 00146 00147 // decant the set into the vector that needs to be returned 00148 indices_of_updated_neighbours.clear(); 00149 for (set<int>::iterator iter = index_set.begin(); 00150 iter != index_set.end(); iter++) { 00151 indices_of_updated_neighbours.push_back(*iter); 00152 } 00153 } 00154 00155 00156 FASTJET_END_NAMESPACE 00157 00158 #endif // DROP_CGAL 00159