FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequenceVoronoiArea.cc 2687 2011-11-14 11:17:51Z soyez $ 00003 // 00004 // Copyright (c) 2006-2007 Matteo Cacciari, Gavin Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of a simple command-line handling environment 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 #include "fastjet/ClusterSequenceVoronoiArea.hh" 00030 #include "fastjet/internal/Voronoi.hh" 00031 #include <list> 00032 #include <cassert> 00033 #include <ostream> 00034 #include <fstream> 00035 #include <iterator> 00036 #include <cmath> 00037 #include <limits> 00038 00039 using namespace std; 00040 00041 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00042 00043 typedef ClusterSequenceVoronoiArea::VoronoiAreaCalc VAC; 00044 00045 /// class for carrying out a voronoi area calculation on a set of 00046 /// initial vectors 00047 class ClusterSequenceVoronoiArea::VoronoiAreaCalc { 00048 public: 00049 /// constructor that takes a range of a vector together with the 00050 /// effective radius for the intersection of discs with voronoi 00051 /// cells 00052 VoronoiAreaCalc(const vector<PseudoJet>::const_iterator &, 00053 const vector<PseudoJet>::const_iterator &, 00054 double effective_R); 00055 00056 /// return the area of the particle associated with the given 00057 /// index 00058 inline double area (int index) const {return _areas[index];}; 00059 00060 private: 00061 std::vector<double> _areas; ///< areas, numbered as jets 00062 double _effective_R; ///< effective radius 00063 double _effective_R_squared; ///< effective radius squared 00064 00065 /** 00066 * compute the intersection of one triangle with the circle 00067 * the area is returned 00068 */ 00069 double edge_circle_intersection(const VPoint &p0, 00070 const GraphEdge &edge); 00071 00072 /// get the area of a circle of radius R centred on the point 0 with 00073 /// 1 and 2 on each "side" of the arc. dij is the distance between 00074 /// point i and point j and all distances are squared 00075 inline double circle_area(const double d12_2, double d01_2, double d02_2){ 00076 return 0.5*_effective_R_squared 00077 *acos(min(1.0,(d01_2+d02_2-d12_2)/(2*sqrt(d01_2*d02_2)))); 00078 } 00079 }; 00080 00081 00082 /** 00083 * compute the intersection of one triangle with the circle 00084 * the area is returned 00085 */ 00086 double VAC::edge_circle_intersection(const VPoint &p0, 00087 const GraphEdge &edge){ 00088 VPoint p1(edge.x1-p0.x, edge.y1-p0.y); 00089 VPoint p2(edge.x2-p0.x, edge.y2-p0.y); 00090 VPoint pdiff = p2-p1; 00091 00092 //fprintf(stdout, "\tpt(%f,%f)\n", p0.x, p0.y); 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 // compute intersections between edge line and circle 00100 double delta = d12_2*_effective_R_squared - cross*cross; 00101 00102 // if no intersection, area=area_circle 00103 if (delta<=0){ 00104 return circle_area(d12_2, d01_2, d02_2); 00105 } 00106 00107 // we'll only need delta's sqrt now 00108 delta = sqrt(delta); 00109 00110 // b is the projection of 01 onto 12 00111 double b = scalar_product(pdiff, p1); 00112 00113 // intersections with the circle: 00114 // we compute the "coordinate along the line" of the intersection 00115 // with t=0 (1) corresponding to p1 (p2) 00116 // points with 0<t<1 are within the circle others are outside 00117 00118 // positive intersection 00119 double tp = (delta-b)/d12_2; 00120 00121 // if tp is negative, tm also => inters = circle 00122 if (tp<0) 00123 return circle_area(d12_2, d01_2, d02_2); 00124 00125 // we need the second intersection 00126 double tm = -(delta+b)/d12_2; 00127 00128 // if tp<1, it lies in the circle 00129 if (tp<1){ 00130 // if tm<0, the segment has one intersection 00131 // with the circle at p (t=tp) 00132 // the area is a triangle from 1 to p 00133 // then a circle from p to 2 00134 // several tricks can be used: 00135 // - the area of the triangle is tp*area triangle 00136 // - the lenght for the circle are easily obtained 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 // now, 0 < tm < tp < 1 00142 // the segment intersects twice the circle 00143 // area = 2 cirles at ends + a triangle in the middle 00144 // again, simplifications are staightforward 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 // now, we have tp>1 00151 00152 // if in addition tm>1, intersectino is a circle 00153 if (tm>1) 00154 return circle_area(d12_2, d01_2, d02_2); 00155 00156 // if tm<0, the triangle is inside the circle 00157 if (tm<0) 00158 return 0.5*fabs(cross); 00159 00160 // otherwise, only the "tm point" is on the segment 00161 // area = circle from 1 to m and triangle from m to 2 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 // the constructor... 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<VPoint> 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 // loop over jets and create the triangulation, as well as cross-referencing 00188 // info 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 // generate the corresponding point 00194 double rap = jet_it->rap(), phi = jet_it->phi(); 00195 voronoi_particles.push_back(VPoint(rap, phi)); 00196 voronoi_indices.push_back(n_tot); 00197 n_added++; 00198 00199 // insert a copy of the point if it falls within 2*_R_effective 00200 // of the 0,2pi borders (because we are interested in any 00201 // voronoi edge within _R_effective of the other border) 00202 if (phi < 2*_effective_R) { 00203 voronoi_particles.push_back(VPoint(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(VPoint(rap,phi-twopi)); 00208 voronoi_indices.push_back(-1); 00209 n_added++; 00210 } 00211 00212 // track the rapidity range 00213 maxrap = max(maxrap,rap); 00214 minrap = min(minrap,rap); 00215 } 00216 n_tot++; 00217 } 00218 00219 // allow for 0-particle case in graceful way 00220 if (n_added == 0) return; 00221 // assert(n_added > 0); // old (pre 2.4) non-graceful exit 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(VPoint(0.5*(minrap+maxrap)-max_extend, pi)); 00226 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap)+max_extend, pi)); 00227 voronoi_particles.push_back(VPoint(0.5*(minrap+maxrap), pi-max_extend)); 00228 voronoi_particles.push_back(VPoint(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 } 00265 00266 00267 //---------------------------------------------------------------------- 00268 /// 00269 void ClusterSequenceVoronoiArea::_initializeVA () { 00270 // run the VAC on our original particles 00271 _pa_calc = new VAC(_jets.begin(), 00272 _jets.begin()+n_particles(), 00273 _effective_Rfact*_jet_def.R()); 00274 00275 // transfer the areas to our local structure 00276 // -- first the initial ones 00277 _voronoi_area.reserve(2*n_particles()); 00278 for (unsigned int i=0; i<n_particles(); i++) { 00279 _voronoi_area.push_back(_pa_calc->area(i)); 00280 // make a stab at a 4-vector area 00281 if (_jets[i].perp2() > 0) { 00282 _voronoi_area_4vector.push_back((_pa_calc->area(i)/_jets[i].perp()) 00283 * _jets[i]); 00284 } else { 00285 // not sure what to do here -- just put zero (it won't be meaningful 00286 // anyway) 00287 _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0)); 00288 } 00289 } 00290 00291 // -- then the combined areas that arise from the clustering 00292 for (unsigned int i = n_particles(); i < _history.size(); i++) { 00293 double area_local; 00294 PseudoJet area_4vect; 00295 if (_history[i].parent2 >= 0) { 00296 area_local = _voronoi_area[_history[i].parent1] + 00297 _voronoi_area[_history[i].parent2]; 00298 area_4vect = _voronoi_area_4vector[_history[i].parent1] + 00299 _voronoi_area_4vector[_history[i].parent2]; 00300 } else { 00301 area_local = _voronoi_area[_history[i].parent1]; 00302 area_4vect = _voronoi_area_4vector[_history[i].parent1]; 00303 } 00304 _voronoi_area.push_back(area_local); 00305 _voronoi_area_4vector.push_back(area_4vect); 00306 } 00307 00308 } 00309 00310 //---------------------------------------------------------------------- 00311 ClusterSequenceVoronoiArea::~ClusterSequenceVoronoiArea() { 00312 delete _pa_calc; 00313 } 00314 00315 FASTJET_END_NAMESPACE