FastJet 3.0.0
Dnn2piCylinder.cc
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends