fastjet: fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc Class Reference

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

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.

References _areas, _effective_R, _effective_R_squared, edge_circle_intersection(), fastjet::VoronoiDiagramGenerator::generateVoronoi(), fastjet::VoronoiDiagramGenerator::getNext(), fastjet::pi, fastjet::GraphEdge::point1, fastjet::GraphEdge::point2, fastjet::VoronoiDiagramGenerator::resetIterator(), and fastjet::twopi.

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.

References _effective_R_squared, circle_area(), fastjet::norm(), fastjet::scalar_product(), fastjet::vector_product(), fastjet::Point::x, fastjet::GraphEdge::x1, fastjet::GraphEdge::x2, fastjet::Point::y, fastjet::GraphEdge::y1, and fastjet::GraphEdge::y2.

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   }

fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation