FastJet  3.0.6
ClusterSequenceVoronoiArea.cc
1 //STARTHEADER
2 // $Id: ClusterSequenceVoronoiArea.cc 2687 2011-11-14 11:17:51Z soyez $
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
double norm(const VPoint p)
norm of a vector
Definition: Voronoi.hh:150
double scalar_product(const VPoint &p1, const VPoint &p2)
scalar product
Definition: Voronoi.hh:162
double vector_product(const VPoint &p1, const VPoint &p2)
2D vector product
Definition: Voronoi.hh:156