JETAN/fastjet/siscone/spherical/momentum.h

   1 // -*- C++ -*-
   2 ///////////////////////////////////////////////////////////////////////////////
   3 // File: momentum.h                                                          //
   4 // Description: header file for 4-momentum class Cmomentum                   //
   5 // This file is part of the SISCone project.                                 //
   6 // WARNING: this is not the main SISCone trunk but                           //
   7 //          an adaptation to spherical coordinates                           //
   8 // For more details, see http://projects.hepforge.org/siscone                //
   9 //                                                                           //
  10 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez                     //
  11 //                                                                           //
  12 // This program is free software; you can redistribute it and/or modify      //
  13 // it under the terms of the GNU General Public License as published by      //
  14 // the Free Software Foundation; either version 2 of the License, or         //
  15 // (at your option) any later version.                                       //
  16 //                                                                           //
  17 // This program is distributed in the hope that it will be useful,           //
  18 // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
  19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
  20 // GNU General Public License for more details.                              //
  21 //                                                                           //
  22 // You should have received a copy of the GNU General Public License         //
  23 // along with this program; if not, write to the Free Software               //
  24 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
  25 //                                                                           //
  26 // $Revision:: 256                                                          $//
  27 // $Date:: 2008-07-14 13:52:16 +0200 (Mon, 14 Jul 2008)                     $//
  28 ///////////////////////////////////////////////////////////////////////////////
  29
  30 #ifndef __SPH_VECTOR_H__
  31 #define __SPH_VECTOR_H__
  32
  33 #include <vector>
  34 #include <math.h>
  35 #include <siscone/reference.h>
  36 #include "geom_2d.h"
  37 #include <siscone/defines.h>
  38
  39 namespace siscone_spherical{
  40
  41 /**
  42  * \class CSph3vector
  43  * \brief base class for managing the spatial part of Cmomentum (defined after)
  44  *
  45  * This class contains the information for particle or group of
  46  * particles management.
  47  * It is adapted to use spherical geometry, where, for our purposes,
  48  * the only time-consuming operation we need is the computation of
  49  * the norm. To compute it once-and-for-all and store it in a local
  50  * variable, you should call the 'build_norm' method.
  51  * On top of that, the angle phi is computed from the x-axis
  52  * and theta from the "north pole".
  53  */
  54 class CSph3vector{
  55  public:
  56   /// default ctor
  57   CSph3vector();
  58
  59   /// ctor with initialisation
  60   CSph3vector(double _px, double _py, double _pz);
  61
  62   /// default dtor
  63   ~CSph3vector();
  64
  65   /// assignment of vectors
  66   CSph3vector& operator = (const CSph3vector &v);
  67
  68   /// addition of vectors
  69   /// WARNING= norm is not updated
  70   const CSph3vector operator + (const CSph3vector &v);
  71
  72   /// subtraction of vectors
  73   /// WARNING= norm is not updated
  74   const CSph3vector operator - (const CSph3vector &v);
  75
  76   /// division by a constant
  77   /// WARNING= norm is not updated
  78   const CSph3vector operator / (const double &r);
  79
  80   /// incrementation of vectors
  81   /// WARNING= norm is not updated
  82   CSph3vector& operator += (const CSph3vector &v);
  83
  84   /// decrementation of vectors
  85   /// WARNING= norm is not updated
  86   CSph3vector& operator -= (const CSph3vector &v);
  87
  88   /// multiplication by a constant
  89   /// WARNING= norm is not updated
  90   CSph3vector& operator *= (const double &r);
  91
  92   /// division by a constant
  93   /// WARNING= norm is not updated
  94   CSph3vector& operator /= (const double &r);
  95
  96   /// computes pT
  97   inline double perp() const {return sqrt(perp2());}
  98
  99   /// computes pT^2
 100   inline double perp2() const {return px*px+py*py;}
 101
 102   /// 3-vect norm
 103   inline double norm() const {return sqrt(px*px+py*py+pz*pz);}
 104
 105   /// 3-vect norm squared
 106   inline double norm2() const {return px*px+py*py+pz*pz;}
 107
 108   /// 3-vect azimuthal angle
 109   inline double phi() const {return atan2(py, px);}
 110
 111   /// 3-vect polar angle
 112   inline double theta() const {return atan2(perp(),pz);}
 113
 114   /// build the spatial normfrom 4-momentum info
 115   /// !!!                  WARNING                       !!!
 116   /// !!! computing the norm is the only time-consuming  !!!
 117   /// !!! information we need in all computations.       !!!
 118   /// !!! use this whenever you need repeated access     !!!
 119   /// !!! to the norm to store it in the local variable  !!!
 120   void build_norm();
 121
 122   /// just a useful tool to store theta and phi
 123   /// locally (in _theta and _phi) in case you need
 124   /// repeated access
 125   void build_thetaphi();
 126
 127   /// for this direction, compute the two reference directions
 128   /// used to measure angles
 129   void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2);
 130
 131   double px;        ///< x-momentum
 132   double py;        ///< y-momentum
 133   double pz;        ///< z-momentum
 134
 135   double _norm;     ///< particle spatial norm (available ONLY after a call to build_norm)
 136   double _theta;    ///< particle theta angle (available ONLY after a call to build_thetaphi)
 137   double _phi;      ///< particle phi angle   (available ONLY after a call to build_thetaphi)
 138
 139   //////////////////////////////////////////////
 140   // the following part is used for checksums //
 141   //////////////////////////////////////////////
 142   siscone::Creference ref;   ///< reference number for the vector
 143 };
 144
 145 /**
 146  * \class CSphmomentum
 147  * \brief base class for dynamic coordinates management
 148  *
 149  * This class contains the information for particle or group of
 150  * particles management.
 151  * It is adapted to use spherical geometry, where, for our purposes,
 152  * the only time-consuming operation we need is the computation of
 153  * the norm. To compute it once-and-for-all and store it in a local
 154  * variable, you should call the 'build_norm' method.
 155  * On top of that, the angle phi is computed from the x-axis
 156  * and theta from the "north pole".
 157  */
 158 class CSphmomentum : public CSph3vector{
 159  public:
 160   /// default ctor
 161   CSphmomentum();
 162
 163   /// init from a 3-vect
 164   CSphmomentum(CSph3vector &init, double E=0.0);
 165
 166   /// ctor with initialisation
 167   CSphmomentum(double _px, double _py, double _pz, double _E);
 168
 169   /// ctor with detailed initialisation
 170   //CSphmomentum(double _eta, double _phi, siscone::Creference _ref);
 171
 172   /// default dtor
 173   ~CSphmomentum();
 174
 175   /// computes m
 176   inline double mass() const {return sqrt(mass2());}
 177
 178   /// computes m^2
 179   inline double mass2() const {return perpmass2()-perp2();}
 180
 181   /// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2)
 182   inline double perpmass() const {return sqrt((E-pz)*(E+pz));}
 183
 184   /// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2
 185   inline double perpmass2() const {return (E-pz)*(E+pz);}
 186
 187   /// computes transverse energy
 188   inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());}
 189
 190   /// computes transverse energy (squared)
 191   inline double Et2() const {return E*E/(1.0+pz*pz/perp2());}
 192
 193   /// assignment of vectors
 194   CSphmomentum& operator = (const CSphmomentum &v);
 195
 196   /// addition of vectors
 197   /// !!! WARNING !!! no updating of eta and phi !!!
 198   const CSphmomentum operator + (const CSphmomentum &v);
 199
 200   /// incrementation of vectors
 201   /// !!! WARNING !!! no updating of eta and phi !!!
 202   CSphmomentum& operator += (const CSphmomentum &v);
 203
 204   /// decrementation of vectors
 205   /// !!! WARNING !!! no updating of eta and phi !!!
 206   CSphmomentum& operator -= (const CSphmomentum &v);
 207
 208   double E;         ///< energy
 209
 210   int parent_index; ///< particle number in the parent list
 211   int index;        ///< internal particle number
 212 };
 213
 214 /// ordering of two vectors
 215 /// this is by default done w.r.t. their references
 216 bool operator < (const CSphmomentum &v1, const CSphmomentum &v2);
 217
 218 /// ordering of vectors in eta (e.g. used in collinear tests)
 219 bool momentum_theta_less(const CSphmomentum &v1, const CSphmomentum &v2);
 220
 221 /// ordering of vectors in pt
 222 bool momentum_pt_less(const CSphmomentum &v1, const CSphmomentum &v2);
 223
 224
 225 //////////////////////////
 226 // some handy utilities //
 227 //////////////////////////
 228
 229 /// square
 230 inline double sqr(double x){return x*x;}
 231
 232 /// dot product for te spatial 3-vect
 233 /// \param v1    first 4-vect
 234 /// \param v2    second 4-vect
 235 inline double dot_product3(const CSph3vector &v1, const CSph3vector &v2){
 236   //double tmp = v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
 237   //if (!isfinite(tmp)){
 238   //  std::cout << "dot_product inf: " << std::endl;
 239   //  std::cout << "  angles: " << v1._theta << " " << v1._phi << " and " << v2._theta << " " << v2._phi << std::endl;
 240   //  std::cout << "  moms  : " << v1.px << " " << v1.py << " " << v1.pz
 241   //          << " and "      << v2.px << " " << v2.py << " " << v2.pz << std::endl;
 242   //}
 243   return v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
 244 }
 245
 246 /// cross product for the spatial 3-vect
 247 /// \param v1    first 4-vect
 248 /// \param v2    second 4-vect
 249 inline CSph3vector cross_product3(const CSph3vector &v1, const CSph3vector &v2){
 250   //CSph3vector tmp;
 251   //tmp.px = v1.py*v2.pz-v1.pz*v2.py;
 252   //tmp.py = v1.pz*v2.px-v1.px*v2.pz;
 253   //tmp.pz = v1.px*v2.py-v1.py*v2.px;
 254   //return tmp;
 255   return CSph3vector(v1.py*v2.pz-v1.pz*v2.py,
 256                   v1.pz*v2.px-v1.px*v2.pz,
 257                   v1.px*v2.py-v1.py*v2.px);
 258 }
 259
 260 /// squared norm of the cross product for the spatial 3-vect (energy is set to 0)
 261 /// \param v1    first 4-vect
 262 /// \param v2    second 4-vect
 263 inline double norm2_cross_product3(const CSph3vector &v1, const CSph3vector &v2){
 264   return sqr(v1.py*v2.pz-v1.pz*v2.py) + sqr(v1.pz*v2.px-v1.px*v2.pz) + sqr(v1.px*v2.py-v1.py*v2.px);
 265 }
 266
 267 /// get tangent squared of the spherical distance between two vectors
 268 /// \param v1    vector defining the first point
 269 /// \param v2    vector defining the second point
 270 inline double get_tan2_distance(const CSphmomentum &v1, const CSphmomentum &v2){
 271   return norm2_cross_product3(v1,v2)/sqr(dot_product3(v1,v2));
 272 }
 273
 274 /// get spherical distance between to vectors
 275 /// \param v1    vector defining the first point
 276 /// \param v2    vector defining the second point
 277 inline double get_distance(const CSph3vector *v1, const CSph3vector *v2){
 278   return atan2(sqrt(norm2_cross_product3(*v1,*v2)), dot_product3(*v1,*v2));
 279 }
 280
 281 /// return true if the two points are distant by less than get spherical distance between two vectors
 282 /// \param v1      vector defining the first point
 283 /// \param v2      vector defining the second point
 284 /// \param tan2R   tangent squared of the max distance
 285 /// WARNING: using the tangent here is dangerous for R>pi/2.
 286 ///          this never happens per se for "regular R" but
 287 ///          it may in the vicinity computation as we're using
 288 ///          2R there.
 289 inline bool is_closer(const CSph3vector *v1, const CSph3vector *v2, const double tan2R){
 290   double dot = dot_product3(*v1,*v2);
 291   return (dot>=0) && (norm2_cross_product3(*v1,*v2)<=tan2R*dot*dot);
 292 }
 293
 294 /// return true if the two points are distant by less than  get spherical distance between to vectors
 295 /// \param v1      vector defining the first point
 296 /// \param v2      vector defining the second point
 297 /// \param tan2R   tangent squared of the max distance
 298 /// safer version but computes the norm
 299 inline bool is_closer_safer(const CSph3vector *v1, const CSph3vector *v2, const double cosR){
 300   return dot_product3(*v1,*v2)>=cosR*sqrt(v1->norm2()*v2->norm2());
 301   //double dot = dot_product3(*v1,*v2);
 302   //return (dot>=0) && (norm2_cross_product3(*v1,*v2)<tan2R*dot*dot);
 303 }
 304
 305 /// multiply a vector by a constant
 306 /// WARNING: norm not updated
 307 inline CSph3vector operator * (const double &r, const CSph3vector &v){
 308   CSph3vector tmp = v;
 309   return tmp*=r;
 310 }
 311 }
 312 #endif