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 49 of file ClusterSequenceVoronoiArea.cc.
|
||||||||||||||||
|
constructor that takes a range of a vector together with the effective radius for the intersection of discs with voronoi cells
Definition at line 172 of file ClusterSequenceVoronoiArea.cc. References _areas, _effective_R, _effective_R_squared, edge_circle_intersection(), fastjet::VoronoiDiagramGenerator::generateVoronoi(), fastjet::VoronoiDiagramGenerator::getNext(), fastjet::pi, fastjet::VoronoiDiagramGenerator::resetIterator(), and fastjet::twopi. 00174 {
00175
00176 assert(effective_R < 0.5*pi);
00177
00178 vector<Point> voronoi_particles;
00179 vector<int> voronoi_indices;
00180
00181 _effective_R = effective_R;
00182 _effective_R_squared = effective_R*effective_R;
00183
00184 double minrap = numeric_limits<double>::max();
00185 double maxrap = -minrap;
00186
00187 unsigned int n_tot = 0, n_added = 0;
00188
00189 // loop over jets and create the triangulation, as well as cross-referencing
00190 // info
00191 for (vector<PseudoJet>::const_iterator jet_it = jet_begin;
00192 jet_it != jet_end; jet_it++) {
00193 _areas.push_back(0.0);
00194 if ((jet_it->perp2()) != 0.0 || (jet_it->E() != jet_it->pz())){
00195 // generate the corresponding point
00196 double rap = jet_it->rap(), phi = jet_it->phi();
00197 voronoi_particles.push_back(Point(rap, phi));
00198 voronoi_indices.push_back(n_tot);
00199 n_added++;
00200
00201 // insert a copy of the point if it falls within 2*_R_effective
00202 // of the 0,2pi borders (because we are interested in any
00203 // voronoi edge within _R_effective of the other border)
00204 if (phi < 2*_effective_R) {
00205 voronoi_particles.push_back(Point(rap,phi+twopi));
00206 voronoi_indices.push_back(-1);
00207 n_added++;
00208 } else if (twopi-phi < 2*_effective_R) {
00209 voronoi_particles.push_back(Point(rap,phi-twopi));
00210 voronoi_indices.push_back(-1);
00211 n_added++;
00212 }
00213
00214 // track the rapidity range
00215 maxrap = max(maxrap,rap);
00216 minrap = min(minrap,rap);
00217 }
00218 n_tot++;
00219 }
00220
00221 assert(n_added > 0);
00222
00223 // add extreme cases (corner particles):
00224 double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
00225 voronoi_particles.push_back(Point(0.5*(minrap+maxrap)-max_extend, pi));
00226 voronoi_particles.push_back(Point(0.5*(minrap+maxrap)+max_extend, pi));
00227 voronoi_particles.push_back(Point(0.5*(minrap+maxrap), pi-max_extend));
00228 voronoi_particles.push_back(Point(0.5*(minrap+maxrap), pi+max_extend));
00229
00230 // Build the VD
00231 VoronoiDiagramGenerator vdg;
00232 vdg.generateVoronoi(&voronoi_particles,
00233 0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
00234 pi-max_extend, pi+max_extend);
00235
00236 vdg.resetIterator();
00237 GraphEdge *e=NULL;
00238 unsigned int v_index;
00239 int p_index;
00240 vector<PseudoJet>::const_iterator jet;
00241
00242 while(vdg.getNext(&e)){
00243 v_index = e->point1;
00244 if (v_index<n_added){ // this removes the corner particles
00245 p_index = voronoi_indices[v_index];
00246 if (p_index!=-1){ // this removes the copies
00247 jet = jet_begin+voronoi_indices[v_index];
00248 _areas[p_index]+=
00249 edge_circle_intersection(voronoi_particles[v_index], *e);
00250 }
00251 }
00252 v_index = e->point2;
00253 if (v_index<n_added){ // this removes the corner particles
00254 p_index = voronoi_indices[v_index];
00255 if (p_index!=-1){ // this removes the copies
00256 jet = jet_begin+voronoi_indices[v_index];
00257 _areas[p_index]+=
00258 edge_circle_intersection(voronoi_particles[v_index], *e);
00259 }
00260 }
00261 }
00262
00263
00264 }
|
|
|
return the area of the particle associated with the given index
Definition at line 60 of file ClusterSequenceVoronoiArea.cc. Referenced by fastjet::ClusterSequenceVoronoiArea::_initializeVA(). 00060 {return _areas[index];};
|
|
||||||||||||||||
|
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 77 of file ClusterSequenceVoronoiArea.cc. Referenced by edge_circle_intersection(). 00077 {
00078 return 0.5*_effective_R_squared
00079 *acos((d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2)));
00080 }
|
|
||||||||||||
|
compute the intersection of one triangle with the circle the area is returned
Definition at line 88 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(). 00089 {
00090 Point p1(edge.x1-p0.x, edge.y1-p0.y);
00091 Point p2(edge.x2-p0.x, edge.y2-p0.y);
00092 Point pdiff = p2-p1;
00093
00094 //fprintf(stdout, "\tpt(%f,%f)\n", p0.x, p0.y);
00095
00096 double cross = vector_product(p1, p2);
00097 double d12_2 = norm(pdiff);
00098 double d01_2 = norm(p1);
00099 double d02_2 = norm(p2);
00100
00101 // compute intersections between edge line and circle
00102 double delta = d12_2*_effective_R_squared - cross*cross;
00103
00104 // if no intersection, area=area_circle
00105 if (delta<=0){
00106 return circle_area(d12_2, d01_2, d02_2);
00107 }
00108
00109 // we'll only need delta's sqrt now
00110 delta = sqrt(delta);
00111
00112 // b is the projection of 01 onto 12
00113 double b = scalar_product(pdiff, p1);
00114
00115 // intersections with the circle:
00116 // we compute the "coordinate along the line" of the intersection
00117 // with t=0 (1) corresponding to p1 (p2)
00118 // points with 0<t<1 are within the circle others are outside
00119
00120 // positive intersection
00121 double tp = (delta-b)/d12_2;
00122
00123 // if tp is negative, tm also => inters = circle
00124 if (tp<0)
00125 return circle_area(d12_2, d01_2, d02_2);
00126
00127 // we need the second intersection
00128 double tm = -(delta+b)/d12_2;
00129
00130 // if tp<1, it lies in the circle
00131 if (tp<1){
00132 // if tm<0, the segment has one intersection
00133 // with the circle at p (t=tp)
00134 // the area is a triangle from 1 to p
00135 // then a circle from p to 2
00136 // several tricks can be used:
00137 // - the area of the triangle is tp*area triangle
00138 // - the lenght for the circle are easily obtained
00139 if (tm<0)
00140 return tp*0.5*fabs(cross)
00141 +circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
00142
00143 // now, 0 < tm < tp < 1
00144 // the segment intersects twice the circle
00145 // area = 2 cirles at ends + a triangle in the middle
00146 // again, simplifications are staightforward
00147 return (tp-tm)*0.5*fabs(cross)
00148 + circle_area(tm*tm*d12_2, d01_2, _effective_R_squared)
00149 + circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
00150 }
00151
00152 // now, we have tp>1
00153
00154 // if in addition tm>1, intersectino is a circle
00155 if (tm>1)
00156 return circle_area(d12_2, d01_2, d02_2);
00157
00158 // if tm<0, the triangle is inside the circle
00159 if (tm<0)
00160 return 0.5*fabs(cross);
00161
00162 // otherwise, only the "tm point" is on the segment
00163 // area = circle from 1 to m and triangle from m to 2
00164
00165 return (1-tm)*0.5*fabs(cross)
00166 +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
00167 }
|
|
|
areas, numbered as jets
Definition at line 60 of file ClusterSequenceVoronoiArea.cc. Referenced by VoronoiAreaCalc(). |
|
|
effective radius
Definition at line 64 of file ClusterSequenceVoronoiArea.cc. Referenced by VoronoiAreaCalc(). |
|
|
effective radius squared
Definition at line 65 of file ClusterSequenceVoronoiArea.cc. Referenced by edge_circle_intersection(), and VoronoiAreaCalc(). |
1.4.2