FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: Dnn2piCylinder.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/Dnn2piCylinder.hh" 00033 using namespace std; 00034 00035 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00036 00037 //---------------------------------------------------------------------- 00038 /// initialiser... 00039 Dnn2piCylinder::Dnn2piCylinder( 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 vector<int> plane_point_indices(input_points.size()); 00048 //plane_points.reserve(2*input_points.size()); 00049 00050 for (unsigned int i=0; i < input_points.size(); i++) { 00051 _RegisterCylinderPoint(input_points[i], plane_points); 00052 plane_point_indices[i] = i; 00053 } 00054 00055 if (_verbose) cout << "============== Preparing _DNN" << endl; 00056 _DNN = new DnnPlane(plane_points, verbose); 00057 00058 00059 vector<int> updated_point_indices; // we'll not use information from this 00060 _CreateNecessaryMirrorPoints(plane_point_indices,updated_point_indices); 00061 } 00062 00063 00064 //---------------------------------------------------------------------- 00065 /// Actions here are similar to those in the 00066 /// Dnn3piCylinder::_RegisterCylinderPoint case, however here we do 00067 /// NOT create the mirror point -- instead we initialise the structure 00068 /// as if there were no need for the mirror point. 00069 /// 00070 /// ADDITIONALLY push the cylinder_point onto the vector plane_points. 00071 void Dnn2piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point, 00072 vector<EtaPhi> & plane_points) { 00073 double phi = cylinder_point.second; 00074 assert(phi >= 0.0 && phi < 2*pi); 00075 00076 // do main point 00077 MirrorVertexInfo mvi; 00078 mvi.main_index = _cylinder_index_of_plane_vertex.size(); 00079 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); 00080 plane_points.push_back(cylinder_point); 00081 mvi.mirror_index = INEXISTENT_VERTEX; 00082 00083 // 00084 _mirror_info.push_back(mvi); 00085 } 00086 00087 00088 00089 //---------------------------------------------------------------------- 00090 /// For each plane point specified in the vector plane_indices, 00091 /// establish whether there is a need to create a mirror point 00092 /// according to the following criteria: 00093 /// 00094 /// . phi < pi 00095 /// . mirror does not already exist 00096 /// . phi < NearestNeighbourDistance 00097 /// (if this is not true then there is no way that its mirror point 00098 /// could have a nearer neighbour). 00099 /// 00100 /// If conditions all hold, then create the mirror point, insert it 00101 /// into the _DNN structure, adjusting any nearest neighbours, and 00102 /// return the list of plane points whose nearest neighbours have 00103 /// changed (this will include the new neighbours that have just been 00104 /// added) 00105 void Dnn2piCylinder::_CreateNecessaryMirrorPoints( 00106 const vector<int> & plane_indices, 00107 vector<int> & updated_plane_points) { 00108 00109 vector<EtaPhi> new_plane_points; 00110 00111 for (size_t i = 0; i < plane_indices.size(); i++) { 00112 00113 int ip = plane_indices[i]; // plane index 00114 EtaPhi position = _DNN->etaphi(ip); 00115 double phi = position.second; 00116 00117 //BAD // require phi < pi 00118 //BAD if (phi >= pi) {continue;} 00119 00120 // require absence of mirror 00121 int ic = _cylinder_index_of_plane_vertex[ip]; 00122 if (_mirror_info[ic].mirror_index != INEXISTENT_VERTEX) {continue;} 00123 00124 //printf("%3d %3d %10.5f %10.5f %3d\n",ic, ip, phi, 00125 // _DNN->NearestNeighbourDistance(ip),_DNN->NearestNeighbourIndex(ip)); 00126 00127 00128 // check that we are sufficiently close to the border -- 00129 // i.e. closer than nearest neighbour distance. But RECALL: 00130 // nearest neighbourDistance actually returns a squared distance 00131 // (this was really stupid on our part -- caused considerable loss 00132 // of time ... ) 00133 double nndist = _DNN->NearestNeighbourDistance(ip); 00134 if (phi*phi >= nndist && (twopi-phi)*(twopi-phi) >= nndist) {continue;} 00135 00136 // now proceed to prepare the point for addition 00137 new_plane_points.push_back(_remap_phi(position)); 00138 _mirror_info[ic].mirror_index = _cylinder_index_of_plane_vertex.size(); 00139 _cylinder_index_of_plane_vertex.push_back(ic); 00140 } 00141 00142 vector<int> indices_to_remove; // empty 00143 vector<int> indices_added; // will be filled as result of call 00144 _DNN->RemoveAndAddPoints(indices_to_remove,new_plane_points,indices_added, 00145 updated_plane_points); 00146 00147 // occasionally, addition of points might cause a change in the 00148 // nearest neighbour of a point in the 0--pi range? (But should be 00149 // impossible -- we add points beyond 2pi, so they can only be 00150 // nearest neighbours of points in the range pi--2pi; there is one 00151 // exception -- the nearest neighbour of one's self -- but in that 00152 // case this will already have been discovered, so there should be 00153 // nothing left to do). 00154 00155 // To be on the safe side, check to see if we have updated points 00156 // with phi<pi and no current mirror point. BUT: this check, as 00157 // written, only makes sense when the mirror image was created only 00158 // beyond 2pi, which is no longer the case. Only apparent 00159 // alternative is to run separate insertions for beyond 2pi and 00160 // below phi=0, with separate checks in the two cases. But, given 00161 // that across a large number of recombinations and events in the 00162 // old method (single mirror), we never ran into this problem, it's 00163 // probably safe to assume that the arguments given above are OK! So 00164 // comment out the check... 00165 //NOTNEEDED for (size_t i = 0; i < updated_plane_points.size(); i++) { 00166 //NOTNEEDED int ip = updated_plane_points[i]; 00167 //NOTNEEDED double phi = _DNN->phi(ip); 00168 //NOTNEEDED int ic = _cylinder_index_of_plane_vertex[ip]; 00169 //NOTNEEDED assert (!(phi < pi && _mirror_info[ic].mirror_index == INEXISTENT_VERTEX)); 00170 //NOTNEEDED } 00171 // alternative recursive code 00172 //vector<int> extra_updated_plane_points; 00173 //_CreateNecessaryMirrorPoints(updated_plane_points,extra_updated_plane_points) 00174 //updated_plane_points.push_back(extra_updated_plane_points); 00175 } 00176 00177 00178 00179 //---------------------------------------------------------------------- 00180 /// insertion and removal of points 00181 void Dnn2piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove, 00182 const vector<EtaPhi> & points_to_add, 00183 vector<int> & indices_added, 00184 vector<int> & indices_of_updated_neighbours) { 00185 00186 // translate from "cylinder" indices of points to remove to the 00187 // plane indices of points to remove, bearing in mind that sometimes 00188 // there are multple plane points to remove. 00189 vector<int> plane_indices_to_remove; 00190 for (unsigned int i=0; i < indices_to_remove.size(); i++) { 00191 MirrorVertexInfo * mvi; 00192 mvi = & _mirror_info[indices_to_remove[i]]; 00193 plane_indices_to_remove.push_back(mvi->main_index); 00194 if (mvi->mirror_index != INEXISTENT_VERTEX) { 00195 plane_indices_to_remove.push_back(mvi->mirror_index); 00196 } 00197 } 00198 00199 // given "cylinder" points to add get hold of the list of 00200 // plane-points to add. 00201 vector<EtaPhi> plane_points_to_add; 00202 indices_added.clear(); 00203 for (unsigned int i=0; i < points_to_add.size(); i++) { 00204 indices_added.push_back(_mirror_info.size()); 00205 _RegisterCylinderPoint(points_to_add[i], plane_points_to_add); 00206 } 00207 00208 // now get the hard work done (note that we need to supply the 00209 // plane_indices_added vector but that we will not actually check 00210 // its contents in any way -- the indices_added that is actually 00211 // returned has been calculated above). 00212 vector<int> updated_plane_neighbours, plane_indices_added; 00213 _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add, 00214 plane_indices_added, updated_plane_neighbours); 00215 00216 vector<int> extra_updated_plane_neighbours; 00217 _CreateNecessaryMirrorPoints(updated_plane_neighbours, 00218 extra_updated_plane_neighbours); 00219 00220 // extract, from the updated_plane_neighbours, and 00221 // extra_updated_plane_neighbours, the set of cylinder neighbours 00222 // that have changed 00223 set<int> index_set; 00224 unsigned int i; 00225 for (i=0; i < updated_plane_neighbours.size(); i++) { 00226 index_set.insert( 00227 _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);} 00228 for (i=0; i < extra_updated_plane_neighbours.size(); i++) { 00229 index_set.insert( 00230 _cylinder_index_of_plane_vertex[extra_updated_plane_neighbours[i]]);} 00231 00232 // decant the set into the vector that needs to be returned 00233 indices_of_updated_neighbours.clear(); 00234 for (set<int>::iterator iter = index_set.begin(); 00235 iter != index_set.end(); iter++) { 00236 indices_of_updated_neighbours.push_back(*iter); 00237 } 00238 } 00239 00240 FASTJET_END_NAMESPACE 00241 00242 #endif // DROP_CGAL