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