sseBasic3DVector.h

Go to the documentation of this file.

#ifndef GeometryVector_newBasic3DVector_h
#define GeometryVector_newBasic3DVector_h

#include "DataFormats/GeometryVector/interface/Basic2DVector.h"
#include "DataFormats/GeometryVector/interface/Theta.h"
#include "DataFormats/GeometryVector/interface/Phi.h"
#include "DataFormats/GeometryVector/interface/PreciseFloatType.h"
#include "DataFormats/GeometryVector/interface/CoordinateSets.h"
#include "DataFormats/Math/interface/SSEVec.h"
#include <iosfwd>
#include <cmath>

namespace detailsBasic3DVector {
  inline float __attribute__((always_inline)) __attribute__ ((pure))
  eta(float x, float y, float z) { float t(z/std::sqrt(x*x+y*y)); return ::asinhf(t);} 
  inline double __attribute__((always_inline)) __attribute__ ((pure))
  eta(double x, double y, double z) { double t(z/std::sqrt(x*x+y*y)); return ::asinh(t);} 
  inline long double __attribute__((always_inline)) __attribute__ ((pure))
  eta(long double x, long double y, long double z) { long double t(z/std::sqrt(x*x+y*y)); return ::asinhl(t);} 
}


template < typename T> 
class Basic3DVector {
public:

  typedef T                                   ScalarType;
  typedef mathSSE::Vec4<T>                    VectorType;
  typedef mathSSE::Vec4<T>                    MathVector;
  typedef Geom::Cylindrical2Cartesian<T>      Cylindrical;
  typedef Geom::Spherical2Cartesian<T>        Spherical;
  typedef Spherical                           Polar; // synonym
    
  Basic3DVector() {}

  Basic3DVector( const Basic3DVector & p) : 
    v(p.v) {}

  template <class U>
  Basic3DVector( const Basic3DVector<U> & p) : 
    v(p.v) {}


  Basic3DVector( const Basic2DVector<T> & p) : 
    v(p.x(),p.y(),0) {}

 
  template <class OtherPoint> 
  explicit Basic3DVector( const OtherPoint& p) : 
        v(p.x(),p.y(),p.z()) {}


  // constructor from Vec4
  template<class U>
  Basic3DVector(mathSSE::Vec4<U> const& iv) : v(iv){}

  Basic3DVector( const T& x, const T& y, const T& z, const T&w=0) : 
    v(x,y,z,w){}

  template <typename U>
  Basic3DVector( const Geom::Theta<U>& theta, 
         const Geom::Phi<U>& phi, const T& r) {
    Polar p( theta.value(), phi.value(), r);
    v.o.theX = p.x(); v.o.theY = p.y(); v.o.theZ = p.z();
  }

  MathVector const & mathVector() const { return v;}
  MathVector & mathVector() { return v;}

  T operator[](int i) const { return v[i];}
  T & operator[](int i) { return v[i];}

  T x() const { return v.o.theX;}

  T y() const { return v.o.theY;}

  T z() const { return v.o.theZ;}

  T w() const { return v.o.theW;}

  Basic2DVector<T> xy() const { return v.xy();}

  // equality
  bool operator==(const Basic3DVector& rh) const {
    return v==rh.v;
  }

  T mag2() const { return  ::dot(v,v);}

  T mag() const  { return std::sqrt( mag2());}

  T perp2() const { return ::dotxy(v,v);}

  T perp() const { return std::sqrt( perp2());}

  T transverse() const { return perp();}

  T barePhi() const {return std::atan2(y(),x());}
  Geom::Phi<T> phi() const {return Geom::Phi<T>(barePhi());}

  T bareTheta() const {return std::atan2(perp(),z());}
  Geom::Theta<T> theta() const {return Geom::Theta<T>(std::atan2(perp(),z()));}

  // T eta() const { return -log( tan( theta()/2.));} 
  T eta() const { return detailsBasic3DVector::eta(x(),y(),z());} // correct 

  Basic3DVector unit() const {
    T my_mag = mag2();
    return (0!=my_mag) ? (*this)*(T(1)/std::sqrt(my_mag)) : *this;
  }

  template <class U> 
  Basic3DVector& operator+= ( const Basic3DVector<U>& p) {
    v = v + p.v;
    return *this;
  } 

  template <class U> 
  Basic3DVector& operator-= ( const Basic3DVector<U>& p) {
    v = v - p.v;
    return *this;
  } 

  Basic3DVector operator-() const { return Basic3DVector(-v);}

  Basic3DVector& operator*= ( T t) {
    v = t*v;
    return *this;
  } 

  Basic3DVector& operator/= ( T t) {
    //t = T(1)/t;
    v = v/t;
    return *this;
  } 

  T dot( const Basic3DVector& rh) const { 
    return ::dot(v,rh.v);
  }

  template <class U> 
  typename PreciseFloatType<T,U>::Type dot( const Basic3DVector<U>& lh) const { 
    return Basic3DVector<typename PreciseFloatType<T,U>::Type>(*this)
      .dot(Basic3DVector<typename PreciseFloatType<T,U>::Type>(lh));
  }

  Basic3DVector cross( const Basic3DVector& lh) const {
    return ::cross(v,lh.v);
  }


  template <class U> 
  Basic3DVector<typename PreciseFloatType<T,U>::Type> 
  cross( const Basic3DVector<U>& lh) const {
    return Basic3DVector<typename PreciseFloatType<T,U>::Type>(*this)
      .cross(Basic3DVector<typename PreciseFloatType<T,U>::Type>(lh));
  }

public:
  mathSSE::Vec4<T> v;
}  __attribute__ ((aligned (16)));


namespace geometryDetails {
  std::ostream & print3D(std::ostream& s, double x, double y, double z);
}

template <class T>
inline std::ostream & operator<<( std::ostream& s, const Basic3DVector<T>& v) {
  return geometryDetails::print3D(s, v.x(),v.y(), v.z());
}


template <class T>
inline Basic3DVector<T>
operator+( const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
  return a.v+b.v;
}
template <class T>
inline Basic3DVector<T>
operator-( const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
  return a.v-b.v;
}

template <class T, class U>
inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
operator+( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
  typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
  return RT(a).v+RT(b).v;
}

template <class T, class U>
inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
operator-( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
  typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
  return RT(a).v-RT(b).v;
}

template <class T>
inline T operator*( const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
  return v1.dot(v2);
}

template <class T, class U>
inline typename PreciseFloatType<T,U>::Type operator*( const Basic3DVector<T>& v1, 
                               const Basic3DVector<U>& v2) {
  return  v1.dot(v2);
}

template <class T>
inline Basic3DVector<T> operator*( const Basic3DVector<T>& v, T t) {
  return v.v*t;
}

template <class T>
inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
  return v.v*t;
}



template <class T, typename S>
inline Basic3DVector<T> operator*(S t,  const Basic3DVector<T>& v) {
  return static_cast<T>(t)*v;
}

template <class T, typename S>
inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
  return static_cast<T>(t)*v;
}


template <class T>
inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, T t) {
  return v.v/t;
}

template <class T, typename S>
inline Basic3DVector<T> operator/( const Basic3DVector<T>& v, S s) {
  //  T t = S(1)/s; return v*t;
  T t = s;
  return v/t;
}


typedef Basic3DVector<float> Basic3DVectorF;
typedef Basic3DVector<double> Basic3DVectorD;


//  add long double specialization
#include "Basic3DVectorLD.h"

#endif // GeometryVector_Basic3DVector_h