BasicVector3D.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BasicVector3D.h

// -*- C++ -*-
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
//
// History:
// 12.06.01 E.Chernyaev - CLHEP-1.7: initial  version
// 14.03.03 E.Chernyaev - CLHEP-1.9: template version
//

#ifndef BASIC_VECTOR3D_H
#define BASIC_VECTOR3D_H

#include <iosfwd>
#include <type_traits>
#include "CLHEP/Vector/ThreeVector.h"

namespace HepGeom {
  template<class T> class BasicVector3D {
  protected:
    T v_[3];

    BasicVector3D() { v_[0] = 0; v_[1] = 0; v_[2] = 0; }

  public:
    enum {
      X = 0,                 
      Y = 1,                 
      Z = 2,                 
      NUM_COORDINATES = 3,   
      SIZE = NUM_COORDINATES 
    };

    BasicVector3D(T x1, T y1, T z1) { v_[0] = x1; v_[1] = y1; v_[2] = z1; }

    BasicVector3D(const BasicVector3D<T> &) = default;

    template<typename U = T,
             typename = typename std::enable_if<!std::is_same<U,float>::value >::type>
    BasicVector3D(const BasicVector3D<float> & v) {
      v_[0] = v.x(); v_[1] = v.y(); v_[2] = v.z();
    }

    BasicVector3D(BasicVector3D<T> &&) = default;

    virtual ~BasicVector3D() = default;

    // -------------------------
    // Interface to "good old C"
    // -------------------------

    operator T * () { return v_; }

    operator const T * () const { return v_; }

    operator CLHEP::Hep3Vector () const { return CLHEP::Hep3Vector(x(),y(),z()); }

    // -----------------------------
    // General arithmetic operations
    // -----------------------------

    BasicVector3D<T> & operator= (const BasicVector3D<T> &) = default;
    BasicVector3D<T> & operator= (BasicVector3D<T> &&) = default;
    BasicVector3D<T> & operator+=(const BasicVector3D<T> & v) {
      v_[0] += v.v_[0]; v_[1] += v.v_[1]; v_[2] += v.v_[2]; return *this;
    }
    BasicVector3D<T> & operator-=(const BasicVector3D<T> & v) {
      v_[0] -= v.v_[0]; v_[1] -= v.v_[1]; v_[2] -= v.v_[2]; return *this;
    }
    BasicVector3D<T> & operator*=(double a) {
      v_[0] *= a; v_[1] *= a; v_[2] *= a; return *this;
    }
    BasicVector3D<T> & operator/=(double a) {
      v_[0] /= a; v_[1] /= a; v_[2] /= a; return *this;
    }

    // ------------
    // Subscripting
    // ------------

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

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

    // ------------------------------------
    // Cartesian coordinate system: x, y, z
    // ------------------------------------

    T x() const { return v_[0]; }
    T y() const { return v_[1]; }
    T z() const { return v_[2]; }

    void setX(T a) { v_[0] = a; }
    void setY(T a) { v_[1] = a; }
    void setZ(T a) { v_[2] = a; }

    void set(T x1, T y1, T z1) { v_[0] = x1; v_[1] = y1; v_[2] = z1; }

    // ------------------------------------------
    // Cylindrical coordinate system: rho, phi, z
    // ------------------------------------------

    T perp2() const { return x()*x()+y()*y(); }
    T perp() const { return std::sqrt(perp2()); }
    T rho() const { return perp(); }

    void setPerp(T rh) {
      T factor = perp();
      if (factor > 0) {
  factor = rh/factor; v_[0] *= factor; v_[1] *= factor;
      }
    }

    // ------------------------------------------
    // Spherical coordinate system: r, phi, theta
    // ------------------------------------------

    T mag2() const { return x()*x()+y()*y()+z()*z(); }
    T mag() const { return std::sqrt(mag2()); }
    T r() const { return mag(); }
    T phi() const {
      return x() == 0 && y() == 0 ? 0 : std::atan2(y(),x());
    }
    T theta() const {
      return x() == 0 && y() == 0 && z() == 0 ? 0 : std::atan2(perp(),z());
    }
    T cosTheta() const { T ma = mag(); return ma == 0 ? 1 : z()/ma; }

    T getR() const { return r(); }
    T getPhi() const { return phi(); }
    T getTheta() const { return theta(); }

    void setMag(T ma) {
      T factor = mag();
      if (factor > 0) {
  factor = ma/factor; v_[0] *= factor; v_[1] *= factor; v_[2] *= factor;
      }
    }
    void setR(T ma) { setMag(ma); }
    void setPhi(T ph) { T xy = perp(); setX(xy*std::cos(ph)); setY(xy*std::sin(ph)); }
    void setTheta(T th) {
      T ma = mag();
      T ph = phi();
      set(ma*std::sin(th)*std::cos(ph), ma*std::sin(th)*std::sin(ph), ma*std::cos(th));
    }

    // ---------------
    // Pseudo rapidity
    // ---------------

    T pseudoRapidity() const;
    T eta() const { return pseudoRapidity(); }
    T getEta() const { return pseudoRapidity(); }

    void setEta(T a);

    // -------------------
    // Combine two vectors
    // -------------------

    T dot(const BasicVector3D<T> & v) const {
      return x()*v.x()+y()*v.y()+z()*v.z();
    }

    BasicVector3D<T> cross(const BasicVector3D<T> & v) const {
      return BasicVector3D<T>(y()*v.z()-v.y()*z(),
            z()*v.x()-v.z()*x(),
            x()*v.y()-v.x()*y());
    }

    T perp2(const BasicVector3D<T> & v) const {
      T tot = v.mag2(), s = dot(v);
      return tot > 0 ? mag2()-s*s/tot : mag2();
    }

    T perp(const BasicVector3D<T> & v) const {
      return std::sqrt(perp2(v));
    }

    T angle(const BasicVector3D<T> & v) const;

    // ---------------
    // Related vectors
    // ---------------

    BasicVector3D<T> unit() const {
      T len = mag();
      return (len > 0) ?
  BasicVector3D<T>(x()/len, y()/len, z()/len) : BasicVector3D<T>();
    }

    BasicVector3D<T> orthogonal() const {
      T dx = x() < 0 ? -x() : x();
      T dy = y() < 0 ? -y() : y();
      T dz = z() < 0 ? -z() : z();
      if (dx < dy) {
  return dx < dz ?
    BasicVector3D<T>(0,z(),-y()) : BasicVector3D<T>(y(),-x(),0);
      }else{
  return dy < dz ?
    BasicVector3D<T>(-z(),0,x()) : BasicVector3D<T>(y(),-x(),0);
      }
    }

    // ---------
    // Rotations
    // ---------

    BasicVector3D<T> & rotateX(T a);
    BasicVector3D<T> & rotateY(T a);
    BasicVector3D<T> & rotateZ(T a);
    BasicVector3D<T> & rotate(T a, const BasicVector3D<T> & v);
  };

  /*************************************************************************
   *                                                                       *
   * Non-member functions for BasicVector3D<float>                         *
   *                                                                       *
   *************************************************************************/

  std::ostream &
  operator<<(std::ostream &, const BasicVector3D<float> &);

  std::istream &
  operator>>(std::istream &, BasicVector3D<float> &);

  inline BasicVector3D<float>
  operator+(const BasicVector3D<float> & v) { return v; }

  inline BasicVector3D<float>
  operator+(const BasicVector3D<float> & a, const BasicVector3D<float> & b) {
    return BasicVector3D<float>(a.x()+b.x(), a.y()+b.y(), a.z()+b.z());
  }

  inline BasicVector3D<float>
  operator-(const BasicVector3D<float> & v) {
    return BasicVector3D<float>(-v.x(), -v.y(), -v.z());
  }

  inline BasicVector3D<float>
  operator-(const BasicVector3D<float> & a, const BasicVector3D<float> & b) {
    return BasicVector3D<float>(a.x()-b.x(), a.y()-b.y(), a.z()-b.z());
  }

  inline BasicVector3D<float>
  operator*(const BasicVector3D<float> & v, double a) {
    return BasicVector3D<float>(v.x()*static_cast<float>(a), v.y()*static_cast<float>(a), v.z()*static_cast<float>(a));
  }

  inline float
  operator*(const BasicVector3D<float> & a, const BasicVector3D<float> & b) {
    return a.dot(b);
  }

  inline BasicVector3D<float>
  operator*(double a, const BasicVector3D<float> & v) {
    return BasicVector3D<float>(static_cast<float>(a)*v.x(), static_cast<float>(a)*v.y(), static_cast<float>(a)*v.z());
  }

  inline BasicVector3D<float>
  operator/(const BasicVector3D<float> & v, double a) {
    return BasicVector3D<float>(v.x()/static_cast<float>(a), v.y()/static_cast<float>(a), v.z()/static_cast<float>(a));
  }

  inline bool
  operator==(const BasicVector3D<float> & a, const BasicVector3D<float> & b) {
    return (a.x()==b.x() && a.y()==b.y() && a.z()==b.z());
  }

  inline bool
  operator!=(const BasicVector3D<float> & a, const BasicVector3D<float> & b) {
    return (a.x()!=b.x() || a.y()!=b.y() || a.z()!=b.z());
  }

  /*************************************************************************
   *                                                                       *
   * Non-member functions for BasicVector3D<double>                        *
   *                                                                       *
   *************************************************************************/

  std::ostream &
  operator<<(std::ostream &, const BasicVector3D<double> &);

  std::istream &
  operator>>(std::istream &, BasicVector3D<double> &);

  inline BasicVector3D<double>
  operator+(const BasicVector3D<double> & v) { return v; }

  inline BasicVector3D<double>
  operator+(const BasicVector3D<double> & a,const BasicVector3D<double> & b) {
    return BasicVector3D<double>(a.x()+b.x(), a.y()+b.y(), a.z()+b.z());
  }

  inline BasicVector3D<double>
  operator-(const BasicVector3D<double> & v) {
    return BasicVector3D<double>(-v.x(), -v.y(), -v.z());
  }

  inline BasicVector3D<double>
  operator-(const BasicVector3D<double> & a,const BasicVector3D<double> & b) {
    return BasicVector3D<double>(a.x()-b.x(), a.y()-b.y(), a.z()-b.z());
  }

  inline BasicVector3D<double>
  operator*(const BasicVector3D<double> & v, double a) {
    return BasicVector3D<double>(v.x()*a, v.y()*a, v.z()*a);
  }

  inline double
  operator*(const BasicVector3D<double> & a,const BasicVector3D<double> & b) {
    return a.dot(b);
  }

  inline BasicVector3D<double>
  operator*(double a, const BasicVector3D<double> & v) {
    return BasicVector3D<double>(a*v.x(), a*v.y(), a*v.z());
  }

  inline BasicVector3D<double>
  operator/(const BasicVector3D<double> & v, double a) {
    return BasicVector3D<double>(v.x()/a, v.y()/a, v.z()/a);
  }

  inline bool
  operator==(const BasicVector3D<double> & a, const BasicVector3D<double> & b)
  {
    return (a.x()==b.x() && a.y()==b.y() && a.z()==b.z());
  }

  inline bool
  operator!=(const BasicVector3D<double> & a, const BasicVector3D<double> & b)
  {
    return (a.x()!=b.x() || a.y()!=b.y() || a.z()!=b.z());
  }
} /* namespace HepGeom */

#endif /* BASIC_VECTOR3D_H */