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     @sa Overview of the @ref GenVector "physics vector library"
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    PxPyPzM4D() : fX(0.0), fY(0.0), fZ(0.0), fM(0.0) { }
0062
0063
0064    /**
0065       Constructor  from x, y , z , m values
0066    */
0067    PxPyPzM4D(Scalar px, Scalar py, Scalar pz, Scalar m) :
0068       fX(px), fY(py), fZ(pz), fM(m) {
0069
0070       if (fM < 0) RestrictNegMass();
0071    }
0072
0073    /**
0074       construct from any 4D  coordinate system class
0075       implementing X(), Y(), X() and M()
0076    */
0077    template <class CoordSystem>
0078    explicit constexpr PxPyPzM4D(const CoordSystem & v) :
0079       fX( v.X() ), fY( v.Y() ), fZ( v.Z() ), fM( v.M() )
0080    { }
0081
0082    // for g++  3.2 and 3.4 on 32 bits found that the compiler generated copy ctor and assignment are much slower
0083    // so we decided to re-implement them ( there is no no need to have them with g++4)
0084    /**
0085       copy constructor
0086     */
0087    PxPyPzM4D(const PxPyPzM4D & v) :
0088       fX(v.fX), fY(v.fY), fZ(v.fZ), fM(v.fM) { }
0089
0090    /**
0091       assignment operator
0092     */
0093    PxPyPzM4D & operator = (const PxPyPzM4D & v) {
0094       fX = v.fX;
0095       fY = v.fY;
0096       fZ = v.fZ;
0097       fM = v.fM;
0098       return *this;
0099    }
0100
0101
0102    /**
0103       construct from any 4D  coordinate system class
0104       implementing X(), Y(), X() and M()
0105    */
0106    template <class AnyCoordSystem>
0107    PxPyPzM4D & operator = (const AnyCoordSystem & v) {
0108       fX = v.X();
0109       fY = v.Y();
0110       fZ = v.Z();
0111       fM = v.M();
0112       return *this;
0113    }
0114
0115    /**
0116       Set internal data based on an array of 4 Scalar numbers
0117    */
0118    void SetCoordinates( const Scalar src[] ) {
0119       fX=src[0]; fY=src[1]; fZ=src[2]; fM=src[3];
0120       if (fM < 0) RestrictNegMass();
0121    }
0122
0123    /**
0124       get internal data into an array of 4 Scalar numbers
0125    */
0126    void GetCoordinates( Scalar dest[] ) const
0127    { dest[0] = fX; dest[1] = fY; dest[2] = fZ; dest[3] = fM; }
0128
0129    /**
0130       Set internal data based on 4 Scalar numbers
0131    */
0132    void SetCoordinates(Scalar  px, Scalar  py, Scalar  pz, Scalar m) {
0133       fX=px; fY=py; fZ=pz; fM=m;
0134       if (fM < 0) RestrictNegMass();
0135    }
0136
0137    /**
0138       get internal data into 4 Scalar numbers
0139    */
0140    void GetCoordinates(Scalar& px, Scalar& py, Scalar& pz, Scalar& m) const
0141    { px=fX; py=fY; pz=fZ; m=fM;}
0142
0143    // --------- Coordinates and Coordinate-like Scalar properties -------------
0144
0145    // cartesian (Minkowski)coordinate accessors
0146
0147    Scalar Px() const { return fX;}
0148    Scalar Py() const { return fY;}
0149    Scalar Pz() const { return fZ;}
0150    Scalar M() const  { return fM; }
0151
0152    Scalar X() const { return fX;}
0153    Scalar Y() const { return fY;}
0154    Scalar Z() const { return fZ;}
0155
0156    // other coordinate representation
0157    /**
0158       Energy
0159     */
0160    Scalar E() const { using std::sqrt; return sqrt(E2()); }
0161
0162    Scalar T() const { return E();}
0163
0164    /**
0165       squared magnitude of spatial components
0166    */
0167    Scalar P2() const { return fX*fX + fY*fY + fZ*fZ; }
0168
0169    /**
0170       magnitude of spatial components (magnitude of 3-momentum)
0171    */
0172    Scalar P() const { using std::sqrt; return sqrt(P2()); }
0173    Scalar R() const { return P(); }
0174
0175    /**
0176       vector magnitude squared (or mass squared)
0177       In case of negative mass (spacelike particles return negative values)
0178    */
0179    Scalar M2() const   {
0180       return ( fM  >= 0 ) ?  fM*fM :  -fM*fM;
0181    }
0182    Scalar Mag2() const { return M2(); }
0183
0184    Scalar Mag() const    { return M(); }
0185
0186    /**
0187       energy squared
0188    */
0189    Scalar E2() const {
0190       Scalar e2 =  P2() + M2();
0191       // protect against numerical errors when M2() is negative
0192       return e2 > 0 ? e2 : 0;
0193    }
0194
0195    /**
0196        transverse spatial component squared
0197    */
0198    Scalar Pt2()   const { return fX*fX + fY*fY;}
0199    Scalar Perp2() const { return Pt2();}
0200
0201    /**
0202       Transverse spatial component (P_perp or rho)
0203    */
0204    Scalar Pt() const { using std::sqrt; return sqrt(Perp2()); }
0205    Scalar Perp() const { return Pt();}
0206    Scalar Rho()  const { return Pt();}
0207
0208    /**
0209        transverse mass squared
0210    */
0211    Scalar Mt2() const { return E2() - fZ*fZ; }
0212
0213    /**
0214       transverse mass
0215    */
0216    Scalar Mt() const {
0217       const Scalar mm = Mt2();
0218       if (mm >= 0) {
0219          using std::sqrt;
0220          return sqrt(mm);
0221       } else {
0222          GenVector::Throw ("PxPyPzM4D::Mt() - Tachyonic:\n"
0223                            "    Pz^2 > E^2 so the transverse mass would be imaginary");
0224          using std::sqrt;
0225          return -sqrt(-mm);
0226       }
0227    }
0228
0229    /**
0230        transverse energy squared
0231    */
0232    Scalar Et2() const {  // is (E^2 * pt ^2) / p^2
0233       // but it is faster to form p^2 from pt^2
0234       Scalar pt2 = Pt2();
0235       return pt2 == 0 ? 0 : E2() * pt2/( pt2 + fZ*fZ );
0236    }
0237
0238    /**
0239       transverse energy
0240    */
0241    Scalar Et() const {
0242       const Scalar etet = Et2();
0243       using std::sqrt;
0244       return sqrt(etet);
0245    }
0246
0247    /**
0248       azimuthal angle
0249    */
0250    Scalar Phi() const { using std::atan2; return (fX == 0.0 && fY == 0.0) ? 0.0 : atan2(fY, fX); }
0251
0252    /**
0253       polar angle
0254    */
0255    Scalar Theta() const { using std::atan2; return (fX == 0.0 && fY == 0.0 && fZ == 0.0) ? 0 : atan2(Pt(), fZ); }
0256
0257    /**
0258        pseudorapidity
0259    */
0260    Scalar Eta() const {
0261       return Impl::Eta_FromRhoZ ( Pt(), fZ);
0262    }
0263
0264    // --------- Set Coordinates of this system  ---------------
0265
0266
0267    /**
0268       set X value
0269    */
0270    void SetPx( Scalar  px) {
0271       fX = px;
0272    }
0273    /**
0274       set Y value
0275    */
0276    void SetPy( Scalar  py) {
0277       fY = py;
0278    }
0279    /**
0280       set Z value
0281    */
0282    void SetPz( Scalar  pz) {
0283       fZ = pz;
0284    }
0285    /**
0286       set T value
0287    */
0288    void SetM( Scalar  m) {
0289       fM = m;
0290       if (fM < 0) RestrictNegMass();
0291    }
0292
0293    /**
0294        set all values
0295    */
0296    void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
0297
0298    // ------ Manipulations -------------
0299
0300    /**
0301       negate the 4-vector -  Note that the energy cannot be negate (would need an additional data member)
0302       therefore negate will work only on the spatial components.
0303       One would need to use negate only with vectors having the energy as data members
0304    */
0305    void Negate( ) {
0306       fX = -fX;
0307       fY = -fY;
0308       fZ = -fZ;
0309       GenVector::Throw ("PxPyPzM4D::Negate - cannot negate the energy - can negate only the spatial components");
0310    }
0311
0312    /**
0313       scale coordinate values by a scalar quantity a
0314    */
0315    void Scale( const Scalar & a) {
0316       fX *= a;
0317       fY *= a;
0318       fZ *= a;
0319       fM *= a;
0320    }
0321
0322
0323    /**
0324       Exact equality
0325    */
0326    bool operator == (const PxPyPzM4D & rhs) const {
0327       return fX == rhs.fX && fY == rhs.fY && fZ == rhs.fZ && fM == rhs.fM;
0328    }
0329    bool operator != (const PxPyPzM4D & rhs) const {return !(operator==(rhs));}
0330
0331
0332    // ============= Compatibility section ==================
0333
0334    // The following make this coordinate system look enough like a CLHEP
0335    // vector that an assignment member template can work with either
0336    Scalar x() const { return X(); }
0337    Scalar y() const { return Y(); }
0338    Scalar z() const { return Z(); }
0339    Scalar t() const { return E(); }
0340
0341
0342
0343 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0344
0345    // ====== Set member functions for coordinates in other systems =======
0346
0347    void SetPt(Scalar pt);
0348
0349    void SetEta(Scalar eta);
0350
0351    void SetPhi(Scalar phi);
0352
0353    void SetE(Scalar t);
0354
0355 #endif
0356
0357 private:
0358
0359    // restrict the value of negative mass to avoid unphysical negative E2 values
0360    // M2 must be less than P2 for the tachionic particles - otherwise use positive values
0361    inline void RestrictNegMass() {
0362       if ( fM >=0 ) return;
0363       if ( P2() - fM*fM  < 0 ) {
0364          GenVector::Throw("PxPyPzM4D::unphysical value of mass, set to closest physical value");
0365          fM = - P();
0366       }
0367       return;
0368    }
0369
0370
0371    /**
0372       (contiguous) data containing the coordinate values x,y,z,t
0373    */
0374
0375    ScalarType fX;
0376    ScalarType fY;
0377    ScalarType fZ;
0378    ScalarType fM;
0379
0380 };
0381
0382 } // end namespace Math
0383 } // end namespace ROOT
0384
0385
0386 // move implementations here to avoid circle dependencies
0387
0388 #include "Math/GenVector/PxPyPzE4D.h"
0389 #include "Math/GenVector/PtEtaPhiM4D.h"
0390
0391 namespace ROOT {
0392
0393 namespace Math {
0394
0395 template <class ScalarType>
0396 inline void PxPyPzM4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
0397    *this = PxPyPzE4D<Scalar> (px, py, pz, e);
0398 }
0399
0400
0401 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0402
0403   // ====== Set member functions for coordinates in other systems =======
0404
0405   // ====== Set member functions for coordinates in other systems =======
0406
0407 template <class ScalarType>
0408 inline void PxPyPzM4D<ScalarType>::SetPt(ScalarType pt) {
0409    GenVector_exception e("PxPyPzM4D::SetPt() is not supposed to be called");
0410    throw e;
0411    PtEtaPhiE4D<ScalarType> v(*this); v.SetPt(pt); *this = PxPyPzM4D<ScalarType>(v);
0412 }
0413 template <class ScalarType>
0414 inline void PxPyPzM4D<ScalarType>::SetEta(ScalarType eta) {
0415    GenVector_exception e("PxPyPzM4D::SetEta() is not supposed to be called");
0416    throw e;
0417    PtEtaPhiE4D<ScalarType> v(*this); v.SetEta(eta); *this = PxPyPzM4D<ScalarType>(v);
0418 }
0419 template <class ScalarType>
0420 inline void PxPyPzM4D<ScalarType>::SetPhi(ScalarType phi) {
0421    GenVector_exception e("PxPyPzM4D::SetPhi() is not supposed to be called");
0422    throw e;
0423    PtEtaPhiE4D<ScalarType> v(*this); v.SetPhi(phi); *this = PxPyPzM4D<ScalarType>(v);
0424 }
0425 template <class ScalarType>
0426 inline void PxPyPzM4D<ScalarType>::SetE(ScalarType energy) {
0427    GenVector_exception e("PxPyPzM4D::SetE() is not supposed to be called");
0428    throw e;
0429    PxPyPzE4D<ScalarType> v(*this); v.SetE(energy);
0430    *this = PxPyPzM4D<ScalarType>(v);
0431 }
0432
0433
0434 #endif  // endif __MAKE__CINT || G__DICTIONARY
0435
0436 } // end namespace Math
0437
0438 } // end namespace ROOT
0439
0440
0441
0442 #endif // ROOT_Math_GenVector_PxPyPzM4D