Math/GenVector/PxPyPzM4D.h

0001 // @(#)root/mathcore:$Id: 464c29f33a8bbd8462a3e15b7e4c30c6f5b74a30 $
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003
0004 /**********************************************************************
0005 *                                                                    *
0006 * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0007 *                                                                    *
0008 *                                                                    *
0009 **********************************************************************/
0010
0011 // Header file for class PxPyPzM4D
0012 //
0013 // Created by: fischler at Wed Jul 20   2005
0014 //   (starting from PxPyPzM4D by moneta)
0015 //
0016 // Last update: $Id: 464c29f33a8bbd8462a3e15b7e4c30c6f5b74a30 $
0017 //
0018 #ifndef ROOT_Math_GenVector_PxPyPzM4D
0019 #define ROOT_Math_GenVector_PxPyPzM4D  1
0020
0021 #include "Math/GenVector/eta.h"
0022
0023 #include "Math/GenVector/GenVector_exception.h"
0024
0025
0026 #include <cmath>
0027
0028 namespace ROOT {
0029
0030 namespace Math {
0031
0032 //__________________________________________________________________________________________
0033 /**
0034     Class describing a 4D coordinate system
0035     or momentum-energy vectors stored as (Px, Py, Pz, M).
0036     This system is useful to describe ultra-relativistic particles
0037     (like electrons at LHC) to avoid numerical errors evaluating the mass
0038     when E >>> m
0039     The metric used is (-,-,-,+)
0040     Spacelike particles (M2 < 0) are described with negative mass values,
0041     but in this case m2 must always be less than P2 to preserve a positive value of E2
0042
0043     @ingroup GenVector
0044
0045     @see GenVector
0046 */
0047
0048 template <class ScalarType = double>
0049 class PxPyPzM4D {
0050
0051 public :
0052
0053    typedef ScalarType Scalar;
0054    static constexpr unsigned int Dimension = 4U;
0055
0056    // --------- Constructors ---------------
0057
0058    /**
0059       Default constructor  with x=y=z=m=0
0060    */
0061    constexpr PxPyPzM4D() noexcept = default;
0062
0063    /**
0064       Constructor  from x, y , z , m values
0065    */
0066    PxPyPzM4D(Scalar px, Scalar py, Scalar pz, Scalar m) :
0067       fX(px), fY(py), fZ(pz), fM(m) {
0068
0069       if (fM < 0) RestrictNegMass();
0070    }
0071
0072    /**
0073       construct from any 4D  coordinate system class
0074       implementing X(), Y(), X() and M()
0075    */
0076    template <class CoordSystem>
0077    explicit constexpr PxPyPzM4D(const CoordSystem & v) :
0078       fX( v.X() ), fY( v.Y() ), fZ( v.Z() ), fM( v.M() )
0079    { }
0080
0081    /**
0082       construct from any 4D  coordinate system class
0083       implementing X(), Y(), X() and M()
0084    */
0085    template <class AnyCoordSystem>
0086    PxPyPzM4D & operator = (const AnyCoordSystem & v) {
0087       fX = v.X();
0088       fY = v.Y();
0089       fZ = v.Z();
0090       fM = v.M();
0091       return *this;
0092    }
0093
0094    /**
0095       Set internal data based on an array of 4 Scalar numbers
0096    */
0097    void SetCoordinates( const Scalar src[] ) {
0098       fX=src[0]; fY=src[1]; fZ=src[2]; fM=src[3];
0099       if (fM < 0) RestrictNegMass();
0100    }
0101
0102    /**
0103       get internal data into an array of 4 Scalar numbers
0104    */
0105    void GetCoordinates( Scalar dest[] ) const
0106    { dest[0] = fX; dest[1] = fY; dest[2] = fZ; dest[3] = fM; }
0107
0108    /**
0109       Set internal data based on 4 Scalar numbers
0110    */
0111    void SetCoordinates(Scalar  px, Scalar  py, Scalar  pz, Scalar m) {
0112       fX=px; fY=py; fZ=pz; fM=m;
0113       if (fM < 0) RestrictNegMass();
0114    }
0115
0116    /**
0117       get internal data into 4 Scalar numbers
0118    */
0119    void GetCoordinates(Scalar& px, Scalar& py, Scalar& pz, Scalar& m) const
0120    { px=fX; py=fY; pz=fZ; m=fM;}
0121
0122    // --------- Coordinates and Coordinate-like Scalar properties -------------
0123
0124    // cartesian (Minkowski)coordinate accessors
0125
0126    Scalar Px() const { return fX;}
0127    Scalar Py() const { return fY;}
0128    Scalar Pz() const { return fZ;}
0129    Scalar M() const  { return fM; }
0130
0131    Scalar X() const { return fX;}
0132    Scalar Y() const { return fY;}
0133    Scalar Z() const { return fZ;}
0134
0135    // other coordinate representation
0136    /**
0137       Energy
0138     */
0139    Scalar E() const { using std::sqrt; return sqrt(E2()); }
0140
0141    Scalar T() const { return E();}
0142
0143    /**
0144       squared magnitude of spatial components
0145    */
0146    Scalar P2() const { return fX*fX + fY*fY + fZ*fZ; }
0147
0148    /**
0149       magnitude of spatial components (magnitude of 3-momentum)
0150    */
0151    Scalar P() const { using std::sqrt; return sqrt(P2()); }
0152    Scalar R() const { return P(); }
0153
0154    /**
0155       vector magnitude squared (or mass squared)
0156       In case of negative mass (spacelike particles return negative values)
0157    */
0158    Scalar M2() const   {
0159       return ( fM  >= 0 ) ?  fM*fM :  -fM*fM;
0160    }
0161    Scalar Mag2() const { return M2(); }
0162
0163    Scalar Mag() const    { return M(); }
0164
0165    /**
0166       energy squared
0167    */
0168    Scalar E2() const {
0169       Scalar e2 =  P2() + M2();
0170       // protect against numerical errors when M2() is negative
0171       return e2 > 0 ? e2 : 0;
0172    }
0173
0174    /**
0175        transverse spatial component squared
0176    */
0177    Scalar Pt2()   const { return fX*fX + fY*fY;}
0178    Scalar Perp2() const { return Pt2();}
0179
0180    /**
0181       Transverse spatial component (P_perp or rho)
0182    */
0183    Scalar Pt() const { using std::sqrt; return sqrt(Perp2()); }
0184    Scalar Perp() const { return Pt();}
0185    Scalar Rho()  const { return Pt();}
0186
0187    /**
0188        transverse mass squared
0189    */
0190    Scalar Mt2() const { return E2() - fZ*fZ; }
0191
0192    /**
0193       transverse mass
0194    */
0195    Scalar Mt() const {
0196       const Scalar mm = Mt2();
0197       if (mm >= 0) {
0198          using std::sqrt;
0199          return sqrt(mm);
0200       } else {
0201          GenVector_Throw("PxPyPzM4D::Mt() - Tachyonic:\n"
0202                          "    Pz^2 > E^2 so the transverse mass would be imaginary");
0203          using std::sqrt;
0204          return -sqrt(-mm);
0205       }
0206    }
0207
0208    /**
0209        transverse energy squared
0210    */
0211    Scalar Et2() const {  // is (E^2 * pt ^2) / p^2
0212       // but it is faster to form p^2 from pt^2
0213       Scalar pt2 = Pt2();
0214       return pt2 == 0 ? 0 : E2() * pt2/( pt2 + fZ*fZ );
0215    }
0216
0217    /**
0218       transverse energy
0219    */
0220    Scalar Et() const {
0221       const Scalar etet = Et2();
0222       using std::sqrt;
0223       return sqrt(etet);
0224    }
0225
0226    /**
0227       azimuthal angle
0228    */
0229    Scalar Phi() const { using std::atan2; return (fX == 0.0 && fY == 0.0) ? 0.0 : atan2(fY, fX); }
0230
0231    /**
0232       polar angle
0233    */
0234    Scalar Theta() const { using std::atan2; return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : atan2(Pt(), fZ); }
0235
0236    /**
0237        pseudorapidity
0238    */
0239    Scalar Eta() const {
0240       return Impl::Eta_FromRhoZ ( Pt(), fZ);
0241    }
0242
0243    // --------- Set Coordinates of this system  ---------------
0244
0245
0246    /**
0247       set X value
0248    */
0249    void SetPx( Scalar  px) {
0250       fX = px;
0251    }
0252    /**
0253       set Y value
0254    */
0255    void SetPy( Scalar  py) {
0256       fY = py;
0257    }
0258    /**
0259       set Z value
0260    */
0261    void SetPz( Scalar  pz) {
0262       fZ = pz;
0263    }
0264    /**
0265       set T value
0266    */
0267    void SetM( Scalar  m) {
0268       fM = m;
0269       if (fM < 0) RestrictNegMass();
0270    }
0271
0272    /**
0273        set all values
0274    */
0275    void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
0276
0277    // ------ Manipulations -------------
0278
0279    /**
0280       negate the 4-vector -  Note that the energy cannot be negate (would need an additional data member)
0281       therefore negate will work only on the spatial components.
0282       One would need to use negate only with vectors having the energy as data members
0283    */
0284    void Negate( ) {
0285       fX = -fX;
0286       fY = -fY;
0287       fZ = -fZ;
0288       GenVector_Throw("PxPyPzM4D::Negate - cannot negate the energy - can negate only the spatial components");
0289    }
0290
0291    /**
0292       scale coordinate values by a scalar quantity a
0293    */
0294    void Scale( const Scalar & a) {
0295       fX *= a;
0296       fY *= a;
0297       fZ *= a;
0298       fM *= a;
0299    }
0300
0301
0302    /**
0303       Exact equality
0304    */
0305    bool operator == (const PxPyPzM4D & rhs) const {
0306       return fX == rhs.fX && fY == rhs.fY && fZ == rhs.fZ && fM == rhs.fM;
0307    }
0308    bool operator != (const PxPyPzM4D & rhs) const {return !(operator==(rhs));}
0309
0310
0311    // ============= Compatibility section ==================
0312
0313    // The following make this coordinate system look enough like a CLHEP
0314    // vector that an assignment member template can work with either
0315    Scalar x() const { return X(); }
0316    Scalar y() const { return Y(); }
0317    Scalar z() const { return Z(); }
0318    Scalar t() const { return E(); }
0319
0320
0321
0322 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0323
0324    // ====== Set member functions for coordinates in other systems =======
0325
0326    void SetPt(Scalar pt);
0327
0328    void SetEta(Scalar eta);
0329
0330    void SetPhi(Scalar phi);
0331
0332    void SetE(Scalar t);
0333
0334 #endif
0335
0336 private:
0337
0338    // restrict the value of negative mass to avoid unphysical negative E2 values
0339    // M2 must be less than P2 for the tachionic particles - otherwise use positive values
0340    inline void RestrictNegMass() {
0341       if ( fM >=0 ) return;
0342       if ( P2() - fM*fM  < 0 ) {
0343          GenVector_Throw("PxPyPzM4D::unphysical value of mass, set to closest physical value");
0344          fM = - P();
0345       }
0346       return;
0347    }
0348
0349
0350    /**
0351       (contiguous) data containing the coordinate values x,y,z,t
0352    */
0353
0354    ScalarType fX = 0;
0355    ScalarType fY = 0;
0356    ScalarType fZ = 0;
0357    ScalarType fM = 0;
0358 };
0359
0360 } // end namespace Math
0361 } // end namespace ROOT
0362
0363
0364 // move implementations here to avoid circle dependencies
0365
0366 #include "Math/GenVector/PxPyPzE4D.h"
0367 #include "Math/GenVector/PtEtaPhiM4D.h"
0368
0369 namespace ROOT {
0370
0371 namespace Math {
0372
0373 template <class ScalarType>
0374 inline void PxPyPzM4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
0375    *this = PxPyPzE4D<Scalar> (px, py, pz, e);
0376 }
0377
0378
0379 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0380
0381   // ====== Set member functions for coordinates in other systems =======
0382
0383   // ====== Set member functions for coordinates in other systems =======
0384
0385 template <class ScalarType>
0386 inline void PxPyPzM4D<ScalarType>::SetPt(ScalarType pt) {
0387    GenVector_exception e("PxPyPzM4D::SetPt() is not supposed to be called");
0388    throw e;
0389    PtEtaPhiE4D<ScalarType> v(*this); v.SetPt(pt); *this = PxPyPzM4D<ScalarType>(v);
0390 }
0391 template <class ScalarType>
0392 inline void PxPyPzM4D<ScalarType>::SetEta(ScalarType eta) {
0393    GenVector_exception e("PxPyPzM4D::SetEta() is not supposed to be called");
0394    throw e;
0395    PtEtaPhiE4D<ScalarType> v(*this); v.SetEta(eta); *this = PxPyPzM4D<ScalarType>(v);
0396 }
0397 template <class ScalarType>
0398 inline void PxPyPzM4D<ScalarType>::SetPhi(ScalarType phi) {
0399    GenVector_exception e("PxPyPzM4D::SetPhi() is not supposed to be called");
0400    throw e;
0401    PtEtaPhiE4D<ScalarType> v(*this); v.SetPhi(phi); *this = PxPyPzM4D<ScalarType>(v);
0402 }
0403 template <class ScalarType>
0404 inline void PxPyPzM4D<ScalarType>::SetE(ScalarType energy) {
0405    GenVector_exception e("PxPyPzM4D::SetE() is not supposed to be called");
0406    throw e;
0407    PxPyPzE4D<ScalarType> v(*this); v.SetE(energy);
0408    *this = PxPyPzM4D<ScalarType>(v);
0409 }
0410
0411
0412 #endif  // endif __MAKE__CINT || G__DICTIONARY
0413
0414 } // end namespace Math
0415
0416 } // end namespace ROOT
0417
0418
0419
0420 #endif // ROOT_Math_GenVector_PxPyPzM4D