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 &	jet_begin,
		const vector< PseudoJet >::const_iterator &	jet_end,
		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 171 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.

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

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 59 of file ClusterSequenceVoronoiArea.cc.

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

00059 {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 87 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().

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

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 76 of file ClusterSequenceVoronoiArea.cc.

Referenced by edge_circle_intersection().

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

fastjet::ClusterSequenceVoronoiArea::VoronoiAreaCalc Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation