ClusterSequenceVoronoiArea.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequenceVoronoiArea.cc 644 2007-05-16 21:45:21Z 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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
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      // defined in fastjet/internal/base.hh
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   //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<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   // 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(Point(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(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       // track the rapidity range
00213       maxrap = max(maxrap,rap);
00214       minrap = min(minrap,rap);
00215     }
00216     n_tot++;
00217   }
00218 
00219   assert(n_added > 0);
00220 
00221   // add extreme cases:
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   // Build the VD
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   // run the VAC on our original particles
00268   _pa_calc = new VAC(_jets.begin(), 
00269                      _jets.begin()+n_particles(),
00270                      _effective_Rfact*_jet_def.R());
00271 
00272   // transfer the areas to our local structure
00273   //  -- first the initial ones
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     // make a stab at a 4-vector area
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       // not sure what to do here -- just put zero (it won't be meaningful
00283       // anyway)
00284       _voronoi_area_4vector.push_back(PseudoJet(0.0,0.0,0.0,0.0));
00285     }
00286   }
00287            
00288   //  -- then the combined areas that arise from the clustering
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

Generated on Thu Jan 3 19:04:30 2008 for fastjet by  doxygen 1.5.2