#ifndef ROOT_Math_GenVector_PxPyPzM4D
#define ROOT_Math_GenVector_PxPyPzM4D 1
#ifndef ROOT_Math_GenVector_eta
#include "Math/GenVector/eta.h"
#endif
#ifndef ROOT_Math_GenVector_GenVector_exception
#include "Math/GenVector/GenVector_exception.h"
#endif
#include <cmath>
namespace ROOT {
namespace Math {
template <class ScalarType = double>
class PxPyPzM4D {
public :
typedef ScalarType Scalar;
PxPyPzM4D() : fX(0), fY(0), fZ(0), fM(0) { }
PxPyPzM4D(Scalar px, Scalar py, Scalar pz, Scalar m) :
fX(px), fY(py), fZ(pz), fM(m) {
if (fM < 0) RestrictNegMass();
}
template <class CoordSystem>
explicit PxPyPzM4D(const CoordSystem & v) :
fX( v.X() ), fY( v.Y() ), fZ( v.Z() ), fM( v.M() )
{ }
PxPyPzM4D(const PxPyPzM4D & v) :
fX(v.fX), fY(v.fY), fZ(v.fZ), fM(v.fM) { }
PxPyPzM4D & operator = (const PxPyPzM4D & v) {
fX = v.fX;
fY = v.fY;
fZ = v.fZ;
fM = v.fM;
return *this;
}
template <class AnyCoordSystem>
PxPyPzM4D & operator = (const AnyCoordSystem & v) {
fX = v.X();
fY = v.Y();
fZ = v.Z();
fM = v.M();
return *this;
}
void SetCoordinates( const Scalar src[] ) {
fX=src[0]; fY=src[1]; fZ=src[2]; fM=src[3];
if (fM < 0) RestrictNegMass();
}
void GetCoordinates( Scalar dest[] ) const
{ dest[0] = fX; dest[1] = fY; dest[2] = fZ; dest[3] = fM; }
void SetCoordinates(Scalar px, Scalar py, Scalar pz, Scalar m) {
fX=px; fY=py; fZ=pz; fM=m;
if (fM < 0) RestrictNegMass();
}
void GetCoordinates(Scalar& px, Scalar& py, Scalar& pz, Scalar& m) const
{ px=fX; py=fY; pz=fZ; m=fM;}
Scalar Px() const { return fX;}
Scalar Py() const { return fY;}
Scalar Pz() const { return fZ;}
Scalar M() const { return fM; }
Scalar X() const { return fX;}
Scalar Y() const { return fY;}
Scalar Z() const { return fZ;}
Scalar E() const { return std::sqrt(E2() ); }
Scalar T() const { return E();}
Scalar P2() const { return fX*fX + fY*fY + fZ*fZ; }
Scalar P() const { return std::sqrt(P2()); }
Scalar R() const { return P(); }
Scalar M2() const {
return ( fM >= 0 ) ? fM*fM : -fM*fM;
}
Scalar Mag2() const { return M2(); }
Scalar Mag() const { return M(); }
Scalar E2() const {
Scalar e2 = P2() + M2();
return e2 > 0 ? e2 : 0;
}
Scalar Pt2() const { return fX*fX + fY*fY;}
Scalar Perp2() const { return Pt2();}
Scalar Pt() const { return std::sqrt(Perp2());}
Scalar Perp() const { return Pt();}
Scalar Rho() const { return Pt();}
Scalar Mt2() const { return E2() - fZ*fZ; }
Scalar Mt() const {
Scalar mm = Mt2();
if (mm >= 0) {
return std::sqrt(mm);
} else {
GenVector::Throw ("PxPyPzM4D::Mt() - Tachyonic:\n"
" Pz^2 > E^2 so the transverse mass would be imaginary");
return -std::sqrt(-mm);
}
}
Scalar Et2() const {
Scalar pt2 = Pt2();
return pt2 == 0 ? 0 : E2() * pt2/( pt2 + fZ*fZ );
}
Scalar Et() const {
Scalar etet = Et2();
return std::sqrt(etet);
}
Scalar Phi() const {
return (fX == 0.0 && fY == 0.0) ? 0.0 : std::atan2(fY,fX);
}
Scalar Theta() const {
return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : std::atan2(Pt(),fZ);
}
Scalar Eta() const {
return Impl::Eta_FromRhoZ ( Pt(), fZ);
}
void SetPx( Scalar px) {
fX = px;
}
void SetPy( Scalar py) {
fY = py;
}
void SetPz( Scalar pz) {
fZ = pz;
}
void SetM( Scalar m) {
fM = m;
if (fM < 0) RestrictNegMass();
}
void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
void Negate( ) {
fX = -fX;
fY = -fY;
fZ = -fZ;
GenVector::Throw ("PxPyPzM4D::Negate - cannot negate the energy - can negate only the spatial components");
}
void Scale( const Scalar & a) {
fX *= a;
fY *= a;
fZ *= a;
fM *= a;
}
bool operator == (const PxPyPzM4D & rhs) const {
return fX == rhs.fX && fY == rhs.fY && fZ == rhs.fZ && fM == rhs.fM;
}
bool operator != (const PxPyPzM4D & rhs) const {return !(operator==(rhs));}
Scalar x() const { return X(); }
Scalar y() const { return Y(); }
Scalar z() const { return Z(); }
Scalar t() const { return E(); }
#if defined(__MAKECINT__) || defined(G__DICTIONARY)
void SetPt(Scalar pt);
void SetEta(Scalar eta);
void SetPhi(Scalar phi);
void SetE(Scalar t);
#endif
private:
inline void RestrictNegMass() {
if ( fM >=0 ) return;
if ( P2() - fM*fM < 0 ) {
GenVector::Throw("PxPyPzM4D::unphysical value of mass, set to closest physical value");
fM = - P();
}
return;
}
ScalarType fX;
ScalarType fY;
ScalarType fZ;
ScalarType fM;
};
}
}
#ifndef ROOT_Math_GenVector_PxPyPzE4D
#include "Math/GenVector/PxPyPzE4D.h"
#endif
#ifndef ROOT_Math_GenVector_PtEtaPhiM4D
#include "Math/GenVector/PtEtaPhiM4D.h"
#endif
namespace ROOT {
namespace Math {
template <class ScalarType>
inline void PxPyPzM4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
*this = PxPyPzE4D<Scalar> (px, py, pz, e);
}
#if defined(__MAKECINT__) || defined(G__DICTIONARY)
template <class ScalarType>
inline void PxPyPzM4D<ScalarType>::SetPt(ScalarType pt) {
GenVector_exception e("PxPyPzM4D::SetPt() is not supposed to be called");
throw e;
PtEtaPhiE4D<ScalarType> v(*this); v.SetPt(pt); *this = PxPyPzM4D<ScalarType>(v);
}
template <class ScalarType>
inline void PxPyPzM4D<ScalarType>::SetEta(ScalarType eta) {
GenVector_exception e("PxPyPzM4D::SetEta() is not supposed to be called");
throw e;
PtEtaPhiE4D<ScalarType> v(*this); v.SetEta(eta); *this = PxPyPzM4D<ScalarType>(v);
}
template <class ScalarType>
inline void PxPyPzM4D<ScalarType>::SetPhi(ScalarType phi) {
GenVector_exception e("PxPyPzM4D::SetPhi() is not supposed to be called");
throw e;
PtEtaPhiE4D<ScalarType> v(*this); v.SetPhi(phi); *this = PxPyPzM4D<ScalarType>(v);
}
template <class ScalarType>
inline void PxPyPzM4D<ScalarType>::SetE(ScalarType energy) {
GenVector_exception e("PxPyPzM4D::SetE() is not supposed to be called");
throw e;
PxPyPzE4D<ScalarType> v(*this); v.SetE(energy);
*this = PxPyPzM4D<ScalarType>(v);
}
#endif // endif __MAKE__CINT || G__DICTIONARY
}
}
#endif // ROOT_Math_GenVector_PxPyPzM4D