# fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc Class Reference

class for carrying out a voronoi area calculation on a set of initial vectors More...

List of all members.

## Public Member Functions

VoronoiAreaCalc (const vector< PseudoJet >::const_iterator &, const vector< PseudoJet >::const_iterator &, double effective_R)
constructor that takes a range of a vector together with the effective radius for the intersection of discs with voronoi cells
double area (int index) const
return the area of the particle associated with the given index

## Private Member Functions

double edge_circle_intersection (const Point &p0, const GraphEdge &edge)
compute the intersection of one triangle with the circle the area is returned
double circle_area (const double d12_2, double d01_2, double d02_2)
get the area of a circle of radius R centred on the point 0 with 1 and 2 on each "side" of the arc.

## Private Attributes

std::vector< double > _areas
areas, numbered as jets
double _effective_R
effective radius
double _effective_R_squared
effective radius squared

## Detailed Description

class for carrying out a voronoi area calculation on a set of initial vectors

Definition at line 47 of file ClusterSequenceVoronoiArea.cc.

## Constructor & Destructor Documentation

 fastjet::VAC::VoronoiAreaCalc ( const vector< PseudoJet >::const_iterator & , const vector< PseudoJet >::const_iterator & , double effective_R )

constructor that takes a range of a vector together with the effective radius for the intersection of discs with voronoi cells

Definition at line 170 of file ClusterSequenceVoronoiArea.cc.

```00172                                          {
00173
00174   assert(effective_R < 0.5*pi);
00175
00176   vector<Point> voronoi_particles;
00177   vector<int> voronoi_indices;
00178
00179   _effective_R         = effective_R;
00180   _effective_R_squared = effective_R*effective_R;
00181
00182   double minrap = numeric_limits<double>::max();
00183   double maxrap = -minrap;
00184
00185   unsigned int n_tot = 0, n_added = 0;
00186
00187   // loop over jets and create the triangulation, as well as cross-referencing
00188   // info
00189   for (vector<PseudoJet>::const_iterator jet_it = jet_begin;
00190        jet_it != jet_end; jet_it++) {
00191     _areas.push_back(0.0);
00192     if ((jet_it->perp2()) != 0.0 || (jet_it->E() != jet_it->pz())){
00193       // generate the corresponding point
00194       double rap = jet_it->rap(), phi = jet_it->phi();
00195       voronoi_particles.push_back(Point(rap, phi));
00196       voronoi_indices.push_back(n_tot);
00197       n_added++;
00198
00199       // insert a copy of the point if it falls within 2*_R_effective
00200       // of the 0,2pi borders (because we are interested in any
00201       // voronoi edge within _R_effective of the other border)
00202       if (phi < 2*_effective_R) {
00203         voronoi_particles.push_back(Point(rap,phi+twopi));
00204         voronoi_indices.push_back(-1);
00205         n_added++;
00206       } else if (twopi-phi < 2*_effective_R) {
00207         voronoi_particles.push_back(Point(rap,phi-twopi));
00208         voronoi_indices.push_back(-1);
00209         n_added++;
00210       }
00211
00212       // track the rapidity range
00213       maxrap = max(maxrap,rap);
00214       minrap = min(minrap,rap);
00215     }
00216     n_tot++;
00217   }
00218
00219   assert(n_added > 0);
00220
00221   // add extreme cases:
00222   double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
00223   voronoi_particles.push_back(Point(0.5*(minrap+maxrap)-max_extend, M_PI));
00224   voronoi_particles.push_back(Point(0.5*(minrap+maxrap)+max_extend, M_PI));
00225   voronoi_particles.push_back(Point(0.5*(minrap+maxrap), M_PI-max_extend));
00226   voronoi_particles.push_back(Point(0.5*(minrap+maxrap), M_PI+max_extend));
00227
00228   // Build the VD
00229   VoronoiDiagramGenerator vdg;
00230   vdg.generateVoronoi(&voronoi_particles,
00231                       0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
00232                       M_PI-max_extend, M_PI+max_extend);
00233
00234   vdg.resetIterator();
00235   GraphEdge *e=NULL;
00236   unsigned int v_index;
00237   int p_index;
00238   vector<PseudoJet>::const_iterator jet;
00239
00240   while(vdg.getNext(&e)){
00241     v_index = e->point1;
00242     if (v_index<n_added){
00243       p_index = voronoi_indices[v_index];
00244       if (p_index!=-1){
00245         jet = jet_begin+voronoi_indices[v_index];
00246         _areas[p_index]+=
00247           edge_circle_intersection(voronoi_particles[v_index], *e);
00248       }
00249     }
00250     v_index = e->point2;
00251     if (v_index<n_added){
00252       p_index = voronoi_indices[v_index];
00253       if (p_index!=-1){
00254         jet = jet_begin+voronoi_indices[v_index];
00255         _areas[p_index]+=
00256           edge_circle_intersection(voronoi_particles[v_index], *e);
00257       }
00258     }
00259   }
00260
00261 }
```

## Member Function Documentation

 double fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc::area ( int index ) const` [inline]`

return the area of the particle associated with the given index

Definition at line 58 of file ClusterSequenceVoronoiArea.cc.

Referenced by fastjet::ClusterSequenceVoronoiArea::_initializeVA().

```00058 {return _areas[index];};
```

 double fastjet::VAC::edge_circle_intersection ( const Point & p0, const GraphEdge & edge ) ` [private]`

compute the intersection of one triangle with the circle the area is returned

Definition at line 86 of file ClusterSequenceVoronoiArea.cc.

Referenced by VoronoiAreaCalc().

```00087                                                            {
00088   Point p1(edge.x1-p0.x, edge.y1-p0.y);
00089   Point p2(edge.x2-p0.x, edge.y2-p0.y);
00090   Point pdiff = p2-p1;
00091
00092   //fprintf(stdout, "\tpt(%f,%f)\n", p0.x, p0.y);
00093
00094   double cross = vector_product(p1, p2);
00095   double d12_2 = norm(pdiff);
00096   double d01_2 = norm(p1);
00097   double d02_2 = norm(p2);
00098
00099   // compute intersections between edge line and circle
00100   double delta = d12_2*_effective_R_squared - cross*cross;
00101
00102   // if no intersection, area=area_circle
00103   if (delta<=0){
00104     return circle_area(d12_2, d01_2, d02_2);
00105   }
00106
00107   // we'll only need delta's sqrt now
00108   delta = sqrt(delta);
00109
00110   // b is the projection of 01 onto 12
00111   double b = scalar_product(pdiff, p1);
00112
00113   // intersections with the circle:
00114   //   we compute the "coordinate along the line" of the intersection
00115   //   with t=0 (1) corresponding to p1 (p2)
00116   // points with 0<t<1 are within the circle others are outside
00117
00118   // positive intersection
00119   double tp = (delta-b)/d12_2;
00120
00121   // if tp is negative, tm also => inters = circle
00122   if (tp<0)
00123     return circle_area(d12_2, d01_2, d02_2);
00124
00125   // we need the second intersection
00126   double tm = -(delta+b)/d12_2;
00127
00128   // if tp<1, it lies in the circle
00129   if (tp<1){
00130     // if tm<0, the segment has one intersection
00131     // with the circle at p (t=tp)
00132     // the area is a triangle from 1 to p
00133     //        then a circle   from p to 2
00134     // several tricks can be used:
00135     //  - the area of the triangle is tp*area triangle
00136     //  - the lenght for the circle are easily obtained
00137     if (tm<0)
00138       return tp*0.5*fabs(cross)
00139         +circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
00140
00141     // now, 0 < tm < tp < 1
00142     // the segment intersects twice the circle
00143     //   area = 2 cirles at ends + a triangle in the middle
00144     // again, simplifications are staightforward
00145     return (tp-tm)*0.5*fabs(cross)
00146       + circle_area(tm*tm*d12_2, d01_2, _effective_R_squared)
00147       + circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
00148   }
00149
00150   // now, we have tp>1
00151
00152   // if in addition tm>1, intersectino is a circle
00153   if (tm>1)
00154     return circle_area(d12_2, d01_2, d02_2);
00155
00156   // if tm<0, the triangle is inside the circle
00157   if (tm<0)
00158     return 0.5*fabs(cross);
00159
00160   // otherwise, only the "tm point" is on the segment
00161   //   area = circle from 1 to m and triangle from m to 2
00162
00163   return (1-tm)*0.5*fabs(cross)
00164     +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
00165 }
```

 double fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc::circle_area ( const double d12_2, double d01_2, double d02_2 ) ` [inline, private]`

get the area of a circle of radius R centred on the point 0 with 1 and 2 on each "side" of the arc.

dij is the distance between point i and point j and all distances are squared

Definition at line 75 of file ClusterSequenceVoronoiArea.cc.

Referenced by edge_circle_intersection().

```00075                                                                            {
00076     return 0.5*_effective_R_squared
00077       *acos((d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2)));
00078   }
```

## Member Data Documentation

 std::vector fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc::_areas` [private]`

areas, numbered as jets

Definition at line 58 of file ClusterSequenceVoronoiArea.cc.

Referenced by VoronoiAreaCalc().

 double fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc::_effective_R` [private]`

effective radius

Definition at line 62 of file ClusterSequenceVoronoiArea.cc.

Referenced by VoronoiAreaCalc().

 double fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc::_effective_R_squared` [private]`

effective radius squared

Definition at line 63 of file ClusterSequenceVoronoiArea.cc.

Referenced by edge_circle_intersection(), and VoronoiAreaCalc().

The documentation for this class was generated from the following file:
Generated on Tue Dec 18 17:05:52 2007 for fastjet by  1.5.2