LorentzVector.h

Go to the documentation of this file.

// -*- C++ -*-
//
// LorentzVector.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_LorentzVector_H
#define ThePEG_LorentzVector_H

#include "LorentzVector.fh"
#include "ThePEG/Utilities/Direction.h"
#include "ThePEG/Utilities/UnitIO.h"
#include "LorentzRotation.h"
#include "ThreeVector.h"

#ifdef NDEBUG
#define ERROR_IF(condition,message) if (false) {}
#else
#define ERROR_IF(condition,message) \
  if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror)
#endif

namespace ThePEG {

template <typename Value> class LorentzVector; 

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

public:
  LorentzVector() 
    : theX(), theY(), theZ(), theT() {}

  LorentzVector(Value x, Value y, Value z, Value t)
    : theX(x), theY(y), theZ(z), theT(t) {}

  LorentzVector(const ThreeVector<Value> & v, Value t)
    : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}

  template<typename U>
  LorentzVector(const LorentzVector<U> & v)
    : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}

  template <typename ValueB>
  LorentzVector<Value> & operator=(const LorentzVector<ValueB> & b) {
    setX(b.x());
    setY(b.y());
    setZ(b.z());
    setT(b.t());
    return *this;
  }

public:

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


  void setX(Value x)  {  theX = x; }
  void setY(Value y)  {  theY = y; }
  void setZ(Value z)  {  theZ = z; }
  void setT(Value t)  {  theT = t; }
  void setE(Value e)  {  setT(e);  }

public:
  ThreeVector<Value> vect() const {
    return ThreeVector<Value>(x(),y(),z());
  }

  operator ThreeVector<Value>() const { return vect(); }
  
  void setVect(const ThreeVector<Value> & p) {
    theX = p.x();
    theY = p.y();
    theZ = p.z();
  } 

public:
  LorentzVector<Value> conjugate() const 
  {
    return LorentzVector<Value>(conj(x()),conj(y()),conj(z()),conj(t()));
  }

  Value2 m2() const 
  { 
    return (t()-z())*(t()+z()) - sqr(x()) - sqr(y()); 
  }

  Value2 m2(const LorentzVector<Value> & a) const {
    Value tt(a.t()+t()),zz(a.z()+z());
    return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
  }

  Value  m() const 
  {
    Value2 tmp = m2();
    return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
  }

  Value2 mt2()  const { return (t()-z())*(t()+z()); }

  Value  mt()  const 
  { 
    Value2 tmp = mt2();
    return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
  }

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

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

  template <typename U>
  Value2 perp2(const ThreeVector<U> & p) const 
  {
    return vect().perp2(p);
  }

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

  Value2 et2() const 
  {
    Value2 pt2 = vect().perp2();
    return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z());
  }

  Value et() const 
  {
    Value2 etet = et2();
    return e() < Value() ? -sqrt(etet) : sqrt(etet);
  }

  Value2 et2(const ThreeVector<double> & v) const 
  {
    Value2 pt2 = vect().perp2(v);
    Value pv = vect().dot(v.unit());
    return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv);
  }

  Value et(const ThreeVector<double> & v) const 
  {
    Value2 etet = et2(v);
    return e() < Value() ? -sqrt(etet) : sqrt(etet);
  }


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

  Value  rho()   const { return sqrt(rho2()); }

  void setRho(Value newRho) 
  { 
    Value oldRho = rho();
    if (oldRho == Value()) 
      return;
    double factor = newRho / oldRho;
    setX(x()*factor);
    setY(y()*factor);
    setZ(z()*factor);
  }

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

  double cosTheta() const 
  {
    Value ptot = rho();
    assert( ptot > Value() );
    return z() / ptot;
  }

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

  double eta() const {
    Value m = rho();
    if ( m ==  Value() ) return  0.0;
    Value pt = max(Constants::epsilon*m, perp());
    double rap = log((m + abs(z()))/pt);
    return z() > ZERO? rap: -rap;
  }

  double angle(const LorentzVector<Value> & w) const 
  {
    return vect().angle(w.vect());
  }

  double rapidity() const {
    if ( z() == ZERO ) return 0.0;
    ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
    Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2()));
    double rap = log((t() + abs(z()))/pt);
    return z() > ZERO? rap: -rap;
  }

  double rapidity(const Axis & ref) const {
    double r = ref.mag2();
    ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity");
    Value vdotu = vect().dot(ref)/sqrt(r);
    if ( vdotu == ZERO ) return 0.0;
    ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
    Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2()));
    double rap = log((t() + abs(z()))/pt);
    return z() > ZERO? rap: -rap;
  }

  Boost boostVector() const {
    if (t() == Value()) {
      if (rho2() == Value2()) 
        return Boost();
      else 
        ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result");
    }
    // result will make analytic sense but is physically meaningless
    ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector");
    return vect() * (1./t());
  }
  
  Boost findBoostToCM() const 
  {
    return -boostVector();
  }

  Value plus()  const { return t() + z(); }
  Value minus() const { return t() - z(); }

  bool isNear(const LorentzVector<Value> & w, double epsilon) const 
  {
    Value2 limit = abs(vect().dot(w.vect()));
    limit += 0.25 * sqr( t() + w.t() );
    limit *= sqr(epsilon);
    Value2 delta = (vect() - w.vect()).mag2();
    delta +=  sqr( t() - w.t() );
    return (delta <= limit);
  }
  
  LorentzVector<Value> & transform(const SpinOneLorentzRotation & m) 
  {
    return *this = m.operator*(*this);
  }

  LorentzVector<Value> & operator*=(const SpinOneLorentzRotation & m) 
  {
    return transform(m);
  }

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


public:

  LorentzVector<Value> & 
  boost(double bx, double by, double bz, double gamma=-1.) 
  {
    const double b2 = bx*bx + by*by + bz*bz;
    if ( b2 == 0.0 ) return *this;
    if ( gamma < 0.0 ) {
        gamma = 1.0 / sqrt(1.0 - b2);
    }
    const Value bp = bx*x() + by*y() + bz*z();
    const double gamma2 = (gamma - 1.0)/b2;
    
    setX(x() + gamma2*bp*bx + gamma*bx*t());
    setY(y() + gamma2*bp*by + gamma*by*t());
    setZ(z() + gamma2*bp*bz + gamma*bz*t());
    setT(gamma*(t() + bp));
    return *this;
  }
  
  LorentzVector<Value> &  boost(Boost b, double gamma=-1.) {
    return boost(b.x(), b.y(), b.z(),gamma);
  }

  LorentzVector<Value> & rotateX (double phi) {
    double sinphi = sin(phi);
    double cosphi = cos(phi);
    Value  ty = y() * cosphi - z() * sinphi;
    theZ      = z() * cosphi + y() * sinphi;
    theY = ty;
    return *this;
  }
  
  LorentzVector<Value> & rotateY (double phi) {
    double sinphi = sin(phi);
    double cosphi = cos(phi);
    Value  tz = z() * cosphi - x() * sinphi;
    theX      = x() * cosphi + z() * sinphi;
    theZ = tz;
    return *this;
  }
  
  LorentzVector<Value> & rotateZ (double phi) {
    double sinphi = sin(phi);
    double cosphi = cos(phi);
    Value  tx = x() * cosphi - y() * sinphi;
    theY      = y() * cosphi + x() * sinphi;
    theX = tx;
    return *this;
  }
  
  LorentzVector<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;
  }
  
  template <typename U>
  LorentzVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
    if (angle == 0.0) 
      return *this;
    const U ll = axis.mag();
    assert( ll > U() );

    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;
  }




public:

  LorentzVector<Complex> & operator+=(const LorentzVector<complex<QtyDouble> > & a) {
    theX += a.x();
    theY += a.y();
    theZ += a.z();
    theT += a.t();
    return *this;
  }
  
  template <typename ValueB>
  LorentzVector<Value> & operator+=(const LorentzVector<ValueB> & a) {
    theX += a.x();
    theY += a.y();
    theZ += a.z();
    theT += a.t();
    return *this;
  }
  
  LorentzVector<Complex> & operator-=(const LorentzVector<complex<QtyDouble> > & a) {
    theX -= Complex(a.x());
    theY -= Complex(a.y());
    theZ -= Complex(a.z());
    theT -= Complex(a.t());
    return *this;
  }
  
  template <typename ValueB>
  LorentzVector<Value> & operator-=(const LorentzVector<ValueB> & a) {
    theX -= a.x();
    theY -= a.y();
    theZ -= a.z();
    theT -= a.t();
    return *this;
  }

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

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

private:

  Value theX;
  Value theY;
  Value theZ;
  Value theT;
};


template <typename Value>
inline LorentzVector<double>
operator/(const LorentzVector<Value> & v, Value a) {
  return LorentzVector<double>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
}

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

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

template <typename ValueA, typename ValueB>
inline LorentzVector<ValueA>
operator+(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
  return a += b;
}

template <typename ValueA, typename ValueB>
inline LorentzVector<ValueA>
operator-(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
  return a -= b;
}

template <typename Value>
inline LorentzVector<Value>
operator*(const LorentzVector<Value> & a, double b) {
  return LorentzVector<Value>(a.x()*b, a.y()*b, a.z()*b, a.t()*b);
}

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

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

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

template <typename ValueA, typename ValueB>
inline auto operator/(const LorentzVector<ValueA> & v, ValueB b) 
-> LorentzVector<decltype(v.x()/b)>
{
  return {v.x()/b, v.y()/b, v.z()/b, v.t()/b};
}


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


template <typename Value>
inline bool
operator==(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
  return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
}

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

template <typename Value>
inline Value dirPlus(const LorentzVector<Value> & p) {
  return Direction<0>::pos()? p.plus(): p.minus();
}

template <typename Value>
inline Value dirMinus(const LorentzVector<Value> & p) {
  return Direction<0>::neg()? p.plus(): p.minus();
}

template <typename Value>
inline Value dirZ(const LorentzVector<Value> & p) {
  return Direction<0>::dir()*p.z();
}

template <typename Value>
inline double dirTheta(const LorentzVector<Value> & p) {
  return Direction<0>::pos()? p.theta(): Constants::pi - p.theta();
}

template <typename Value>
inline double dirCosTheta(const LorentzVector<Value> & p) {
  return Direction<0>::pos()? p.cosTheta():  -p.cosTheta();
}

template <typename Value>
inline ThreeVector<Value> dirBoostVector(const LorentzVector<Value> & p) {
  ThreeVector<Value> b(p.boostVector());
  if ( Direction<0>::neg() ) b.setZ(-b.z());
  return b;
}

template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus, Value x, Value y) {
  LorentzVector<Value> r(x, y, 0.5*(plus-minus), 0.5*(plus+minus));
  return r;
}

template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus) {
  // g++-3.3 has a problem with using Value() directly
  // gcc-bug c++/3650, fixed in 3.4
  static const Value zero = Value();
  LorentzVector<Value> r(zero, zero,
                         0.5*(plus-minus), 0.5*(plus+minus));
  return r;
}

}


// delayed header inclusion to break inclusion loop:
// LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec
#include "Transverse.h"



namespace ThePEG {

template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus, Transverse<Value> pt) {
  LorentzVector<Value> r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus));
  return r;
}

template <typename Value>
inline LorentzVector<Value>
lightConeDir(Value plus, Value minus,
             Value x = Value(), Value y = Value()) {
  LorentzVector<Value> r(x, y, Direction<0>::dir()*0.5*(plus - minus),
                  0.5*(plus + minus));
  return r;
}

template <typename Value>
inline LorentzVector<Value>
lightConeDir(Value plus, Value minus, Transverse<Value> pt) {
  LorentzVector<Value> r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus),
                  0.5*(plus + minus));
  return r;

}

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

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


}

#undef ERROR_IF
#endif /* ThePEG_LorentzVector_H */