FastJet 3.0.0
ClusterSequenceVoronoiArea.cc
00001 //STARTHEADER
00002 // $Id: ClusterSequenceVoronoiArea.cc 2577 2011-09-13 15:11:38Z salam $
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;
00294     PseudoJet area_4vect;
00295     if (_history[i].parent2 >= 0) {
00296       area = _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 = _voronoi_area[_history[i].parent1];
00302       area_4vect = _voronoi_area_4vector[_history[i].parent1];
00303     }
00304     _voronoi_area.push_back(area);
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends