29 #include "fastjet/ClusterSequenceVoronoiArea.hh"
30 #include "fastjet/internal/Voronoi.hh"
41 FASTJET_BEGIN_NAMESPACE
43 typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC;
47 class ClusterSequenceVoronoiArea::VoronoiAreaCalc {
52 VoronoiAreaCalc(
const vector<PseudoJet>::const_iterator &,
53 const vector<PseudoJet>::const_iterator &,
58 inline double area (
int index)
const {
return _areas[index];};
61 std::vector<double> _areas;
63 double _effective_R_squared;
69 double edge_circle_intersection(
const VPoint &p0,
70 const GraphEdge &edge);
75 inline double circle_area(
const double d12_2,
double d01_2,
double d02_2){
76 return 0.5*_effective_R_squared
77 *acos(min(1.0,(d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2))));
86 double VAC::edge_circle_intersection(
const VPoint &p0,
87 const GraphEdge &edge){
88 VPoint p1(edge.x1-p0.x, edge.y1-p0.y);
89 VPoint p2(edge.x2-p0.x, edge.y2-p0.y);
95 double d12_2 =
norm(pdiff);
96 double d01_2 =
norm(p1);
97 double d02_2 =
norm(p2);
100 double delta = d12_2*_effective_R_squared - cross*cross;
104 return circle_area(d12_2, d01_2, d02_2);
119 double tp = (delta-b)/d12_2;
123 return circle_area(d12_2, d01_2, d02_2);
126 double tm = -(delta+b)/d12_2;
138 return tp*0.5*fabs(cross)
139 +circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
145 return (tp-tm)*0.5*fabs(cross)
146 + circle_area(tm*tm*d12_2, d01_2, _effective_R_squared)
147 + circle_area((1-tp)*(1-tp)*d12_2, _effective_R_squared, d02_2);
154 return circle_area(d12_2, d01_2, d02_2);
158 return 0.5*fabs(cross);
163 return (1-tm)*0.5*fabs(cross)
164 +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
170 VAC::VoronoiAreaCalc(
const vector<PseudoJet>::const_iterator &jet_begin,
171 const vector<PseudoJet>::const_iterator &jet_end,
172 double effective_R) {
174 assert(effective_R < 0.5*pi);
176 vector<VPoint> voronoi_particles;
177 vector<int> voronoi_indices;
179 _effective_R = effective_R;
180 _effective_R_squared = effective_R*effective_R;
182 double minrap = numeric_limits<double>::max();
183 double maxrap = -minrap;
185 unsigned int n_tot = 0, n_added = 0;
189 for (vector<PseudoJet>::const_iterator jet_it = jet_begin;
190 jet_it != jet_end; jet_it++) {
191 _areas.push_back(0.0);
192 if ((jet_it->perp2()) != 0.0 || (jet_it->E() != jet_it->pz())){
194 double rap = jet_it->rap(), phi = jet_it->phi();
195 voronoi_particles.push_back(VPoint(rap, phi));
196 voronoi_indices.push_back(n_tot);
202 if (phi < 2*_effective_R) {
203 voronoi_particles.push_back(VPoint(rap,phi+twopi));
204 voronoi_indices.push_back(-1);
206 }
else if (twopi-phi < 2*_effective_R) {
207 voronoi_particles.push_back(VPoint(rap,phi-twopi));
208 voronoi_indices.push_back(-1);
213 maxrap = max(maxrap,rap);
214 minrap = min(minrap,rap);
220 if (n_added == 0)
return;
224 double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
225 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)-max_extend, pi));
226 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)+max_extend, pi));
227 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi-max_extend));
228 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi+max_extend));
231 VoronoiDiagramGenerator vdg;
232 vdg.generateVoronoi(&voronoi_particles,
233 0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
234 pi-max_extend, pi+max_extend);
238 unsigned int v_index;
240 vector<PseudoJet>::const_iterator jet;
242 while(vdg.getNext(&e)){
244 if (v_index<n_added){
245 p_index = voronoi_indices[v_index];
247 jet = jet_begin+voronoi_indices[v_index];
249 edge_circle_intersection(voronoi_particles[v_index], *e);
253 if (v_index<n_added){
254 p_index = voronoi_indices[v_index];
256 jet = jet_begin+voronoi_indices[v_index];
258 edge_circle_intersection(voronoi_particles[v_index], *e);
269 void ClusterSequenceVoronoiArea::_initializeVA () {
271 _pa_calc =
new VAC(_jets.begin(),
272 _jets.begin()+n_particles(),
273 _effective_Rfact*_jet_def.R());
277 _voronoi_area.reserve(2*n_particles());
278 for (
unsigned int i=0; i<n_particles(); i++) {
279 _voronoi_area.push_back(_pa_calc->area(i));
281 if (_jets[i].perp2() > 0) {
282 _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp())
287 _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
292 for (
unsigned int i = n_particles(); i < _history.size(); i++) {
294 PseudoJet area_4vect;
295 if (_history[i].parent2 >= 0) {
296 area_local = _voronoi_area[_history[i].parent1] +
297 _voronoi_area[_history[i].parent2];
298 area_4vect = _voronoi_area_4vector[_history[i].parent1] +
299 _voronoi_area_4vector[_history[i].parent2];
301 area_local = _voronoi_area[_history[i].parent1];
302 area_4vect = _voronoi_area_4vector[_history[i].parent1];
304 _voronoi_area.push_back(area_local);
305 _voronoi_area_4vector.push_back(area_4vect);
311 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
315 FASTJET_END_NAMESPACE