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.

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   }


Member Data Documentation

std::vector<double> 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  doxygen 1.5.2