Lorentz5Vector.h

// -*- C++ -*-
//
// Lorentz5Vector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 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_Lorentz5Vector_H
#define ThePEG_Lorentz5Vector_H
 
// This is the declaration of the Lorentz5vector class.
 
#include "LorentzVector.h"
#include "Lorentz5Vector.fh"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/Utilities/Direction.h"
#include "ThePEG/Utilities/UnitIO.h"
#include "LorentzRotation.h"
 
namespace ThePEG {
 
template <typename Value>
class Lorentz5Vector: public LorentzVector<Value> {
 
public:
 
  using Value2 = decltype(sqr(std::declval<Value>()));
 
public:
 
  Value x() const { return LorentzVector<Value>::x(); }
  Value y() const { return LorentzVector<Value>::y(); }
  Value z() const { return LorentzVector<Value>::z(); }
  Value t() const { return LorentzVector<Value>::t(); }
 
public:
 
  Lorentz5Vector() : mm() {}
 
  Lorentz5Vector(Value m) 
    : LorentzVector<Value>(Value(), Value(), Value(), m), mm(m) {}
 
  Lorentz5Vector(Value x, Value y, Value z, Value t = Value())
    : LorentzVector<Value>(x, y, z, t) { rescaleMass(); }
 
  Lorentz5Vector(Value x, Value y, Value z, Value t, Value tau)
    : LorentzVector<Value>(x, y, z, t), mm(tau) {}
 
  Lorentz5Vector(const ThreeVector<Value> & p, Value e)
    : LorentzVector<Value>(p, e) { rescaleMass(); }
 
  Lorentz5Vector(Value m, const ThreeVector<Value> & p) 
    : LorentzVector<Value>(p, sqrt(p.mag2() + m*m)), mm(m) {}
 
  Lorentz5Vector(const ThreeVector<Value> & p, Value t, Value tau)
    : LorentzVector<Value>(p, t), mm(tau) {}
 
  Lorentz5Vector(const LorentzVector<Value> & p, Value m) 
    : LorentzVector<Value>(p), mm(m) {}
 
  Lorentz5Vector(const LorentzVector<Value> & p)
    : LorentzVector<Value>(p) { rescaleMass(); } 
 
  template<class U>
  Lorentz5Vector(const Lorentz5Vector<U> & p)
    : LorentzVector<Value>(p), mm(p.m) {}
 
  void setTau(Value a) { mm = a; }
 
  void setMass(Value a) { mm = a; }
 
  Lorentz5Vector & operator=(const LorentzVector<Value> & q) {
    LorentzVector<Value>::operator=(q);
    return *this;
  }
 
  void rescaleEnergy() {
    LorentzVector<Value>::setT(sqrt(LorentzVector<Value>::vect().mag2() + mass2()));
  }
 
  void rescaleRho() {
    LorentzVector<Value>::setRho(sqrt(t()*t() - mass2()));
  }
 
  void rescaleMass() {
    mm = LorentzVector<Value>::m();
  }
 
  double massError() const {
    return sqrt(abs(Math::relativeError(mass2(), 
                                        LorentzVector<Value>::m2())));
  }
 
  double energyError() const {
    return sqrt(abs(Math::relativeError(t()*t(), mass2() 
                                        + LorentzVector<Value>::vect().mag2())));
  }
 
  double rhoError() const {
    return sqrt(abs(Math::relativeError(LorentzVector<Value>::vect().mag2(), 
                                        t()*t() - mass2())));
  }
 
  Value2 mass2() const { return mm > Value() ? mm*mm: -mm*mm; }
 
  Value2 tau2() const { return mass2(); }
 
  Value mass() const { return mm; }
 
 
  Value tau() const { return mass(); }
 
  Value dirPlus() const {
    return Direction<0>::pos() ? 
      LorentzVector<Value>::plus() 
      : 
      LorentzVector<Value>::minus();
  }
 
  Value dirMinus() const {
    return Direction<0>::neg() ? 
      LorentzVector<Value>::plus() 
      : 
      LorentzVector<Value>::minus();
  }
 
  Lorentz5Vector & transform(const LorentzRotation & r) 
  {
    LorentzVector<Value>::transform(r.one());
    return *this;
  }
 
private:
 
  Value mm;
 
};
 
template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const Lorentz5Vector<T> & p, UT & u) {
  os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
     << ounit(p.e(), u) << ounit(p.mass(), u);
}
 
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, Lorentz5Vector<T> & p, UT & u) {
  T x, y, z, e, mass;
  is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u)
     >> iunit(mass, u);
  p = Lorentz5Vector<T>(x, y, z, e, mass);
}
 
 
 
template <typename ValueA, typename ValueB>
inline auto
operator*(const Lorentz5Vector<ValueA> & a, const Lorentz5Vector<ValueB> & b)
-> decltype(a.dot(b))
{
  return a.dot(b);
}
 
template <typename ValueA, typename ValueB>
inline auto
operator*(const LorentzVector<ValueA> & a, const Lorentz5Vector<ValueB> & b)
-> decltype(a.dot(b))
{
  return a.dot(b);
}
 
template <typename ValueA, typename ValueB>
inline auto
operator*(const Lorentz5Vector<ValueA> & a, const LorentzVector<ValueB> & b)
-> decltype(a.dot(b))
{
  return a.dot(b);
}
 
template <typename Value>
inline auto
operator*(const Lorentz5Vector<Value> & a, const Lorentz5Vector<Value> & b)
-> decltype(a.dot(b))
{
  return a.dot(b);
}
}
 
#endif /* ThePEG_Particle_H */