31 #include "fastjet/ClusterSequenceVoronoiArea.hh" 32 #include "fastjet/internal/Voronoi.hh" 43 FASTJET_BEGIN_NAMESPACE
45 typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC;
49 class ClusterSequenceVoronoiArea::VoronoiAreaCalc {
54 VoronoiAreaCalc(
const vector<PseudoJet>::const_iterator &,
55 const vector<PseudoJet>::const_iterator &,
60 inline double area (
int index)
const {
return _areas[index];};
63 std::vector<double> _areas;
65 double _effective_R_squared;
71 double edge_circle_intersection(
const VPoint &p0,
72 const GraphEdge &edge);
77 inline double circle_area(
const double d12_2,
double d01_2,
double d02_2){
78 return 0.5*_effective_R_squared
79 *acos(min(1.0,(d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2))));
88 double VAC::edge_circle_intersection(
const VPoint &p0,
89 const GraphEdge &edge){
90 VPoint p1(edge.x1-p0.x, edge.y1-p0.y);
91 VPoint p2(edge.x2-p0.x, edge.y2-p0.y);
97 double d12_2 =
norm(pdiff);
98 double d01_2 =
norm(p1);
99 double d02_2 =
norm(p2);
102 double delta = d12_2*_effective_R_squared - cross*cross;
106 return circle_area(d12_2, d01_2, d02_2);
121 double tp = (delta-b)/d12_2;
125 return circle_area(d12_2, d01_2, d02_2);
128 double tm = -(delta+b)/d12_2;
140 return tp*0.5*fabs(cross)
141 +circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
147 return (tp-tm)*0.5*fabs(cross)
148 + circle_area(tm*tm*d12_2, d01_2, _effective_R_squared)
149 + circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
156 return circle_area(d12_2, d01_2, d02_2);
160 return 0.5*fabs(cross);
165 return (1-tm)*0.5*fabs(cross)
166 +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
172 VAC::VoronoiAreaCalc(
const vector<PseudoJet>::const_iterator &jet_begin,
173 const vector<PseudoJet>::const_iterator &jet_end,
174 double effective_R) {
176 assert(effective_R < 0.5*pi);
178 vector<VPoint> voronoi_particles;
179 vector<int> voronoi_indices;
181 _effective_R = effective_R;
182 _effective_R_squared = effective_R*effective_R;
184 double minrap = numeric_limits<double>::max();
185 double maxrap = -minrap;
187 unsigned int n_tot = 0, n_added = 0;
191 for (vector<PseudoJet>::const_iterator jet_it = jet_begin;
192 jet_it != jet_end; jet_it++) {
193 _areas.push_back(0.0);
194 if ((jet_it->perp2()) != 0.0 || (jet_it->E() != jet_it->pz())){
196 double rap = jet_it->rap(), phi = jet_it->phi();
197 voronoi_particles.push_back(VPoint(rap, phi));
198 voronoi_indices.push_back(n_tot);
204 if (phi < 2*_effective_R) {
205 voronoi_particles.push_back(VPoint(rap,phi+twopi));
206 voronoi_indices.push_back(-1);
208 }
else if (twopi-phi < 2*_effective_R) {
209 voronoi_particles.push_back(VPoint(rap,phi-twopi));
210 voronoi_indices.push_back(-1);
215 maxrap = max(maxrap,rap);
216 minrap = min(minrap,rap);
222 if (n_added == 0)
return;
226 double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
227 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)-max_extend, pi));
228 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)+max_extend, pi));
229 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi-max_extend));
230 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi+max_extend));
233 VoronoiDiagramGenerator vdg;
234 vdg.generateVoronoi(&voronoi_particles,
235 0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
236 pi-max_extend, pi+max_extend);
240 unsigned int v_index;
242 vector<PseudoJet>::const_iterator jet;
244 while(vdg.getNext(&e)){
246 if (v_index<n_added){
247 p_index = voronoi_indices[v_index];
249 jet = jet_begin+voronoi_indices[v_index];
251 edge_circle_intersection(voronoi_particles[v_index], *e);
255 if (v_index<n_added){
256 p_index = voronoi_indices[v_index];
258 jet = jet_begin+voronoi_indices[v_index];
260 edge_circle_intersection(voronoi_particles[v_index], *e);
271 void ClusterSequenceVoronoiArea::_initializeVA () {
273 _pa_calc =
new VAC(_jets.begin(),
274 _jets.begin()+n_particles(),
275 _effective_Rfact*_jet_def.R());
279 _voronoi_area.reserve(2*n_particles());
280 _voronoi_area_4vector.reserve(2*n_particles());
281 for (
unsigned int i=0; i<n_particles(); i++) {
282 _voronoi_area.push_back(_pa_calc->area(i));
284 if (_jets[i].perp2() > 0) {
285 _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp())
290 _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
295 for (
unsigned int i = n_particles(); i < _history.size(); i++) {
297 PseudoJet area_4vect;
298 if (_history[i].parent2 >= 0) {
299 area_local = _voronoi_area[_history[i].parent1] +
300 _voronoi_area[_history[i].parent2];
301 area_4vect = _voronoi_area_4vector[_history[i].parent1] +
302 _voronoi_area_4vector[_history[i].parent2];
304 area_local = _voronoi_area[_history[i].parent1];
305 area_4vect = _voronoi_area_4vector[_history[i].parent1];
307 _voronoi_area.push_back(area_local);
308 _voronoi_area_4vector.push_back(area_4vect);
314 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
318 FASTJET_END_NAMESPACE
double norm(const VPoint p)
norm of a vector
double scalar_product(const VPoint &p1, const VPoint &p2)
scalar product
double vector_product(const VPoint &p1, const VPoint &p2)
2D vector product