ThreeVector.h

Go to the documentation of this file.

// -*- C++ -*-
//
// ThreeVector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2019 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_ThreeVector_H
#define ThePEG_ThreeVector_H

#include "ThreeVector.fh"
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Utilities/UnitIO.h"
#include <cassert>
#include <cmath>

namespace ThePEG {

template <typename Value>
class ThreeVector 
{  
private:
  using Value2 = decltype(sqr(std::declval<Value>()));

public:
  ThreeVector() 
    : theX(), theY(), theZ() {}

  ThreeVector(Value x, Value y, Value z)
    : theX(x), theY(y), theZ(z) {}

  template<typename ValueB>
  ThreeVector(const ThreeVector<ValueB> & v)
    : theX(v.x()), theY(v.y()), theZ(v.z()) {}

public:

  Value x() const { return theX; }
  Value y() const { return theY; }
  Value z() const { return theZ; }


  void setX(Value x)  {  theX = x; }
  void setY(Value y)  {  theY = y; }
  void setZ(Value z)  {  theZ = z; }

public:
  Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }

  Value  mag() const { return sqrt(mag2()); }

  Value2 perp2() const { return sqr(x()) + sqr(y()); }

  Value  perp()  const { return sqrt(perp2()); }

  template <typename U>
  auto dot(const ThreeVector<U> & a) const 
  -> decltype(this->x()*a.x())
  {
    return x()*a.x() + y()*a.y() + z()*a.z();
  }

  template <typename U>
  Value2 perp2(const ThreeVector<U> & p) const {
    const auto pMag2 = p.mag2();
    assert( pMag2 > ZERO );
    auto ss = this->dot(p);
    Value2 ret = mag2() - sqr(ss)/pMag2;
    if ( ret <= ZERO )
      ret = ZERO;
    return ret;
  }

  template <typename U>
  Value perp(const ThreeVector<U> & p) const {
    return sqrt(perp2(p));
  }


  double theta() const {
    assert(!(x() == ZERO && y() == ZERO && z() == ZERO));
    return atan2(perp(),z());
  }

  double phi()   const {
    return atan2(y(),x());
  }

  void setTheta(double th) {
    double ma  = mag();
    double ph  = phi();
    setX(ma*sin(th)*cos(ph));
    setY(ma*sin(th)*sin(ph));
    setZ(ma*cos(th));
  }

  void setPhi(double ph) {
    double xy = perp();
    setX(xy*cos(ph));
    setY(xy*sin(ph));
  }

  ThreeVector<double> unit() const {
    Value2 mg2 = mag2();
    assert(mg2 > ZERO);
    Value mg = sqrt(mg2);
    return {x()/mg, y()/mg, z()/mg};
  }
  
  ThreeVector<Value> orthogonal() const {
    Value xx = abs(x());
    Value yy = abs(y());
    Value zz = abs(z());
    using TVec = ThreeVector<Value>;
    if (xx < yy) {
      return xx < zz ? TVec{ZERO,z(),-y()} : TVec{y(),-x(),ZERO};
    } else {
      return yy < zz ? TVec{-z(),ZERO,x()} : TVec{y(),-x(),ZERO};
    }
  }

  template <typename U>
  double deltaPhi  (const ThreeVector<U> & v2) const {
    double dphi = v2.phi() - phi();
    if ( dphi > Constants::pi ) {
      dphi -= Constants::twopi;
    } else if ( dphi <= -Constants::pi ) {
      dphi += Constants::twopi;
    }
    return dphi;
  } 

  template <typename U>
  ThreeVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
    if (angle == 0.0) 
      return *this;
    const U ll = axis.mag();
    assert( ll > ZERO );

    const double sa = sin(angle), ca = cos(angle);
    const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
    const Value  xx  = x(), yy = y(), zz = z(); 

    setX((ca+(1-ca)*dx*dx)     * xx
         +((1-ca)*dx*dy-sa*dz) * yy
         +((1-ca)*dx*dz+sa*dy) * zz
         );
    setY(((1-ca)*dy*dx+sa*dz)  * xx
         +(ca+(1-ca)*dy*dy)    * yy
         +((1-ca)*dy*dz-sa*dx) * zz
         );
    setZ(((1-ca)*dz*dx-sa*dy)  * xx
         +((1-ca)*dz*dy+sa*dx) * yy
         +(ca+(1-ca)*dz*dz)    * zz
         );
    return *this;
  }


  ThreeVector<Value> & rotateUz (const Axis & axis) {
    Axis ax = axis.unit();
    double u1 = ax.x();
    double u2 = ax.y();
    double u3 = ax.z();
    double up = u1*u1 + u2*u2;
    if (up>0) {
      up = sqrt(up);
      Value px = x(),  py = y(),  pz = z();
      setX( (u1*u3*px - u2*py)/up + u1*pz );
      setY( (u2*u3*px + u1*py)/up + u2*pz );
      setZ(    -up*px +             u3*pz );
    }
    else if (u3 < 0.) {
      setX(-x());
      setZ(-z()); 
    }
    return *this;
  }

  ThreeVector<Value> & rotateUzBack (const Axis & axis) {
    Axis ax = axis.unit();
    double u1 = ax.x();
    double u2 = ax.y();
    double u3 = ax.z();
    double up = u1*u1 + u2*u2;
    if (up>0) {
      up = sqrt(up);
      Value px = x(),  py = y(),  pz = z();
      setX( ( u1*u3*px + u2*u3*py)/up - up*pz );
      setY( (-u2*px    + u1*py)/up );
      setZ(   u1*px    + u2*py        + u3*pz );
    }
    else if (u3 < 0.) {
      setX(-x());
      setZ(-z()); 
    }
    return *this;
  }

  template <typename U>
  auto cross(const ThreeVector<U> & a) const 
  -> ThreeVector<decltype(this->y()*a.z())>
  {
    return { y()*a.z()-z()*a.y(),
            -x()*a.z()+z()*a.x(),
             x()*a.y()-y()*a.x() };
  }
  
  public:  

  bool operator==(const ThreeVector<Value> & a) const {
    return (theX == a.x() && theY == a.y() && theZ == a.z());
  }
  bool operator!=(const ThreeVector<Value> & a) const {
    return !(*this == a);
  }
  bool almostEqual(const ThreeVector<Value> & a, double threshold = 1e-04) const {
    return ((std::abs(theX - a.x()) < threshold) && (std::abs(theY - a.y()) < threshold) && (std::abs(theZ - a.z()) < threshold));
  }
  bool almostUnequal(const ThreeVector<Value> & a, double threshold = 1e-04) const {
    return ! this->almostEqual(a, threshold);
  }
  
public:  

  ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
    theX += a.x();
    theY += a.y();
    theZ += a.z();
    return *this;
  }

  ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
    theX -= a.x();
    theY -= a.y();
    theZ -= a.z();
    return *this;
  }

  ThreeVector<Value> & operator*=(double a) {
    theX *= a;
    theY *= a;
    theZ *= a;
    return *this;
  }

  ThreeVector<Value> & operator/=(double a) {
    theX /= a;
    theY /= a;
    theZ /= a;
    return *this;
  }
  
  template <typename U>
  double cosTheta(const ThreeVector<U> & q) const {
    auto ptot = mag()*q.mag();
    assert( ptot > ZERO );
    double arg = dot(q)/ptot;
    if     (arg >  1.0) arg =  1.0;
    else if(arg < -1.0) arg = -1.0;
    return arg;
  }
  
  template <typename U>
  double angle(const ThreeVector<U> & v) const {
    return acos(cosTheta(v));
  }

private:

  Value theX;
  Value theY;
  Value theZ;
};

inline ostream & 
operator<< (ostream & os, const ThreeVector<double> & v)
{
  return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
}


template <typename Value>
inline ThreeVector<Value>
operator+(ThreeVector<Value> a, 
          const ThreeVector<Value> & b)
{
  return a += b;
}

template <typename Value>
inline ThreeVector<Value>
operator-(ThreeVector<Value> a, 
          const ThreeVector<Value> & b)
{
  return a -= b;
}

template <typename Value>
inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
  return {-v.x(),-v.y(),-v.z()};
}

template <typename Value>
inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
  return v *= a;
}

template <typename Value>
inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
  return v *= a;
}

template <typename ValueA, typename ValueB>
inline auto operator*(ValueB a, ThreeVector<ValueA> v) 
-> ThreeVector<decltype(a*v.x())>
{
  return {a*v.x(), a*v.y(), a*v.z()};
}

template <typename ValueA, typename ValueB>
inline auto operator*(ThreeVector<ValueA> v, ValueB a) 
-> ThreeVector<decltype(v.x()*a)>
{
  return {v.x()*a, v.y()*a, v.z()*a};
}

template <typename ValueA, typename ValueB>
inline auto operator*(const ThreeVector<ValueA> & a, 
                      const ThreeVector<ValueB> & b)
-> decltype(a.x()*b.x())
{
  return a.dot(b);
}

template <typename Value>
ThreeVector<double> unitVector(const ThreeVector<Value> & v) {
  return v.unit();
}


template <typename OStream, typename UT, typename Value>
void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
  os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
}

template <typename IStream, typename UT, typename Value>
void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
  Value x, y, z;
  is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
  p = ThreeVector<Value>(x, y, z);
}

}

#endif /* ThePEG_ThreeVector_H */