1 | //STARTHEADER
|
---|
2 | // $Id$
|
---|
3 | //
|
---|
4 | // Copyright (c) 2006-2007 Matteo Cacciari, Gavin Salam and Gregory Soyez
|
---|
5 | //
|
---|
6 | //----------------------------------------------------------------------
|
---|
7 | // This file is part of a simple command-line handling environment
|
---|
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 | #include "fastjet/ClusterSequenceVoronoiArea.hh"
|
---|
30 | #include "fastjet/internal/Voronoi.hh"
|
---|
31 | #include <list>
|
---|
32 | #include <cassert>
|
---|
33 | #include <ostream>
|
---|
34 | #include <fstream>
|
---|
35 | #include <iterator>
|
---|
36 | #include <cmath>
|
---|
37 | #include <limits>
|
---|
38 |
|
---|
39 | using namespace std;
|
---|
40 |
|
---|
41 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
---|
42 |
|
---|
43 | typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC;
|
---|
44 |
|
---|
45 | /// class for carrying out a voronoi area calculation on a set of
|
---|
46 | /// initial vectors
|
---|
47 | class ClusterSequenceVoronoiArea::VoronoiAreaCalc {
|
---|
48 | public:
|
---|
49 | /// constructor that takes a range of a vector together with the
|
---|
50 | /// effective radius for the intersection of discs with voronoi
|
---|
51 | /// cells
|
---|
52 | VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &,
|
---|
53 | const vector<PseudoJet>::const_iterator &,
|
---|
54 | double effective_R);
|
---|
55 |
|
---|
56 | /// return the area of the particle associated with the given
|
---|
57 | /// index
|
---|
58 | inline double area (int index) const {return _areas[index];};
|
---|
59 |
|
---|
60 | private:
|
---|
61 | std::vector<double> _areas; ///< areas, numbered as jets
|
---|
62 | double _effective_R; ///< effective radius
|
---|
63 | double _effective_R_squared; ///< effective radius squared
|
---|
64 |
|
---|
65 | /**
|
---|
66 | * compute the intersection of one triangle with the circle
|
---|
67 | * the area is returned
|
---|
68 | */
|
---|
69 | double edge_circle_intersection(const VPoint &p0,
|
---|
70 | const GraphEdge &edge);
|
---|
71 |
|
---|
72 | /// get the area of a circle of radius R centred on the point 0 with
|
---|
73 | /// 1 and 2 on each "side" of the arc. dij is the distance between
|
---|
74 | /// point i and point j and all distances are squared
|
---|
75 | inline double circle_area(const double d12_2, double d01_2, double d02_2){
|
---|
76 | return 0.5*_effective_R_squared
|
---|
77 | *acos(min(1.0,(d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2))));
|
---|
78 | }
|
---|
79 | };
|
---|
80 |
|
---|
81 |
|
---|
82 | /**
|
---|
83 | * compute the intersection of one triangle with the circle
|
---|
84 | * the area is returned
|
---|
85 | */
|
---|
86 | double VAC::edge_circle_intersection(const VPoint &p0,
|
---|
87 | const GraphEdge &edge){
|
---|
88 | VPoint p1(edge.x1-p0.x, edge.y1-p0.y);
|
---|
89 | VPoint p2(edge.x2-p0.x, edge.y2-p0.y);
|
---|
90 | VPoint pdiff = p2-p1;
|
---|
91 |
|
---|
92 | //fprintf(stdout, "\tpt(%f,%f)\n", p0.x, p0.y);
|
---|
93 |
|
---|
94 | double cross = vector_product(p1, p2);
|
---|
95 | double d12_2 = norm(pdiff);
|
---|
96 | double d01_2 = norm(p1);
|
---|
97 | double d02_2 = norm(p2);
|
---|
98 |
|
---|
99 | // compute intersections between edge line and circle
|
---|
100 | double delta = d12_2*_effective_R_squared - cross*cross;
|
---|
101 |
|
---|
102 | // if no intersection, area=area_circle
|
---|
103 | if (delta<=0){
|
---|
104 | return circle_area(d12_2, d01_2, d02_2);
|
---|
105 | }
|
---|
106 |
|
---|
107 | // we'll only need delta's sqrt now
|
---|
108 | delta = sqrt(delta);
|
---|
109 |
|
---|
110 | // b is the projection of 01 onto 12
|
---|
111 | double b = scalar_product(pdiff, p1);
|
---|
112 |
|
---|
113 | // intersections with the circle:
|
---|
114 | // we compute the "coordinate along the line" of the intersection
|
---|
115 | // with t=0 (1) corresponding to p1 (p2)
|
---|
116 | // points with 0<t<1 are within the circle others are outside
|
---|
117 |
|
---|
118 | // positive intersection
|
---|
119 | double tp = (delta-b)/d12_2;
|
---|
120 |
|
---|
121 | // if tp is negative, tm also => inters = circle
|
---|
122 | if (tp<0)
|
---|
123 | return circle_area(d12_2, d01_2, d02_2);
|
---|
124 |
|
---|
125 | // we need the second intersection
|
---|
126 | double tm = -(delta+b)/d12_2;
|
---|
127 |
|
---|
128 | // if tp<1, it lies in the circle
|
---|
129 | if (tp<1){
|
---|
130 | // if tm<0, the segment has one intersection
|
---|
131 | // with the circle at p (t=tp)
|
---|
132 | // the area is a triangle from 1 to p
|
---|
133 | // then a circle from p to 2
|
---|
134 | // several tricks can be used:
|
---|
135 | // - the area of the triangle is tp*area triangle
|
---|
136 | // - the lenght for the circle are easily obtained
|
---|
137 | if (tm<0)
|
---|
138 | return tp*0.5*fabs(cross)
|
---|
139 | +circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
|
---|
140 |
|
---|
141 | // now, 0 < tm < tp < 1
|
---|
142 | // the segment intersects twice the circle
|
---|
143 | // area = 2 cirles at ends + a triangle in the middle
|
---|
144 | // again, simplifications are staightforward
|
---|
145 | return (tp-tm)*0.5*fabs(cross)
|
---|
146 | + circle_area(tm*tm*d12_2, d01_2, _effective_R_squared)
|
---|
147 | + circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
|
---|
148 | }
|
---|
149 |
|
---|
150 | // now, we have tp>1
|
---|
151 |
|
---|
152 | // if in addition tm>1, intersectino is a circle
|
---|
153 | if (tm>1)
|
---|
154 | return circle_area(d12_2, d01_2, d02_2);
|
---|
155 |
|
---|
156 | // if tm<0, the triangle is inside the circle
|
---|
157 | if (tm<0)
|
---|
158 | return 0.5*fabs(cross);
|
---|
159 |
|
---|
160 | // otherwise, only the "tm point" is on the segment
|
---|
161 | // area = circle from 1 to m and triangle from m to 2
|
---|
162 |
|
---|
163 | return (1-tm)*0.5*fabs(cross)
|
---|
164 | +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
|
---|
165 | }
|
---|
166 |
|
---|
167 |
|
---|
168 | // the constructor...
|
---|
169 | //----------------------------------------------------------------------
|
---|
170 | VAC::VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &jet_begin,
|
---|
171 | const vector<PseudoJet>::const_iterator &jet_end,
|
---|
172 | double effective_R) {
|
---|
173 |
|
---|
174 | assert(effective_R < 0.5*pi);
|
---|
175 |
|
---|
176 | vector<VPoint> voronoi_particles;
|
---|
177 | vector<int> voronoi_indices;
|
---|
178 |
|
---|
179 | _effective_R = effective_R;
|
---|
180 | _effective_R_squared = effective_R*effective_R;
|
---|
181 |
|
---|
182 | double minrap = numeric_limits<double>::max();
|
---|
183 | double maxrap = -minrap;
|
---|
184 |
|
---|
185 | unsigned int n_tot = 0, n_added = 0;
|
---|
186 |
|
---|
187 | // loop over jets and create the triangulation, as well as cross-referencing
|
---|
188 | // info
|
---|
189 | for (vector<PseudoJet>::const_iterator jet_it = jet_begin;
|
---|
190 | jet_it != jet_end; jet_it++) {
|
---|
191 | _areas.push_back(0.0);
|
---|
192 | if ((jet_it->perp2()) != 0.0 || (jet_it->E() != jet_it->pz())){
|
---|
193 | // generate the corresponding point
|
---|
194 | double rap = jet_it->rap(), phi = jet_it->phi();
|
---|
195 | voronoi_particles.push_back(VPoint(rap, phi));
|
---|
196 | voronoi_indices.push_back(n_tot);
|
---|
197 | n_added++;
|
---|
198 |
|
---|
199 | // insert a copy of the point if it falls within 2*_R_effective
|
---|
200 | // of the 0,2pi borders (because we are interested in any
|
---|
201 | // voronoi edge within _R_effective of the other border)
|
---|
202 | if (phi < 2*_effective_R) {
|
---|
203 | voronoi_particles.push_back(VPoint(rap,phi+twopi));
|
---|
204 | voronoi_indices.push_back(-1);
|
---|
205 | n_added++;
|
---|
206 | } else if (twopi-phi < 2*_effective_R) {
|
---|
207 | voronoi_particles.push_back(VPoint(rap,phi-twopi));
|
---|
208 | voronoi_indices.push_back(-1);
|
---|
209 | n_added++;
|
---|
210 | }
|
---|
211 |
|
---|
212 | // track the rapidity range
|
---|
213 | maxrap = max(maxrap,rap);
|
---|
214 | minrap = min(minrap,rap);
|
---|
215 | }
|
---|
216 | n_tot++;
|
---|
217 | }
|
---|
218 |
|
---|
219 | // allow for 0-particle case in graceful way
|
---|
220 | if (n_added == 0) return;
|
---|
221 | // assert(n_added > 0); // old (pre 2.4) non-graceful exit
|
---|
222 |
|
---|
223 | // add extreme cases (corner particles):
|
---|
224 | double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
|
---|
225 | voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)-max_extend, pi));
|
---|
226 | voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)+max_extend, pi));
|
---|
227 | voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi-max_extend));
|
---|
228 | voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi+max_extend));
|
---|
229 |
|
---|
230 | // Build the VD
|
---|
231 | VoronoiDiagramGenerator vdg;
|
---|
232 | vdg.generateVoronoi(&voronoi_particles,
|
---|
233 | 0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
|
---|
234 | pi-max_extend, pi+max_extend);
|
---|
235 |
|
---|
236 | vdg.resetIterator();
|
---|
237 | GraphEdge *e=NULL;
|
---|
238 | unsigned int v_index;
|
---|
239 | int p_index;
|
---|
240 | vector<PseudoJet>::const_iterator jet;
|
---|
241 |
|
---|
242 | while(vdg.getNext(&e)){
|
---|
243 | v_index = e->point1;
|
---|
244 | if (v_index<n_added){ // this removes the corner particles
|
---|
245 | p_index = voronoi_indices[v_index];
|
---|
246 | if (p_index!=-1){ // this removes the copies
|
---|
247 | jet = jet_begin+voronoi_indices[v_index];
|
---|
248 | _areas[p_index]+=
|
---|
249 | edge_circle_intersection(voronoi_particles[v_index], *e);
|
---|
250 | }
|
---|
251 | }
|
---|
252 | v_index = e->point2;
|
---|
253 | if (v_index<n_added){ // this removes the corner particles
|
---|
254 | p_index = voronoi_indices[v_index];
|
---|
255 | if (p_index!=-1){ // this removes the copies
|
---|
256 | jet = jet_begin+voronoi_indices[v_index];
|
---|
257 | _areas[p_index]+=
|
---|
258 | edge_circle_intersection(voronoi_particles[v_index], *e);
|
---|
259 | }
|
---|
260 | }
|
---|
261 | }
|
---|
262 |
|
---|
263 |
|
---|
264 | }
|
---|
265 |
|
---|
266 |
|
---|
267 | //----------------------------------------------------------------------
|
---|
268 | ///
|
---|
269 | void ClusterSequenceVoronoiArea::_initializeVA () {
|
---|
270 | // run the VAC on our original particles
|
---|
271 | _pa_calc = new VAC(_jets.begin(),
|
---|
272 | _jets.begin()+n_particles(),
|
---|
273 | _effective_Rfact*_jet_def.R());
|
---|
274 |
|
---|
275 | // transfer the areas to our local structure
|
---|
276 | // -- first the initial ones
|
---|
277 | _voronoi_area.reserve(2*n_particles());
|
---|
278 | for (unsigned int i=0; i<n_particles(); i++) {
|
---|
279 | _voronoi_area.push_back(_pa_calc->area(i));
|
---|
280 | // make a stab at a 4-vector area
|
---|
281 | if (_jets[i].perp2() > 0) {
|
---|
282 | _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp())
|
---|
283 | * _jets[i]);
|
---|
284 | } else {
|
---|
285 | // not sure what to do here -- just put zero (it won't be meaningful
|
---|
286 | // anyway)
|
---|
287 | _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
|
---|
288 | }
|
---|
289 | }
|
---|
290 |
|
---|
291 | // -- then the combined areas that arise from the clustering
|
---|
292 | for (unsigned int i = n_particles(); i < _history.size(); i++) {
|
---|
293 | double area_local;
|
---|
294 | PseudoJet area_4vect;
|
---|
295 | if (_history[i].parent2 >= 0) {
|
---|
296 | area_local = _voronoi_area[_history[i].parent1] +
|
---|
297 | _voronoi_area[_history[i].parent2];
|
---|
298 | area_4vect = _voronoi_area_4vector[_history[i].parent1] +
|
---|
299 | _voronoi_area_4vector[_history[i].parent2];
|
---|
300 | } else {
|
---|
301 | area_local = _voronoi_area[_history[i].parent1];
|
---|
302 | area_4vect = _voronoi_area_4vector[_history[i].parent1];
|
---|
303 | }
|
---|
304 | _voronoi_area.push_back(area_local);
|
---|
305 | _voronoi_area_4vector.push_back(area_4vect);
|
---|
306 | }
|
---|
307 |
|
---|
308 | }
|
---|
309 |
|
---|
310 | //----------------------------------------------------------------------
|
---|
311 | ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
|
---|
312 | delete _pa_calc;
|
---|
313 | }
|
---|
314 |
|
---|
315 | FASTJET_END_NAMESPACE
|
---|