PxPyPzM4D.h

Go to the documentation of this file.

// @(#)root/mathcore:$Id: 464c29f33a8bbd8462a3e15b7e4c30c6f5b74a30 $
// Authors: W. Brown, M. Fischler, L. Moneta    2005
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
 *                                                                    *
 *                                                                    *
 **********************************************************************/
 
// Header file for class PxPyPzM4D
//
// Created by: fischler at Wed Jul 20   2005
//   (starting from PxPyPzM4D by moneta)
//
// Last update: $Id: 464c29f33a8bbd8462a3e15b7e4c30c6f5b74a30 $
//
#ifndef ROOT_MathX_GenVectorX_PxPyPzM4D
#define ROOT_MathX_GenVectorX_PxPyPzM4D 1
 
#include "MathX/GenVectorX/eta.h"
 
#include "MathX/GenVectorX/GenVector_exception.h"
 
#include "MathX/GenVectorX/MathHeaders.h"
 
#include "MathX/GenVectorX/AccHeaders.h"
 
using namespace ROOT::ROOT_MATH_ARCH;
 
#include <cmath>
 
namespace ROOT {
 
namespace ROOT_MATH_ARCH {
 
//__________________________________________________________________________________________
/**
    Class describing a 4D coordinate system
    or momentum-energy vectors stored as (Px, Py, Pz, M).
    This system is useful to describe ultra-relativistic particles
    (like electrons at LHC) to avoid numerical errors evaluating the mass
    when E >>> m
    The metric used is (-,-,-,+)
    Spacelike particles (M2 < 0) are described with negative mass values,
    but in this case m2 must always be less than P2 to preserve a positive value of E2
 
    @ingroup GenVectorX
 
    @see GenVectorX
*/
 
template <class ScalarType = double>
class PxPyPzM4D {
 
public:
   typedef ScalarType Scalar;
 
   // --------- Constructors ---------------
 
   /**
      Default constructor  with x=y=z=m=0
   */
   constexpr PxPyPzM4D() noexcept = default;
 
   /**
      Constructor  from x, y , z , m values
   */
   PxPyPzM4D(Scalar px, Scalar py, Scalar pz, Scalar m) : fX(px), fY(py), fZ(pz), fM(m)
   {
 
      if (fM < 0)
         RestrictNegMass();
   }
 
   /**
      construct from any 4D  coordinate system class
      implementing X(), Y(), X() and M()
   */
   template <class CoordSystem>
   explicit PxPyPzM4D(const CoordSystem &v) : fX(v.X()), fY(v.Y()), fZ(v.Z()), fM(v.M())
   {
   }
 
   /**
      construct from any 4D  coordinate system class
      implementing X(), Y(), X() and M()
   */
   template <class AnyCoordSystem>
   PxPyPzM4D &operator=(const AnyCoordSystem &v)
   {
      fX = v.X();
      fY = v.Y();
      fZ = v.Z();
      fM = v.M();
      return *this;
   }
 
   /**
      Set internal data based on an array of 4 Scalar numbers
   */
   void SetCoordinates(const Scalar src[])
   {
      fX = src[0];
      fY = src[1];
      fZ = src[2];
      fM = src[3];
      if (fM < 0)
         RestrictNegMass();
   }
 
   /**
      get internal data into an array of 4 Scalar numbers
   */
   void GetCoordinates(Scalar dest[]) const
   {
      dest[0] = fX;
      dest[1] = fY;
      dest[2] = fZ;
      dest[3] = fM;
   }
 
   /**
      Set internal data based on 4 Scalar numbers
   */
   void SetCoordinates(Scalar px, Scalar py, Scalar pz, Scalar m)
   {
      fX = px;
      fY = py;
      fZ = pz;
      fM = m;
      if (fM < 0)
         RestrictNegMass();
   }
 
   /**
      get internal data into 4 Scalar numbers
   */
   void GetCoordinates(Scalar &px, Scalar &py, Scalar &pz, Scalar &m) const
   {
      px = fX;
      py = fY;
      pz = fZ;
      m = fM;
   }
 
   // --------- Coordinates and Coordinate-like Scalar properties -------------
 
   // cartesian (Minkowski)coordinate accessors
 
   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; }
 
   // other coordinate representation
   /**
      Energy
    */
   Scalar E() const { return math_sqrt(E2()); }
 
   Scalar T() const { return E(); }
 
   /**
      squared magnitude of spatial components
   */
   Scalar P2() const { return fX * fX + fY * fY + fZ * fZ; }
 
   /**
      magnitude of spatial components (magnitude of 3-momentum)
   */
   Scalar P() const { return math_sqrt(P2()); }
   Scalar R() const { return P(); }
 
   /**
      vector magnitude squared (or mass squared)
      In case of negative mass (spacelike particles return negative values)
   */
   Scalar M2() const { return (fM >= 0) ? fM * fM : -fM * fM; }
   Scalar Mag2() const { return M2(); }
 
   Scalar Mag() const { return M(); }
 
   /**
      energy squared
   */
   Scalar E2() const
   {
      Scalar e2 = P2() + M2();
      // protect against numerical errors when M2() is negative
      return e2 > 0 ? e2 : 0;
   }
 
   /**
       transverse spatial component squared
   */
   Scalar Pt2() const { return fX * fX + fY * fY; }
   Scalar Perp2() const { return Pt2(); }
 
   /**
      Transverse spatial component (P_perp or rho)
   */
   Scalar Pt() const { return math_sqrt(Perp2()); }
   Scalar Perp() const { return Pt(); }
   Scalar Rho() const { return Pt(); }
 
   /**
       transverse mass squared
   */
   Scalar Mt2() const { return E2() - fZ * fZ; }
 
   /**
      transverse mass
   */
   Scalar Mt() const
   {
      const Scalar mm = Mt2();
      if (mm >= 0) {
         return math_sqrt(mm);
      } else {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
         GenVector_Throw("PxPyPzM4D::Mt() - Tachyonic:\n"
                         "    Pz^2 > E^2 so the transverse mass would be imaginary");
#endif
         return -math_sqrt(-mm);
      }
   }
 
   /**
       transverse energy squared
   */
   Scalar Et2() const
   { // is (E^2 * pt ^2) / p^2
      // but it is faster to form p^2 from pt^2
      Scalar pt2 = Pt2();
      return pt2 == 0 ? 0 : E2() * pt2 / (pt2 + fZ * fZ);
   }
 
   /**
      transverse energy
   */
   Scalar Et() const
   {
      const Scalar etet = Et2();
      return math_sqrt(etet);
   }
 
   /**
      azimuthal angle
   */
   Scalar Phi() const { return (fX == 0.0 && fY == 0.0) ? 0.0 : math_atan2(fY, fX); }
 
   /**
      polar angle
   */
   Scalar Theta() const { return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : math_atan2(Pt(), fZ); }
 
   /**
       pseudorapidity
   */
   Scalar Eta() const { return Impl::Eta_FromRhoZ(Pt(), fZ); }
 
   // --------- Set Coordinates of this system  ---------------
 
   /**
      set X value
   */
   void SetPx(Scalar px) { fX = px; }
   /**
      set Y value
   */
   void SetPy(Scalar py) { fY = py; }
   /**
      set Z value
   */
   void SetPz(Scalar pz) { fZ = pz; }
   /**
      set T value
   */
   void SetM(Scalar m)
   {
      fM = m;
      if (fM < 0)
         RestrictNegMass();
   }
 
   /**
       set all values
   */
   void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
 
   // ------ Manipulations -------------
 
   /**
      negate the 4-vector -  Note that the energy cannot be negate (would need an additional data member)
      therefore negate will work only on the spatial components.
      One would need to use negate only with vectors having the energy as data members
   */
   void Negate()
   {
      fX = -fX;
      fY = -fY;
      fZ = -fZ;
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("PxPyPzM4D::Negate - cannot negate the energy - can negate only the spatial components");
#endif
   }
 
   /**
      scale coordinate values by a scalar quantity a
   */
   void Scale(const Scalar &a)
   {
      fX *= a;
      fY *= a;
      fZ *= a;
      fM *= a;
   }
 
   /**
      Exact equality
   */
   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)); }
 
   // ============= Compatibility section ==================
 
   // The following make this coordinate system look enough like a CLHEP
   // vector that an assignment member template can work with either
   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)
 
   // ====== Set member functions for coordinates in other systems =======
 
   void SetPt(Scalar pt);
 
   void SetEta(Scalar eta);
 
   void SetPhi(Scalar phi);
 
   void SetE(Scalar t);
 
#endif
 
private:
   // restrict the value of negative mass to avoid unphysical negative E2 values
   // M2 must be less than P2 for the tachionic particles - otherwise use positive values
   inline void RestrictNegMass()
   {
      if (fM >= 0)
         return;
      if (P2() - fM * fM < 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
         GenVector_Throw("PxPyPzM4D::unphysical value of mass, set to closest physical value");
#endif
         fM = -P();
      }
      return;
   }
 
   /**
      (contiguous) data containing the coordinate values x,y,z,t
   */
 
   ScalarType fX = 0;
   ScalarType fY = 0;
   ScalarType fZ = 0;
   ScalarType fM = 0;
};
 
} // end namespace ROOT_MATH_ARCH
} // end namespace ROOT
 
// move implementations here to avoid circle dependencies
 
#include "MathX/GenVectorX/PxPyPzE4D.h"
#include "MathX/GenVectorX/PtEtaPhiM4D.h"
 
namespace ROOT {
 
namespace ROOT_MATH_ARCH {
 
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)
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
 
// ====== Set member functions for coordinates in other systems =======
 
// ====== Set member functions for coordinates in other systems =======
 
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
 
} // end namespace ROOT_MATH_ARCH
 
} // end namespace ROOT
 
#endif // ROOT_MathX_GenVectorX_PxPyPzM4D