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 <iterator>
00037 #include <cmath>
00038
00039 using namespace std;
00040
00041 FASTJET_BEGIN_NAMESPACE
00042
00043 typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC;
00044
00047 class ClusterSequenceVoronoiArea::VoronoiAreaCalc {
00048 public:
00052 VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &,
00053 const vector<PseudoJet>::const_iterator &,
00054 double effective_R);
00055
00058 inline double area (int index) const {return _areas[index];};
00059
00060 private:
00061 std::vector<double> _areas;
00062 double _effective_R;
00063 double _effective_R_squared;
00064
00069 double edge_circle_intersection(const Point &p0,
00070 const GraphEdge &edge);
00071
00075 inline double circle_area(const double d12_2, double d01_2, double d02_2){
00076 return 0.5*_effective_R_squared
00077 *acos((d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2)));
00078 }
00079 };
00080
00081
00086 double VAC::edge_circle_intersection(const Point &p0,
00087 const GraphEdge &edge){
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
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
00100 double delta = d12_2*_effective_R_squared - cross*cross;
00101
00102
00103 if (delta<=0){
00104 return circle_area(d12_2, d01_2, d02_2);
00105 }
00106
00107
00108 delta = sqrt(delta);
00109
00110
00111 double b = scalar_product(pdiff, p1);
00112
00113
00114
00115
00116
00117
00118
00119 double tp = (delta-b)/d12_2;
00120
00121
00122 if (tp<0)
00123 return circle_area(d12_2, d01_2, d02_2);
00124
00125
00126 double tm = -(delta+b)/d12_2;
00127
00128
00129 if (tp<1){
00130
00131
00132
00133
00134
00135
00136
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
00142
00143
00144
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
00151
00152
00153 if (tm>1)
00154 return circle_area(d12_2, d01_2, d02_2);
00155
00156
00157 if (tm<0)
00158 return 0.5*fabs(cross);
00159
00160
00161
00162
00163 return (1-tm)*0.5*fabs(cross)
00164 +circle_area(tm*tm*d12_2, d01_2, _effective_R_squared);
00165 }
00166
00167
00168
00169
00170 VAC::VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &jet_begin,
00171 const vector<PseudoJet>::const_iterator &jet_end,
00172 double effective_R) {
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
00188
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
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
00200
00201
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
00213 maxrap = max(maxrap,rap);
00214 minrap = min(minrap,rap);
00215 }
00216 n_tot++;
00217 }
00218
00219 assert(n_added > 0);
00220
00221
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
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 }
00262
00263
00264
00266 void ClusterSequenceVoronoiArea::_initializeVA () {
00267
00268 _pa_calc = new VAC(_jets.begin(),
00269 _jets.begin()+n_particles(),
00270 _effective_Rfact*_jet_def.R());
00271
00272
00273
00274 _voronoi_area.reserve(2*n_particles());
00275 for (unsigned int i=0; i<n_particles(); i++) {
00276 _voronoi_area.push_back(_pa_calc->area(i));
00277
00278 if (_jets[i].perp2() > 0) {
00279 _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp())
00280 * _jets[i]);
00281 } else {
00282
00283
00284 _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
00285 }
00286 }
00287
00288
00289 for (unsigned int i = n_particles(); i < _history.size(); i++) {
00290 double area;
00291 PseudoJet area_4vect;
00292 if (_history[i].parent2 >= 0) {
00293 area = _voronoi_area[_history[i].parent1] +
00294 _voronoi_area[_history[i].parent2];
00295 area_4vect = _voronoi_area_4vector[_history[i].parent1] +
00296 _voronoi_area_4vector[_history[i].parent2];
00297 } else {
00298 area = _voronoi_area[_history[i].parent1];
00299 area_4vect = _voronoi_area_4vector[_history[i].parent1];
00300 }
00301 _voronoi_area.push_back(area);
00302 _voronoi_area_4vector.push_back(area_4vect);
00303 }
00304
00305 }
00306
00307
00308 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() {
00309 delete _pa_calc;
00310 }
00311
00312 FASTJET_END_NAMESPACE