include/DDRec/Vector3D.h

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : F.Gaede
0011 //
0012 //==========================================================================
0013 #ifndef DDREC_VECTOR3D_H
0014 #define DDREC_VECTOR3D_H 1
0015
0016 #include <cmath>
0017 #include <iostream>
0018 #include <cassert>
0019
0020
0021 namespace dd4hep{ namespace rec {
0022
0023   /** Simple three dimensional vector providing the components for cartesian,
0024    *  cylindrical and spherical coordinate  systems - internal reperesentation is
0025    *  cartesian. (copy of original version from gear).
0026    *
0027    *  @author F. Gaede, DESY
0028    *  @version $Id$
0029    *  @date Apr 6 2014
0030    */
0031
0032   class Vector3D{
0033
0034   public:
0035
0036     /** Default c'tor - zero vector */
0037     Vector3D() : _x(0.0),_y(0.0),_z(0.0) {}
0038
0039
0040     /** Copy constructor*/
0041     Vector3D(const Vector3D& v) : _x(v[0]),_y(v[1]),_z(v[2]) {}
0042
0043     /** Constructor for float array.*/
0044     Vector3D(const float* v) : _x(v[0]),_y(v[1]),_z(v[2]) {}
0045
0046     /** Constructor for double array.*/
0047     Vector3D(const double* v) : _x(v[0]),_y(v[1]),_z(v[2]) {}
0048
0049
0050     /** Templated c'tor - allows to have overloaded c'tors for different coordinates */
0051     template <class T>
0052     Vector3D( double x,double y, double z , T(&)() ) ;
0053
0054
0055     /** Default corrdinate system for initialization is cartesian */
0056     Vector3D( double x_val,double y_val, double z_val ) :
0057       _x(x_val),
0058       _y(y_val),
0059       _z(z_val) {
0060     }
0061
0062     // ---- this causes all sorts of template lookup errors ...
0063     // /** Copy c'tor for three vectors from other packages - requires T::x(),T::y(), T::z().
0064     //  */
0065     // template <class T>
0066     // Vector3D( const T& t) :
0067     //   _x( t.x() ) ,
0068     //   _y( t.y() ) ,
0069     //   _z( t.z() ){
0070     // }
0071
0072     //assignment operator
0073     Vector3D& operator=(const Vector3D& v) {
0074       _x = v[0] ;
0075       _y = v[1] ;
0076       _z = v[2] ;
0077       return *this ;
0078     }
0079
0080     /// fill vector from arbitrary class that defines operator[]
0081     template <class T>
0082     inline const Vector3D& fill(  const T& v  ) {
0083
0084       _x = v[0] ; _y = v[1] ; _z = v[2] ;
0085       return *this ;
0086     }
0087
0088     /// fill vector from double array
0089     inline const Vector3D& fill( const double* v) {
0090
0091       _x = v[0] ; _y = v[1] ; _z = v[2] ;
0092       return *this ;
0093     }
0094
0095     /// fill from double values
0096     inline const Vector3D& fill( double x_val, double y_val, double z_val) {
0097       _x = x_val ; _y = y_val ; _z = z_val ;
0098       return *this ;
0099     }
0100
0101
0102     /** Cartesian x coordinate */
0103     inline double x() const { return  _x ; }
0104
0105     /** Cartesian y coordinate */
0106     inline double y() const { return  _y ; }
0107
0108     /** Cartesian cartesian  z coordinate */
0109     inline double z() const { return  _z ; }
0110
0111     /** Assign to  cartesian x coordinate */
0112     inline double& x() { return  _x ; }
0113
0114     /**  Assign to  cartesian y coordinate */
0115     inline double& y() { return  _y ; }
0116
0117     /**  Assign to cartesian z coordinate */
0118     inline double& z() { return  _z ; }
0119
0120
0121    /** Accessing x,y,z with bracket operator */
0122     inline double operator[](int i) const {
0123       switch(i) {
0124       case 0: return _x ; break ;
0125       case 1: return _y ; break ;
0126       case 2: return _z ; break ;
0127       }
0128       return 0.0 ;
0129     }
0130     /** Accessing x,y,z with bracket operator for assignment */
0131     inline double& operator[](int i) {
0132       switch(i) {
0133       case 0: return _x ; break ;
0134       case 1: return _y ; break ;
0135       case 2: return _z ; break ;
0136       }
0137       static double dummy(0.0) ;
0138       return dummy ;
0139     }
0140
0141     /** Azimuthal angle - cylindrical and spherical */
0142     inline double phi() const {
0143
0144       return _x == 0.0 && _y == 0.0 ? 0.0 : atan2(_y,_x);
0145     }
0146
0147     /** Transversal component - cylindrical 'r' */
0148     inline double rho() const {
0149
0150       return trans() ;
0151     }
0152
0153     /** Transversal component */
0154     inline double trans() const {
0155
0156       return sqrt( _x*_x + _y*_y ) ;
0157     }
0158
0159     /** Transversal component squared */
0160     inline double trans2() const {
0161
0162       return  _x*_x + _y*_y  ;
0163     }
0164
0165     /** Spherical r/magnitude */
0166     inline double r() const {
0167
0168       return sqrt( _x*_x + _y*_y + _z*_z ) ;
0169     }
0170
0171
0172     /** Spherical r/magnitude, squared */
0173     inline double r2() const {
0174
0175       return  _x*_x + _y*_y + _z*_z  ;
0176     }
0177
0178     /** Polar angle - spherical */
0179     inline double theta() const {
0180
0181       return _x == 0.0 && _y == 0.0 && _z == 0.0 ? 0.0 : atan2( rho(),_z) ;
0182     }
0183
0184     /** Scalar product */
0185     inline double dot( const Vector3D& v) const {
0186       return _x * v.x() + _y * v.y() + _z * v.z() ;
0187     }
0188
0189
0190     /** Vector product */
0191     inline Vector3D cross( const Vector3D& v) const {
0192
0193       return Vector3D( _y * v.z() - _z * v.y() ,
0194                _z * v.x() - _x * v.z() ,
0195                _x * v.y() - _y * v.x() )  ;
0196     }
0197
0198     /** Parallel unit vector */
0199     inline Vector3D unit() const {
0200
0201       double n = r() ;
0202       return Vector3D( _x / n , _y / n , _z / n ) ;
0203     }
0204
0205
0206     /// direct access to data as const double*
0207     inline operator const double*() const {
0208       return &_x ;
0209     }
0210     /// direct access to data as const double*
0211     inline const double* const_array() const {
0212       return &_x ;
0213     }
0214
0215     /// direct access to data as double* - allows modification
0216     inline double* array() {
0217       return &_x ;
0218     }
0219
0220
0221       /** Component wise comparison of two vectors - true if all components differ less than epsilon */
0222       inline bool isEqual(  const Vector3D& b , double epsilon=1e-6) {
0223
0224       if( fabs( x() - b.x() ) < epsilon &&
0225       fabs( y() - b.y() ) < epsilon &&
0226       fabs( z() - b.z() ) < epsilon )
0227     return true;
0228       else
0229     return false;
0230   }
0231
0232
0233
0234     // this causes template lookup errors on some machines :
0235     //   -> use explicit conversion with to<T>()
0236     // /** Implicit templated conversion to anything that has a c'tor T(x,y,z)
0237     //  *  and accessor functions x(),y(),z(). For safety the result is checked which
0238     //  *  causes a small performance penalty.
0239     //  *  @see to()
0240     //  *
0241     //  */
0242     // template <class T>
0243     // inline operator T() const {
0244
0245     //   T t( _x, _y , _z ) ;
0246
0247     //   assert( t.x()== _x && t.y()== _y && t.z()== _z ) ;
0248
0249     //   return t ;
0250
0251     //   //     return T( _x, _y, _z ) ;
0252     // }
0253
0254
0255     /** Explicit, unchecked conversion to anything that has a c'tor T(x,y,z).
0256      *  Example:
0257      *  CLHEP::Vector3D clhv = v.to< CLHEP::Vector3D>() ;
0258      *  @see operator T()
0259      */
0260     template <class T>
0261     inline T to() const { return T( _x, _y, _z ) ; }
0262
0263
0264   protected:
0265
0266     double _x,_y,_z ;
0267
0268
0269     // helper classes and function to allow
0270     // different c'tors selected at compile time
0271   public:
0272
0273     struct Cartesian   { } ;
0274     struct Cylindrical { } ;
0275     struct Spherical   { } ;
0276
0277     static Cartesian   cartesian()  { return Cartesian() ; }
0278     static Cylindrical cylindrical(){ return Cylindrical() ;}
0279     static Spherical   spherical()  { return Spherical() ; }
0280
0281   } ;
0282
0283   /** Addition of two vectors */
0284   inline Vector3D operator+(  const Vector3D& a, const Vector3D& b ) {
0285
0286     return Vector3D( a.x() + b.x()  , a.y() + b.y(), a.z() + b.z()  ) ;
0287   }
0288   /** Subtraction of two vectors */
0289   inline Vector3D operator-(  const Vector3D& a, const Vector3D& b ) {
0290
0291     return Vector3D( a.x() - b.x()  , a.y() - b.y(), a.z() - b.z()  ) ;
0292   }
0293   /** Exact comparison of two vectors */
0294   inline bool operator==(  const Vector3D& a, const Vector3D& b ) {
0295
0296     if( a.x() == b.x()  &&  a.y() == b.y() && a.z() == b.z()  )
0297       return true;
0298     else
0299       return false;
0300   }
0301
0302   /** Multiplication with scalar */
0303   inline Vector3D operator*( double s , const Vector3D& v ) {
0304
0305     return Vector3D( s * v.x()  , s * v.y()  ,  s * v.z() ) ;
0306   }
0307
0308  /** Negative vector */
0309   inline Vector3D operator-(  const Vector3D& v) {
0310
0311     return Vector3D( -v.x(), - v.y(), - v.z()  ) ;
0312   }
0313
0314   /// operator for scalar product
0315   inline double operator*( const Vector3D& v0, const Vector3D& v1 ){
0316     return v0.dot( v1 ) ;
0317   }
0318
0319
0320   // template specializations for constructors of  different coordinate systems
0321
0322   /** Cartesian c'tor  - example: <br>
0323    *  Vector3D  v( x, y, c , Vector3D::cartesian ) ;
0324    */
0325   template <>
0326   inline Vector3D::Vector3D( double x_val,double y_val, double z_val, Vector3D::Cartesian (&)() ) :
0327     _x(x_val),
0328     _y(y_val),
0329     _z(z_val) {
0330   }
0331
0332   /** Cylindrical c'tor  - example: <br>
0333    *  Vector3D  v( rho, phi, z , Vector3D::cylindrical ) ;
0334    */
0335   template <>
0336   inline Vector3D::Vector3D( double rho_val,double phi_val, double z_val, Vector3D::Cylindrical (&)() ) : _z(z_val)  {
0337
0338     _x = rho_val * cos( phi_val ) ;
0339     _y = rho_val * sin( phi_val ) ;
0340   }
0341
0342
0343   /** Spherical c'tor  - example: <br>
0344    *  Vector3D  v( r, phi, theta , Vector3D::spherical ) ;
0345    */
0346   template <>
0347   inline Vector3D::Vector3D( double r_val,double phi_val, double theta_val, Vector3D::Spherical (&)() ) {
0348     double rst =  r_val * sin( theta_val ) ;
0349     _x = rst * cos( phi_val ) ;
0350     _y = rst * sin( phi_val ) ;
0351     _z = r_val * cos( theta_val ) ;
0352   }
0353
0354
0355
0356   /** Output operator */
0357   inline std::ostream & operator << (std::ostream & os, const Vector3D &v) {
0358
0359     //    os << "( " << v[0] << ", " << v[1] << ", " << v[2] << " )" ;
0360     os << "  ( " << v[0]
0361        << ", " << v[1]
0362        << ", " << v[2]
0363        << " ) -  [ phi: " << v.phi()
0364        << " , rho: " << v.rho() << " ] "
0365        << "  [ theta: " << v.theta()
0366        << " , r: " << v.r() << " ] " ;
0367
0368     return os ;
0369   }
0370
0371
0372
0373 }} // namespace
0374
0375
0376
0377
0378
0379
0380 #endif