FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClusterSequenceVoronoiArea.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequenceVoronoiArea.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2006-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 #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
60  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
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  */
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 //----------------------------------------------------------------------
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 ///
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 //----------------------------------------------------------------------
314 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
315  delete _pa_calc;
316 }
317 
318 FASTJET_END_NAMESPACE
double norm(const VPoint p)
norm of a vector
Definition: Voronoi.hh:155
double scalar_product(const VPoint &p1, const VPoint &p2)
scalar product
Definition: Voronoi.hh:167
double vector_product(const VPoint &p1, const VPoint &p2)
2D vector product
Definition: Voronoi.hh:161