extBasic3DVector.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/ExtVec.h"
#include "FWCore/Utilities/interface/Likely.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);
  }
}  // namespace detailsBasic3DVector

template <typename T>
class Basic3DVector {
public:
  typedef T ScalarType;
  typedef Vec4<T> VectorType;
  typedef Vec4<T> MathVector;
  typedef Geom::Cylindrical2Cartesian<T> Cylindrical;
  typedef Geom::Spherical2Cartesian<T> Spherical;
  typedef Spherical Polar;  // synonym

  Basic3DVector() : v{0, 0, 0, 0} {}

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

  Basic3DVector& operator=(const Basic3DVector&) = default;

  template <class U>
  Basic3DVector(const Basic3DVector<U>& p) : v{T(p.v[0]), T(p.v[1]), T(p.v[2]), T(p.v[3])} {}

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

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

  // constructor from Vec4
  Basic3DVector(MathVector const& iv) : v(iv) {}

  template <class U>
  Basic3DVector(Vec4<U> const& iv) : v{T(iv[0]), T(iv[1]), T(iv[2]), T(iv[3])} {}

  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[0] = p.x();
    v[1] = p.y();
    v[2] = p.z();
  }

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

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

  T x() const { return v[0]; }

  T y() const { return v[1]; }

  T z() const { return v[2]; }

  T w() const { return v[3]; }

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

  // equality
  bool operator==(const Basic3DVector& rh) const {
    auto res = v == rh.v;
    return res[0] & res[1] & res[2] & res[3];
  }

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

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

  T perp2() const { return ::dot2(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 LIKELY(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 ::cross3(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:
  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