Fork me on GitHub

source: git/external/fastjet/Dnn3piCylinder.cc@ 75f49a7

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 75f49a7 was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 5.8 KB
Line 
1//FJSTARTHEADER
2// $Id: Dnn3piCylinder.cc 3433 2014-07-23 08:17:03Z salam $
3//
4// Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
16// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17// FastJet as part of work towards a scientific publication, please
18// quote the version you use and include a citation to the manual and
19// optionally also to hep-ph/0512210.
20//
21// FastJet is distributed in the hope that it will be useful,
22// but WITHOUT ANY WARRANTY; without even the implied warranty of
23// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24// GNU General Public License for more details.
25//
26// You should have received a copy of the GNU General Public License
27// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28//----------------------------------------------------------------------
29//FJENDHEADER
30
31
32#ifndef DROP_CGAL // in case we do not have the code for CGAL
33#include <set>
34#include "fastjet/internal/Dnn3piCylinder.hh"
35using namespace std;
36
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39//----------------------------------------------------------------------
40/// initialiser...
41Dnn3piCylinder::Dnn3piCylinder(
42 const vector<EtaPhi> & input_points,
43 const bool & ignore_nearest_is_mirror,
44 const bool & verbose) {
45
46 _verbose = verbose;
47 _ignore_nearest_is_mirror = ignore_nearest_is_mirror;
48 vector<EtaPhi> plane_points;
49 //plane_points.reserve(2*input_points.size());
50
51 for (unsigned int i=0; i < input_points.size(); i++) {
52 _RegisterCylinderPoint(input_points[i], plane_points);
53 }
54
55 if (_verbose) cout << "============== Preparing _DNN" << endl;
56 _DNN = new DnnPlane(plane_points, verbose);
57}
58
59
60//----------------------------------------------------------------------
61/// What on earth does this do?
62///
63/// Example: last true "cylinder" index was 15
64/// last plane index was 23
65///
66/// Then: _cylinder_index_of_plane_vertex.size() = 24 and
67/// _mirror_info.size() = 16
68///
69/// IF cylinder_point's phi < pi then
70/// create: _mirror_info[16] = (main_index = 24, mirror_index=25)
71/// _cylinder_index_of_plane_vertex[24] = 16
72/// _cylinder_index_of_plane_vertex[25] = 16
73/// ELSE
74/// create: _mirror_info[16] = (main_index = 24, mirror_index=INEXISTENT..)
75/// _cylinder_index_of_plane_vertex[24] = 16
76///
77/// ADDITIONALLY push the cylinder_point (and if it exists the mirror
78/// copy) onto the vector plane_points.
79void Dnn3piCylinder::_RegisterCylinderPoint (const EtaPhi & cylinder_point,
80 vector<EtaPhi> & plane_points) {
81 double phi = cylinder_point.second;
82 assert(phi >= 0.0 && phi < 2*pi);
83
84 // do main point
85 MirrorVertexInfo mvi;
86 mvi.main_index = _cylinder_index_of_plane_vertex.size();
87 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size());
88 plane_points.push_back(cylinder_point);
89
90 // do mirror point if need be
91 if (phi < pi) {
92 mvi.mirror_index = _cylinder_index_of_plane_vertex.size();
93 _cylinder_index_of_plane_vertex.push_back(_mirror_info.size());
94 plane_points.push_back(_remap_phi(cylinder_point));
95 } else {
96 mvi.mirror_index = INEXISTENT_VERTEX;
97 }
98
99 //
100 _mirror_info.push_back(mvi);
101}
102
103
104//----------------------------------------------------------------------
105/// insertion and removal of points
106void Dnn3piCylinder::RemoveAndAddPoints(const vector<int> & indices_to_remove,
107 const vector<EtaPhi> & points_to_add,
108 vector<int> & indices_added,
109 vector<int> & indices_of_updated_neighbours) {
110
111 // translate from "cylinder" indices of points to remove to the
112 // plane indices of points to remove, bearing in mind that sometimes
113 // there are multple plane points to remove.
114 vector<int> plane_indices_to_remove;
115 for (unsigned int i=0; i < indices_to_remove.size(); i++) {
116 MirrorVertexInfo * mvi;
117 mvi = & _mirror_info[indices_to_remove[i]];
118 plane_indices_to_remove.push_back(mvi->main_index);
119 if (mvi->mirror_index != INEXISTENT_VERTEX) {
120 plane_indices_to_remove.push_back(mvi->mirror_index);
121 }
122 }
123
124 // given "cylinder" points to add get hold of the list of
125 // plane-points to add.
126 vector<EtaPhi> plane_points_to_add;
127 indices_added.clear();
128 for (unsigned int i=0; i < points_to_add.size(); i++) {
129 indices_added.push_back(_mirror_info.size());
130 _RegisterCylinderPoint(points_to_add[i], plane_points_to_add);
131 }
132
133 // now get the hard work done (note that we need to supply the
134 // plane_indices_added vector but that we will not actually check
135 // its contents in any way -- the indices_added that is actually
136 // returned has been calculated above).
137 vector<int> updated_plane_neighbours, plane_indices_added;
138 _DNN->RemoveAndAddPoints(plane_indices_to_remove, plane_points_to_add,
139 plane_indices_added, updated_plane_neighbours);
140
141 // extract, from the updated_plane_neighbours, the set of cylinder
142 // neighbours that have changed
143 set<int> index_set;
144 unsigned int i;
145 for (i=0; i < updated_plane_neighbours.size(); i++) {
146 index_set.insert(
147 _cylinder_index_of_plane_vertex[updated_plane_neighbours[i]]);}
148
149 // decant the set into the vector that needs to be returned
150 indices_of_updated_neighbours.clear();
151 for (set<int>::iterator iter = index_set.begin();
152 iter != index_set.end(); iter++) {
153 indices_of_updated_neighbours.push_back(*iter);
154 }
155}
156
157
158FASTJET_END_NAMESPACE
159
160#endif // DROP_CGAL
161
Note: See TracBrowser for help on using the repository browser.