[1] | 1 | //STARTHEADER |
---|
| 2 | // $Id: Dnn3piCylinder.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
| 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/Dnn3piCylinder.hh" |
---|
| 33 | using namespace std; |
---|
| 34 | |
---|
| 35 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
| 36 | |
---|
| 37 | //---------------------------------------------------------------------- |
---|
| 38 | /// initialiser... |
---|
| 39 | Dnn3piCylinder::Dnn3piCylinder( |
---|
| 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 | //plane_points.reserve(2*input_points.size()); |
---|
| 48 | |
---|
| 49 | for (unsigned int i=0; i < input_points.size(); i++) { |
---|
| 50 | _RegisterCylinderPoint(input_points[i], plane_points); |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | if (_verbose) cout << "============== Preparing _DNN" << endl; |
---|
| 54 | _DNN = new DnnPlane(plane_points, verbose); |
---|
| 55 | } |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | //---------------------------------------------------------------------- |
---|
| 59 | /// What on earth does this do? |
---|
| 60 | /// |
---|
| 61 | /// Example: last true "cylinder" index was 15 |
---|
| 62 | /// last plane index was 23 |
---|
| 63 | /// |
---|
| 64 | /// Then: _cylinder_index_of_plane_vertex.size() = 24 and |
---|
| 65 | /// _mirror_info.size() = 16 |
---|
| 66 | /// |
---|
| 67 | /// IF cylinder_point's phi < pi then |
---|
| 68 | /// create: _mirror_info[16] = (main_index = 24, mirror_index=25) |
---|
| 69 | /// _cylinder_index_of_plane_vertex[24] = 16 |
---|
| 70 | /// _cylinder_index_of_plane_vertex[25] = 16 |
---|
| 71 | /// ELSE |
---|
| 72 | /// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..) |
---|
| 73 | /// _cylinder_index_of_plane_vertex[24] = 16 |
---|
| 74 | /// |
---|
| 75 | /// ADDITIONALLY push the cylinder_point (and if it exists the mirror |
---|
| 76 | /// copy) onto the vector plane_points. |
---|
| 77 | void Dnn3piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point, |
---|
| 78 | vector<EtaPhi> & plane_points) { |
---|
| 79 | double phi = cylinder_point.second; |
---|
| 80 | assert(phi >= 0.0 && phi < 2*pi); |
---|
| 81 | |
---|
| 82 | // do main point |
---|
| 83 | MirrorVertexInfo mvi; |
---|
| 84 | mvi.main_index = _cylinder_index_of_plane_vertex.size(); |
---|
| 85 | _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); |
---|
| 86 | plane_points.push_back(cylinder_point); |
---|
| 87 | |
---|
| 88 | // do mirror point if need be |
---|
| 89 | if (phi < pi) { |
---|
| 90 | mvi.mirror_index = _cylinder_index_of_plane_vertex.size(); |
---|
| 91 | _cylinder_index_of_plane_vertex.push_back(_mirror_info.size()); |
---|
| 92 | plane_points.push_back(_remap_phi(cylinder_point)); |
---|
| 93 | } else { |
---|
| 94 | mvi.mirror_index = INEXISTENT_VERTEX; |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | // |
---|
| 98 | _mirror_info.push_back(mvi); |
---|
| 99 | } |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | //---------------------------------------------------------------------- |
---|
| 103 | /// insertion and removal of points |
---|
| 104 | void Dnn3piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove, |
---|
| 105 | const vector<EtaPhi> & points_to_add, |
---|
| 106 | vector<int> & indices_added, |
---|
| 107 | vector<int> & indices_of_updated_neighbours) { |
---|
| 108 | |
---|
| 109 | // translate from "cylinder" indices of points to remove to the |
---|
| 110 | // plane indices of points to remove, bearing in mind that sometimes |
---|
| 111 | // there are multple plane points to remove. |
---|
| 112 | vector<int> plane_indices_to_remove; |
---|
| 113 | for (unsigned int i=0; i < indices_to_remove.size(); i++) { |
---|
| 114 | MirrorVertexInfo * mvi; |
---|
| 115 | mvi = & _mirror_info[indices_to_remove[i]]; |
---|
| 116 | plane_indices_to_remove.push_back(mvi->main_index); |
---|
| 117 | if (mvi->mirror_index != INEXISTENT_VERTEX) { |
---|
| 118 | plane_indices_to_remove.push_back(mvi->mirror_index); |
---|
| 119 | } |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | // given "cylinder" points to add get hold of the list of |
---|
| 123 | // plane-points to add. |
---|
| 124 | vector<EtaPhi> plane_points_to_add; |
---|
| 125 | indices_added.clear(); |
---|
| 126 | for (unsigned int i=0; i < points_to_add.size(); i++) { |
---|
| 127 | indices_added.push_back(_mirror_info.size()); |
---|
| 128 | _RegisterCylinderPoint(points_to_add[i], plane_points_to_add); |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | // now get the hard work done (note that we need to supply the |
---|
| 132 | // plane_indices_added vector but that we will not actually check |
---|
| 133 | // its contents in any way -- the indices_added that is actually |
---|
| 134 | // returned has been calculated above). |
---|
| 135 | vector<int> updated_plane_neighbours, plane_indices_added; |
---|
| 136 | _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add, |
---|
| 137 | plane_indices_added, updated_plane_neighbours); |
---|
| 138 | |
---|
| 139 | // extract, from the updated_plane_neighbours, the set of cylinder |
---|
| 140 | // neighbours that have changed |
---|
| 141 | set<int> index_set; |
---|
| 142 | unsigned int i; |
---|
| 143 | for (i=0; i < updated_plane_neighbours.size(); i++) { |
---|
| 144 | index_set.insert( |
---|
| 145 | _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);} |
---|
| 146 | |
---|
| 147 | // decant the set into the vector that needs to be returned |
---|
| 148 | indices_of_updated_neighbours.clear(); |
---|
| 149 | for (set<int>::iterator iter = index_set.begin(); |
---|
| 150 | iter != index_set.end(); iter++) { |
---|
| 151 | indices_of_updated_neighbours.push_back(*iter); |
---|
| 152 | } |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | FASTJET_END_NAMESPACE |
---|
| 157 | |
---|
| 158 | #endif // DROP_CGAL |
---|
| 159 | |
---|