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 |
Definition at line 47 of file ClusterSequenceVoronoiArea.cc.
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 }
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 }
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().
effective radius squared
Definition at line 63 of file ClusterSequenceVoronoiArea.cc.
Referenced by edge_circle_intersection(), and VoronoiAreaCalc().