00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/ClusterSequenceVoronoiArea.hh"
00032 #include "fastjet/internal/Voronoi.hh"
00033 #include <list>
00034 #include <cassert>
00035 #include <ostream>
00036 #include <fstream>
00037 #include <iterator>
00038 #include <cmath>
00039 #include <limits>
00040
00041 using namespace std;
00042
00043 FASTJET_BEGIN_NAMESPACE
00044
00045 typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC;
00046
00049 class ClusterSequenceVoronoiArea::VoronoiAreaCalc {
00050 public:
00054 VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &,
00055 const vector<PseudoJet>::const_iterator &,
00056 double effective_R);
00057
00060 inline double area (int index) const {return _areas[index];};
00061
00062 private:
00063 std::vector<double> _areas;
00064 double _effective_R;
00065 double _effective_R_squared;
00066
00071 double edge_circle_intersection(const Point &p0,
00072 const GraphEdge &edge);
00073
00077 inline double circle_area(const double d12_2, double d01_2, double d02_2){
00078 return 0.5*_effective_R_squared
00079 *acos((d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2)));
00080 }
00081 };
00082
00083
00088 double VAC::edge_circle_intersection(const Point &p0,
00089 const GraphEdge &edge){
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
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
00102 double delta = d12_2*_effective_R_squared - cross*cross;
00103
00104
00105 if (delta<=0){
00106 return circle_area(d12_2, d01_2, d02_2);
00107 }
00108
00109
00110 delta = sqrt(delta);
00111
00112
00113 double b = scalar_product(pdiff, p1);
00114
00115
00116
00117
00118
00119
00120
00121 double tp = (delta-b)/d12_2;
00122
00123
00124 if (tp<0)
00125 return circle_area(d12_2, d01_2, d02_2);
00126
00127
00128 double tm = -(delta+b)/d12_2;
00129
00130
00131 if (tp<1){
00132
00133
00134
00135
00136
00137
00138
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
00144
00145
00146
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
00153
00154
00155 if (tm>1)
00156 return circle_area(d12_2, d01_2, d02_2);
00157
00158
00159 if (tm<0)
00160 return 0.5*fabs(cross);
00161
00162
00163
00164
00165 return (1-tm)*0.5*fabs(cross)
00166 +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
00167 }
00168
00169
00170
00171
00172 VAC::VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &jet_begin,
00173 const vector<PseudoJet>::const_iterator &jet_end,
00174 double effective_R) {
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
00190
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
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
00202
00203
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
00215 maxrap = max(maxrap,rap);
00216 minrap = min(minrap,rap);
00217 }
00218 n_tot++;
00219 }
00220
00221
00222 if (n_added == 0) return;
00223
00224
00225
00226 double max_extend = 2*max(maxrap-minrap+4*_effective_R, twopi+8*_effective_R);
00227 voronoi_particles.push_back(Point(0.5*(minrap+maxrap)-max_extend, pi));
00228 voronoi_particles.push_back(Point(0.5*(minrap+maxrap)+max_extend, pi));
00229 voronoi_particles.push_back(Point(0.5*(minrap+maxrap), pi-max_extend));
00230 voronoi_particles.push_back(Point(0.5*(minrap+maxrap), pi+max_extend));
00231
00232
00233 VoronoiDiagramGenerator vdg;
00234 vdg.generateVoronoi(&voronoi_particles,
00235 0.5*(minrap+maxrap)-max_extend, 0.5*(minrap+maxrap)+max_extend,
00236 pi-max_extend, pi+max_extend);
00237
00238 vdg.resetIterator();
00239 GraphEdge *e=NULL;
00240 unsigned int v_index;
00241 int p_index;
00242 vector<PseudoJet>::const_iterator jet;
00243
00244 while(vdg.getNext(&e)){
00245 v_index = e->point1;
00246 if (v_index<n_added){
00247 p_index = voronoi_indices[v_index];
00248 if (p_index!=-1){
00249 jet = jet_begin+voronoi_indices[v_index];
00250 _areas[p_index]+=
00251 edge_circle_intersection(voronoi_particles[v_index], *e);
00252 }
00253 }
00254 v_index = e->point2;
00255 if (v_index<n_added){
00256 p_index = voronoi_indices[v_index];
00257 if (p_index!=-1){
00258 jet = jet_begin+voronoi_indices[v_index];
00259 _areas[p_index]+=
00260 edge_circle_intersection(voronoi_particles[v_index], *e);
00261 }
00262 }
00263 }
00264
00265
00266 }
00267
00268
00269
00271 void ClusterSequenceVoronoiArea::_initializeVA () {
00272
00273 _pa_calc = new VAC(_jets.begin(),
00274 _jets.begin()+n_particles(),
00275 _effective_Rfact*_jet_def.R());
00276
00277
00278
00279 _voronoi_area.reserve(2*n_particles());
00280 for (unsigned int i=0; i<n_particles(); i++) {
00281 _voronoi_area.push_back(_pa_calc->area(i));
00282
00283 if (_jets[i].perp2() > 0) {
00284 _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp())
00285 * _jets[i]);
00286 } else {
00287
00288
00289 _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
00290 }
00291 }
00292
00293
00294 for (unsigned int i = n_particles(); i < _history.size(); i++) {
00295 double area;
00296 PseudoJet area_4vect;
00297 if (_history[i].parent2 >= 0) {
00298 area = _voronoi_area[_history[i].parent1] +
00299 _voronoi_area[_history[i].parent2];
00300 area_4vect = _voronoi_area_4vector[_history[i].parent1] +
00301 _voronoi_area_4vector[_history[i].parent2];
00302 } else {
00303 area = _voronoi_area[_history[i].parent1];
00304 area_4vect = _voronoi_area_4vector[_history[i].parent1];
00305 }
00306 _voronoi_area.push_back(area);
00307 _voronoi_area_4vector.push_back(area_4vect);
00308 }
00309
00310 }
00311
00312
00313 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
00314 delete _pa_calc;
00315 }
00316
00317 FASTJET_END_NAMESPACE