Math/GenVector/LorentzVector.h

0001 // @(#)root/mathcore:$Id$
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 LorentzVector
0012 //
0013 // Created by:    moneta   at Tue May 31 17:06:09 2005
0014 // Major mods by: fischler at Wed Jul 20   2005
0015 //
0016 // Last update: $Id$
0017 //
0018 #ifndef ROOT_Math_GenVector_LorentzVector
0019 #define ROOT_Math_GenVector_LorentzVector  1
0020
0021 #include "Math/GenVector/PxPyPzE4D.h"
0022 #include "Math/GenVector/DisplacementVector3D.h"
0023 #include "Math/GenVector/GenVectorIO.h"
0024 #include "Math/Vector2D.h"
0025
0026 #include "TMath.h"
0027
0028 #include <cmath>
0029 #include <string>
0030
0031 namespace ROOT {
0032
0033   namespace Math {
0034
0035 //__________________________________________________________________________________________
0036 /** @ingroup GenVector
0037
0038 Class describing a generic LorentzVector in the 4D space-time,
0039 using the specified coordinate system for the spatial vector part.
0040 The metric used for the LorentzVector is (-,-,-,+).
0041 In the case of LorentzVector we don't distinguish the concepts
0042 of points and displacement vectors as in the 3D case,
0043 since the main use case for 4D Vectors is to describe the kinematics of
0044 relativistic particles. A LorentzVector behaves like a
0045 DisplacementVector in 4D. The Minkowski components could be viewed as
0046 v and t, or for kinematic 4-vectors, as p and E.
0047
0048 ROOT provides specialisations (and aliases to them) of the ROOT::Math::LorentzVector template:
0049 - ROOT::Math::PtEtaPhiMVector based on pt (rho),eta,phi and M (t) coordinates in double precision
0050 - ROOT::Math::PtEtaPhiEVector based on pt (rho),eta,phi and E (t) coordinates in double precision
0051 - ROOT::Math::PxPyPzMVector based on px,py,pz and M coordinates in double precision
0052 - ROOT::Math::PxPyPzEVector based on px,py,pz and E coordinates in double precision
0053 - ROOT::Math::XYZTVector based on x,y,z,t coordinates (cartesian) in double precision (same as PxPyPzEVector)
0054 - ROOT::Math::XYZTVectorF based on x,y,z,t coordinates (cartesian) in float precision (same as PxPyPzEVector but float)
0055 Pt (or rho) refers to transverse momentum, whereas eta refers to pseudorapidity. M is mass, E is energy, p is momentum.
0056
0057 @see GenVector
0058 */
0059
0060     template< class CoordSystem >
0061     class LorentzVector {
0062
0063     public:
0064
0065        // ------ ctors ------
0066
0067        typedef typename CoordSystem::Scalar Scalar;
0068        typedef CoordSystem CoordinateType;
0069
0070        /**
0071           default constructor of an empty vector (Px = Py = Pz = E = 0 )
0072        */
0073        LorentzVector ( ) : fCoordinates() { }
0074
0075        /**
0076           generic constructors from four scalar values.
0077           The association between values and coordinate depends on the
0078           coordinate system.  For PxPyPzE4D,
0079           \param a scalar value (Px)
0080           \param b scalar value (Py)
0081           \param c scalar value (Pz)
0082           \param d scalar value (E)
0083        */
0084        LorentzVector(const Scalar & a,
0085                      const Scalar & b,
0086                      const Scalar & c,
0087                      const Scalar & d) :
0088           fCoordinates(a , b,  c, d)  { }
0089
0090        /**
0091           constructor from a LorentzVector expressed in different
0092           coordinates, or using a different Scalar type
0093        */
0094        template< class Coords >
0095        explicit constexpr LorentzVector(const LorentzVector<Coords> & v ) :
0096           fCoordinates( v.Coordinates() ) { }
0097
0098        /**
0099           Construct from a foreign 4D vector type, for example, HepLorentzVector
0100           Precondition: v must implement methods x(), y(), z(), and t()
0101        */
0102        template<class ForeignLorentzVector,
0103                 typename = decltype(std::declval<ForeignLorentzVector>().x()
0104                                     + std::declval<ForeignLorentzVector>().y()
0105                                     + std::declval<ForeignLorentzVector>().z()
0106                                     + std::declval<ForeignLorentzVector>().t())>
0107        explicit constexpr LorentzVector( const ForeignLorentzVector & v) :
0108           fCoordinates(PxPyPzE4D<Scalar>( v.x(), v.y(), v.z(), v.t()  ) ) { }
0109
0110 #ifdef LATER
0111        /**
0112           construct from a generic linear algebra  vector implementing operator []
0113           and with a size of at least 4. This could be also a C array
0114           In this case v[0] is the first data member
0115           ( Px for a PxPyPzE4D base)
0116           \param v LA vector
0117           \param index0 index of first vector element (Px)
0118        */
0119        template< class LAVector >
0120        explicit constexpr LorentzVector(const LAVector & v, size_t index0 ) {
0121           fCoordinates = CoordSystem ( v[index0], v[index0+1], v[index0+2], v[index0+3] );
0122        }
0123 #endif
0124
0125
0126        // ------ assignment ------
0127
0128        /**
0129           Assignment operator from a lorentz vector of arbitrary type
0130        */
0131        template< class OtherCoords >
0132        LorentzVector & operator= ( const LorentzVector<OtherCoords> & v) {
0133           fCoordinates = v.Coordinates();
0134           return *this;
0135        }
0136
0137        /**
0138           assignment from any other Lorentz vector  implementing
0139           x(), y(), z() and t()
0140        */
0141        template<class ForeignLorentzVector,
0142                 typename = decltype(std::declval<ForeignLorentzVector>().x()
0143                                     + std::declval<ForeignLorentzVector>().y()
0144                                     + std::declval<ForeignLorentzVector>().z()
0145                                     + std::declval<ForeignLorentzVector>().t())>
0146        LorentzVector & operator = ( const ForeignLorentzVector & v) {
0147           SetXYZT( v.x(), v.y(), v.z(), v.t() );
0148           return *this;
0149        }
0150
0151 #ifdef LATER
0152        /**
0153           assign from a generic linear algebra  vector implementing operator []
0154           and with a size of at least 4
0155           In this case v[0] is the first data member
0156           ( Px for a PxPyPzE4D base)
0157           \param v LA vector
0158           \param index0 index of first vector element (Px)
0159        */
0160        template< class LAVector >
0161        LorentzVector & AssignFrom(const LAVector & v, size_t index0=0 ) {
0162           fCoordinates.SetCoordinates( v[index0], v[index0+1], v[index0+2], v[index0+3] );
0163           return *this;
0164        }
0165 #endif
0166
0167        // ------ Set, Get, and access coordinate data ------
0168
0169        /**
0170           Retrieve a const reference to  the coordinates object
0171        */
0172        const CoordSystem & Coordinates() const {
0173           return fCoordinates;
0174        }
0175
0176        /**
0177           Set internal data based on an array of 4 Scalar numbers
0178        */
0179        LorentzVector<CoordSystem>& SetCoordinates( const Scalar src[] ) {
0180           fCoordinates.SetCoordinates(src);
0181           return *this;
0182        }
0183
0184        /**
0185           Set internal data based on 4 Scalar numbers
0186        */
0187        LorentzVector<CoordSystem>& SetCoordinates( Scalar a, Scalar b, Scalar c, Scalar d ) {
0188           fCoordinates.SetCoordinates(a, b, c, d);
0189           return *this;
0190        }
0191
0192        /**
0193           Set internal data based on 4 Scalars at *begin to *end
0194        */
0195        template< class IT >
0196        LorentzVector<CoordSystem>& SetCoordinates( IT begin, IT end  ) {
0197           IT a = begin; IT b = ++begin; IT c = ++begin; IT d = ++begin;
0198           (void)end;
0199           assert (++begin==end);
0200           SetCoordinates (*a,*b,*c,*d);
0201           return *this;
0202        }
0203
0204        /**
0205           get internal data into 4 Scalar numbers
0206           \note Alternatively, you may use structured bindings: `auto const [a, b, c, d] = v`.
0207        */
0208        void GetCoordinates( Scalar& a, Scalar& b, Scalar& c, Scalar & d ) const
0209        { fCoordinates.GetCoordinates(a, b, c, d);  }
0210
0211        /**
0212           get internal data into an array of 4 Scalar numbers
0213        */
0214        void GetCoordinates( Scalar dest[] ) const
0215        { fCoordinates.GetCoordinates(dest);  }
0216
0217        /**
0218           get internal data into 4 Scalars at *begin to *end
0219        */
0220        template <class IT>
0221        void GetCoordinates( IT begin, IT end ) const
0222        { IT a = begin; IT b = ++begin; IT c = ++begin; IT d = ++begin;
0223        (void)end;
0224        assert (++begin==end);
0225        GetCoordinates (*a,*b,*c,*d);
0226        }
0227
0228        /**
0229           get internal data into 4 Scalars at *begin
0230        */
0231        template <class IT>
0232        void GetCoordinates( IT begin ) const {
0233           Scalar a,b,c,d = 0;
0234           GetCoordinates (a,b,c,d);
0235           *begin++ = a;
0236           *begin++ = b;
0237           *begin++ = c;
0238           *begin   = d;
0239        }
0240
0241        /**
0242           set the values of the vector from the cartesian components (x,y,z,t)
0243           (if the vector is held in another coordinates, like (Pt,eta,phi,m)
0244           then (x, y, z, t) are converted to that form)
0245        */
0246        LorentzVector<CoordSystem>& SetXYZT (Scalar xx, Scalar yy, Scalar zz, Scalar tt) {
0247           fCoordinates.SetPxPyPzE(xx,yy,zz,tt);
0248           return *this;
0249        }
0250        LorentzVector<CoordSystem>& SetPxPyPzE (Scalar xx, Scalar yy, Scalar zz, Scalar ee) {
0251           fCoordinates.SetPxPyPzE(xx,yy,zz,ee);
0252           return *this;
0253        }
0254
0255        // ------------------- Equality -----------------
0256
0257        /**
0258           Exact equality
0259        */
0260        bool operator==(const LorentzVector & rhs) const {
0261           return fCoordinates==rhs.fCoordinates;
0262        }
0263        bool operator!= (const LorentzVector & rhs) const {
0264           return !(operator==(rhs));
0265        }
0266
0267        // ------ Individual element access, in various coordinate systems ------
0268
0269        /**
0270           dimension
0271        */
0272        unsigned int Dimension() const
0273        {
0274           return fDimension;
0275        };
0276
0277        // individual coordinate accessors in various coordinate systems
0278
0279        /**
0280           spatial X component
0281        */
0282        Scalar Px() const  { return fCoordinates.Px(); }
0283        Scalar X()  const  { return fCoordinates.Px(); }
0284        /**
0285           spatial Y component
0286        */
0287        Scalar Py() const { return fCoordinates.Py(); }
0288        Scalar Y()  const { return fCoordinates.Py(); }
0289        /**
0290           spatial Z component
0291        */
0292        Scalar Pz() const { return fCoordinates.Pz(); }
0293        Scalar Z()  const { return fCoordinates.Pz(); }
0294        /**
0295           return 4-th component (time, or energy for a 4-momentum vector)
0296        */
0297        Scalar E()  const { return fCoordinates.E(); }
0298        Scalar T()  const { return fCoordinates.E(); }
0299        /**
0300           return magnitude (mass) squared  M2 = T**2 - X**2 - Y**2 - Z**2
0301           (we use -,-,-,+ metric)
0302        */
0303        Scalar M2()   const { return fCoordinates.M2(); }
0304        /**
0305           return magnitude (mass) using the  (-,-,-,+)  metric.
0306           If M2 is negative (space-like vector) a GenVector_exception
0307           is suggested and if continuing, - sqrt( -M2) is returned
0308        */
0309        Scalar M() const    { return fCoordinates.M();}
0310        /**
0311           return the spatial (3D) magnitude ( sqrt(X**2 + Y**2 + Z**2) )
0312        */
0313        Scalar R() const { return fCoordinates.R(); }
0314        Scalar P() const { return fCoordinates.R(); }
0315        /**
0316           return the square of the spatial (3D) magnitude ( X**2 + Y**2 + Z**2 )
0317        */
0318        Scalar P2() const { return P() * P(); }
0319        /**
0320           return the square of the transverse spatial component ( X**2 + Y**2 )
0321        */
0322        Scalar Perp2( ) const { return fCoordinates.Perp2();}
0323
0324        /**
0325           return the  transverse spatial component sqrt ( X**2 + Y**2 )
0326        */
0327        Scalar Pt()  const { return fCoordinates.Pt(); }
0328        Scalar Rho() const { return fCoordinates.Pt(); }
0329
0330        /**
0331           return the transverse mass squared
0332           \f[ m_t^2 = E^2 - p{_z}^2 \f]
0333        */
0334        Scalar Mt2() const { return fCoordinates.Mt2(); }
0335
0336        /**
0337           return the transverse mass
0338           \f[ \sqrt{ m_t^2 = E^2 - p{_z}^2} X sign(E^ - p{_z}^2) \f]
0339        */
0340        Scalar Mt() const { return fCoordinates.Mt(); }
0341
0342        /**
0343           return the transverse energy squared
0344           \f[ e_t = \frac{E^2 p_{\perp}^2 }{ |p|^2 } \f]
0345        */
0346        Scalar Et2() const { return fCoordinates.Et2(); }
0347
0348        /**
0349           return the transverse energy
0350           \f[ e_t = \sqrt{ \frac{E^2 p_{\perp}^2 }{ |p|^2 } } X sign(E) \f]
0351        */
0352        Scalar Et() const { return fCoordinates.Et(); }
0353
0354        /**
0355           azimuthal  Angle
0356        */
0357        Scalar Phi() const  { return fCoordinates.Phi();}
0358
0359        /**
0360           polar Angle
0361        */
0362        Scalar Theta() const { return fCoordinates.Theta(); }
0363
0364        /**
0365           pseudorapidity
0366           \f[ \eta = - \ln { \tan { \frac { \theta} {2} } } \f]
0367        */
0368        Scalar Eta() const { return fCoordinates.Eta(); }
0369
0370        /**
0371           deltaRapidity between this and vector v
0372           \f[ \Delta R = \sqrt { \Delta \eta ^2 + \Delta \phi ^2 } \f]
0373           \param useRapidity true to use Rapidity(), false to use Eta()
0374        */
0375        template <class OtherLorentzVector>
0376        Scalar DeltaR(const OtherLorentzVector &v, const bool useRapidity = false) const
0377        {
0378           const double delta = useRapidity ? Rapidity() - v.Rapidity() : Eta() - v.Eta();
0379           double dphi = Phi() - v.Phi();
0380           // convert dphi angle to the interval (-PI,PI]
0381           if (dphi > TMath::Pi() || dphi <= -TMath::Pi()) {
0382              if (dphi > 0) {
0383                 int n = static_cast<int>(dphi / TMath::TwoPi() + 0.5);
0384                 dphi -= TMath::TwoPi() * n;
0385              } else {
0386                 int n = static_cast<int>(0.5 - dphi / TMath::TwoPi());
0387                 dphi += TMath::TwoPi() * n;
0388              }
0389           }
0390           return std::sqrt(delta * delta + dphi * dphi);
0391        }
0392
0393        /**
0394           get the spatial components of the Vector in a
0395           DisplacementVector based on Cartesian Coordinates
0396        */
0397        ::ROOT::Math::DisplacementVector3D<Cartesian3D<Scalar> > Vect() const {
0398           return ::ROOT::Math::DisplacementVector3D<Cartesian3D<Scalar> >( X(), Y(), Z() );
0399        }
0400
0401        // ------ Operations combining two Lorentz vectors ------
0402
0403        /**
0404           scalar (Dot) product of two LorentzVector vectors (metric is -,-,-,+)
0405           Enable the product using any other LorentzVector implementing
0406           the x(), y() , y() and t() member functions
0407           \param  q  any LorentzVector implementing the x(), y() , z() and t()
0408           member functions
0409           \return the result of v.q of type according to the base scalar type of v
0410        */
0411
0412        template< class OtherLorentzVector >
0413        Scalar Dot(const OtherLorentzVector & q) const {
0414           return t()*q.t() - x()*q.x() - y()*q.y() - z()*q.z();
0415        }
0416
0417        /**
0418           Self addition with another Vector ( v+= q )
0419           Enable the addition with any other LorentzVector
0420           \param  q  any LorentzVector implementing the x(), y() , z() and t()
0421           member functions
0422        */
0423       template< class OtherLorentzVector >
0424       inline LorentzVector & operator += ( const OtherLorentzVector & q)
0425        {
0426           SetXYZT( x() + q.x(), y() + q.y(), z() + q.z(), t() + q.t()  );
0427           return *this;
0428        }
0429
0430        /**
0431           Self subtraction of another Vector from this ( v-= q )
0432           Enable the addition with any other LorentzVector
0433           \param  q  any LorentzVector implementing the x(), y() , z() and t()
0434           member functions
0435        */
0436        template< class OtherLorentzVector >
0437        LorentzVector & operator -= ( const OtherLorentzVector & q) {
0438           SetXYZT( x() - q.x(), y() - q.y(), z() - q.z(), t() - q.t()  );
0439           return *this;
0440        }
0441
0442        /**
0443           addition of two LorentzVectors (v3 = v1 + v2)
0444           Enable the addition with any other LorentzVector
0445           \param  v2  any LorentzVector implementing the x(), y() , z() and t()
0446           member functions
0447           \return a new LorentzVector of the same type as v1
0448        */
0449        template<class OtherLorentzVector>
0450        LorentzVector  operator +  ( const OtherLorentzVector & v2) const
0451        {
0452           LorentzVector<CoordinateType> v3(*this);
0453           v3 += v2;
0454           return v3;
0455        }
0456
0457        /**
0458           subtraction of two LorentzVectors (v3 = v1 - v2)
0459           Enable the subtraction of any other LorentzVector
0460           \param  v2  any LorentzVector implementing the x(), y() , z() and t()
0461           member functions
0462           \return a new LorentzVector of the same type as v1
0463        */
0464        template<class OtherLorentzVector>
0465        LorentzVector  operator -  ( const OtherLorentzVector & v2) const {
0466           LorentzVector<CoordinateType> v3(*this);
0467           v3 -= v2;
0468           return v3;
0469        }
0470
0471        //--- scaling operations ------
0472
0473        /**
0474           multiplication by a scalar quantity v *= a
0475        */
0476        LorentzVector & operator *= (Scalar a) {
0477           fCoordinates.Scale(a);
0478           return *this;
0479        }
0480
0481        /**
0482           division by a scalar quantity v /= a
0483        */
0484        LorentzVector & operator /= (Scalar a) {
0485           fCoordinates.Scale(1/a);
0486           return *this;
0487        }
0488
0489        /**
0490           product of a LorentzVector by a scalar quantity
0491           \param a  scalar quantity of type a
0492           \return a new mathcoreLorentzVector q = v * a same type as v
0493        */
0494        LorentzVector operator * ( const Scalar & a) const {
0495           LorentzVector tmp(*this);
0496           tmp *= a;
0497           return tmp;
0498        }
0499
0500        /**
0501           Divide a LorentzVector by a scalar quantity
0502           \param a  scalar quantity of type a
0503           \return a new mathcoreLorentzVector q = v / a same type as v
0504        */
0505        LorentzVector<CoordSystem> operator / ( const Scalar & a) const {
0506           LorentzVector<CoordSystem> tmp(*this);
0507           tmp /= a;
0508           return tmp;
0509        }
0510
0511        /**
0512           Negative of a LorentzVector (q = - v )
0513           \return a new LorentzVector with opposite direction and time
0514        */
0515        LorentzVector operator - () const {
0516           //LorentzVector<CoordinateType> v(*this);
0517           //v.Negate();
0518           return operator*( Scalar(-1) );
0519        }
0520        LorentzVector operator + () const {
0521           return *this;
0522        }
0523
0524        // ---- Relativistic Properties ----
0525
0526        /**
0527           Rapidity relative to the Z axis:  .5 log [(E+Pz)/(E-Pz)]
0528        */
0529        Scalar Rapidity() const {
0530           // TODO - It would be good to check that E > Pz and use the Throw()
0531           //        mechanism or at least load a NAN if not.
0532           //        We should then move the code to a .cpp file.
0533           const Scalar ee  = E();
0534           const Scalar ppz = Pz();
0535           using std::log;
0536           return Scalar(0.5) * log((ee + ppz) / (ee - ppz));
0537        }
0538
0539        /**
0540           Rapidity in the direction of travel: atanh (|P|/E)=.5 log[(E+P)/(E-P)]
0541        */
0542        Scalar ColinearRapidity() const {
0543           // TODO - It would be good to check that E > P and use the Throw()
0544           //        mechanism or at least load a NAN if not.
0545           const Scalar ee = E();
0546           const Scalar pp = P();
0547           using std::log;
0548           return Scalar(0.5) * log((ee + pp) / (ee - pp));
0549        }
0550
0551        /**
0552           Determine if momentum-energy can represent a physical massive particle
0553        */
0554        bool isTimelike( ) const {
0555           Scalar ee = E(); Scalar pp = P(); return ee*ee > pp*pp;
0556        }
0557
0558        /**
0559           Determine if momentum-energy can represent a massless particle
0560        */
0561        bool isLightlike( Scalar tolerance
0562                          = 100*std::numeric_limits<Scalar>::epsilon() ) const {
0563           Scalar ee = E(); Scalar pp = P(); Scalar delta = ee-pp;
0564           if ( ee==0 ) return pp==0;
0565           return delta*delta < tolerance * ee*ee;
0566        }
0567
0568        /**
0569           Determine if momentum-energy is spacelike, and represents a tachyon
0570        */
0571        bool isSpacelike( ) const {
0572           Scalar ee = E(); Scalar pp = P(); return ee*ee < pp*pp;
0573        }
0574
0575        typedef DisplacementVector3D< Cartesian3D<Scalar> > BetaVector;
0576
0577        /**
0578           The beta vector for the boost that would bring this vector into
0579           its center of mass frame (zero momentum)
0580        */
0581        BetaVector BoostToCM( ) const {
0582           if (E() == 0) {
0583              if (P() == 0) {
0584                 return BetaVector();
0585              } else {
0586                 // TODO - should attempt to Throw with msg about
0587                 // boostVector computed for LorentzVector with t=0
0588                 return -Vect()/E();
0589              }
0590           }
0591           if (M2() <= 0) {
0592              // TODO - should attempt to Throw with msg about
0593              // boostVector computed for a non-timelike LorentzVector
0594           }
0595           return -Vect()/E();
0596        }
0597
0598        /**
0599           The beta vector for the boost that would bring this vector into
0600           its center of mass frame (zero momentum)
0601        */
0602        template <class Other4Vector>
0603        BetaVector BoostToCM(const Other4Vector& v ) const {
0604           Scalar eSum = E() + v.E();
0605           DisplacementVector3D< Cartesian3D<Scalar> > vecSum = Vect() + v.Vect();
0606           if (eSum == 0) {
0607              if (vecSum.Mag2() == 0) {
0608                 return BetaVector();
0609              } else {
0610                 // TODO - should attempt to Throw with msg about
0611                 // boostToCM computed for two 4-vectors with combined t=0
0612                 return BetaVector(vecSum/eSum);
0613              }
0614              // TODO - should attempt to Throw with msg about
0615              // boostToCM computed for two 4-vectors with combined e=0
0616           }
0617           return BetaVector (vecSum * (-1./eSum));
0618        }
0619
0620        //beta and gamma
0621
0622        /**
0623            Return beta scalar value
0624        */
0625        Scalar Beta() const {
0626           if ( E() == 0 ) {
0627              if ( P2() == 0)
0628                 // to avoid Nan
0629                 return 0;
0630              else {
0631                 GenVector_Throw(
0632                    "LorentzVector::Beta() - beta computed for LorentzVector with t = 0. Return an Infinite result");
0633                 return 1./E();
0634              }
0635           }
0636           if ( M2() <= 0 ) {
0637              GenVector_Throw("LorentzVector::Beta() - beta computed for non-timelike LorentzVector . Result is "
0638                              "physically meaningless");
0639           }
0640           return P() / E();
0641        }
0642        /**
0643            Return Gamma scalar value
0644        */
0645        Scalar Gamma() const {
0646           const Scalar v2 = P2();
0647           const Scalar t2 = E() * E();
0648           if (E() == 0) {
0649              if ( P2() == 0) {
0650                 return 1;
0651              } else {
0652                 GenVector_Throw(
0653                    "LorentzVector::Gamma() - gamma computed for LorentzVector with t = 0. Return a zero result");
0654              }
0655           }
0656           if ( t2 < v2 ) {
0657              GenVector_Throw("LorentzVector::Gamma() - gamma computed for a spacelike LorentzVector. Imaginary result");
0658              return 0;
0659           }
0660           else if ( t2 == v2 ) {
0661              GenVector_Throw("LorentzVector::Gamma() - gamma computed for a lightlike LorentzVector. Infinite result");
0662           }
0663           using std::sqrt;
0664           return Scalar(1) / sqrt(Scalar(1) - v2 / t2);
0665        } /* gamma */
0666
0667
0668        // Method providing limited backward name compatibility with CLHEP ----
0669
0670        Scalar x()     const { return fCoordinates.Px();     }
0671        Scalar y()     const { return fCoordinates.Py();     }
0672        Scalar z()     const { return fCoordinates.Pz();     }
0673        Scalar t()     const { return fCoordinates.E();      }
0674        Scalar px()    const { return fCoordinates.Px();     }
0675        Scalar py()    const { return fCoordinates.Py();     }
0676        Scalar pz()    const { return fCoordinates.Pz();     }
0677        Scalar e()     const { return fCoordinates.E();      }
0678        Scalar r()     const { return fCoordinates.R();      }
0679        Scalar theta() const { return fCoordinates.Theta();  }
0680        Scalar phi()   const { return fCoordinates.Phi();    }
0681        Scalar rho()   const { return fCoordinates.Rho();    }
0682        Scalar eta()   const { return fCoordinates.Eta();    }
0683        Scalar pt()    const { return fCoordinates.Pt();     }
0684        Scalar perp2() const { return fCoordinates.Perp2();  }
0685        Scalar mag2()  const { return fCoordinates.M2();     }
0686        Scalar mag()   const { return fCoordinates.M();      }
0687        Scalar mt()    const { return fCoordinates.Mt();     }
0688        Scalar mt2()   const { return fCoordinates.Mt2();    }
0689
0690
0691        // Methods  requested by CMS ---
0692        Scalar energy() const { return fCoordinates.E();      }
0693        Scalar mass()   const { return fCoordinates.M();      }
0694        Scalar mass2()  const { return fCoordinates.M2();     }
0695
0696
0697        /**
0698           Methods setting a Single-component
0699           Work only if the component is one of which the vector is represented.
0700           For example SetE will work for a PxPyPzE Vector but not for a PxPyPzM Vector.
0701        */
0702        LorentzVector<CoordSystem>& SetE  ( Scalar a )  { fCoordinates.SetE  (a); return *this; }
0703        LorentzVector<CoordSystem>& SetEta( Scalar a )  { fCoordinates.SetEta(a); return *this; }
0704        LorentzVector<CoordSystem>& SetM  ( Scalar a )  { fCoordinates.SetM  (a); return *this; }
0705        LorentzVector<CoordSystem>& SetPhi( Scalar a )  { fCoordinates.SetPhi(a); return *this; }
0706        LorentzVector<CoordSystem>& SetPt ( Scalar a )  { fCoordinates.SetPt (a); return *this; }
0707        LorentzVector<CoordSystem>& SetPx ( Scalar a )  { fCoordinates.SetPx (a); return *this; }
0708        LorentzVector<CoordSystem>& SetPy ( Scalar a )  { fCoordinates.SetPy (a); return *this; }
0709        LorentzVector<CoordSystem>& SetPz ( Scalar a )  { fCoordinates.SetPz (a); return *this; }
0710        LorentzVector<CoordSystem> &SetT(Scalar a) { return SetE(a); }
0711        LorentzVector<CoordSystem> &SetX(Scalar a) { return SetPx(a); }
0712        LorentzVector<CoordSystem> &SetY(Scalar a) { return SetPy(a); }
0713        LorentzVector<CoordSystem> &SetZ(Scalar a) { return SetPz(a); }
0714
0715     private:
0716
0717        CoordSystem  fCoordinates;    // internal coordinate system
0718        static constexpr unsigned int fDimension = CoordinateType::Dimension;
0719
0720     };  // LorentzVector<>
0721
0722
0723
0724   // global methods
0725
0726     /**
0727        Scale of a LorentzVector with a scalar quantity a
0728        \param a  scalar quantity of type a
0729        \param v  LorentzVector based on any coordinate system
0730        \return a new mathcoreLorentzVector q = v * a same type as v
0731      */
0732     template <class CoordSystem>
0733     inline LorentzVector<CoordSystem>
0734     operator*(const typename LorentzVector<CoordSystem>::Scalar &a, const LorentzVector<CoordSystem> &v)
0735     {
0736        LorentzVector<CoordSystem> tmp(v);
0737        tmp *= a;
0738        return tmp;
0739     }
0740
0741     /**
0742        pair (p+ p-) acoplanarity `alpha = 1 - |phi+ - phi-|/pi`.
0743        \param pp p+, LorentzVector based on any coordinate system
0744        \param pm p-, LorentzVector based on any coordinate system
0745        \return a scalar
0746        \ingroup GenVector
0747        \see http://doi.org/10.1103/PhysRevLett.121.212301
0748      */
0749     template <class CoordSystem>
0750     typename LorentzVector<CoordSystem>::Scalar
0751     Acoplanarity(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0752     {
0753        auto dphi = pp.Phi() - pm.Phi();
0754        // convert dphi angle to the interval (-PI,PI]
0755        if (dphi > TMath::Pi() || dphi <= -TMath::Pi()) {
0756           if (dphi > 0) {
0757              int n = static_cast<int>(dphi / TMath::TwoPi() + 0.5);
0758              dphi -= TMath::TwoPi() * n;
0759           } else {
0760              int n = static_cast<int>(0.5 - dphi / TMath::TwoPi());
0761              dphi += TMath::TwoPi() * n;
0762           }
0763        }
0764        return 1 - std::abs(dphi) / TMath::Pi();
0765     }
0766
0767     /**
0768        pair (p+ p-) vectorial asymmetry `Av = |Pt+ - Pt-|/|Pt+ + Pt-|`.
0769        In an experimental setting, it reflects a convolution of the experimental resolutions
0770        of particle energy and azimuthal angle measurement.
0771        \param pp p+, LorentzVector based on any coordinate system
0772        \param pm p-, LorentzVector based on any coordinate system
0773        \return a scalar. Returns -1 if both momenta are exactly mirrored.
0774        \ingroup GenVector
0775        \see http://doi.org/10.1103/PhysRevLett.121.212301, https://doi.org/10.1103/PhysRevD.99.093013
0776      */
0777     template <class CoordSystem>
0778     typename LorentzVector<CoordSystem>::Scalar
0779     AsymmetryVectorial(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0780     {
0781        ROOT::Math::XYVector vp(pp.Px(), pp.Py());
0782        ROOT::Math::XYVector vm(pm.Px(), pm.Py());
0783        auto denom = (vp + vm).R();
0784        if (denom == 0.)
0785           return -1;
0786        return (vp - vm).R() / denom;
0787     }
0788
0789     /**
0790        pair (p+ p-) scalar asymmetry `As = ||Pt+| - |Pt-|/||Pt+| + |Pt-||`.
0791        Measures the relative difference in transverse momentum of the pair, e.g. two photons,
0792        and would be ideally zero for two back-to-back photons.
0793        \param pp p+, LorentzVector based on any coordinate system
0794        \param pm p-, LorentzVector based on any coordinate system
0795        \return a scalar. Returns 0 if both transverse momenta are zero
0796        \ingroup GenVector
0797        \see http://doi.org/10.1103/PhysRevLett.121.212301, https://doi.org/10.1103/PhysRevD.99.093013
0798      */
0799     template< class CoordSystem >
0800     typename LorentzVector<CoordSystem>::Scalar AsymmetryScalar(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0801     {
0802        auto ptp = pp.Pt();
0803        auto ptm = pm.Pt();
0804        if (ptp == 0 && ptm == 0)
0805           return 0;
0806        return std::abs(ptp - ptm) / (ptp + ptm);
0807     }
0808
0809     // ------------- I/O to/from streams -------------
0810
0811     template< class char_t, class traits_t, class Coords >
0812     inline
0813     std::basic_ostream<char_t,traits_t> &
0814     operator << ( std::basic_ostream<char_t,traits_t> & os
0815                   , LorentzVector<Coords> const & v
0816        )
0817     {
0818        if( !os )  return os;
0819
0820        typename Coords::Scalar a, b, c, d;
0821        v.GetCoordinates(a, b, c, d);
0822
0823        if( detail::get_manip( os, detail::bitforbit ) )  {
0824         detail::set_manip( os, detail::bitforbit, '\00' );
0825         // TODO: call MF's bitwise-accurate functions on each of a, b, c, d
0826        }
0827        else  {
0828           os << detail::get_manip( os, detail::open  ) << a
0829              << detail::get_manip( os, detail::sep   ) << b
0830              << detail::get_manip( os, detail::sep   ) << c
0831              << detail::get_manip( os, detail::sep   ) << d
0832              << detail::get_manip( os, detail::close );
0833        }
0834
0835        return os;
0836
0837     }  // op<< <>()
0838
0839
0840      template< class char_t, class traits_t, class Coords >
0841      inline
0842      std::basic_istream<char_t,traits_t> &
0843      operator >> ( std::basic_istream<char_t,traits_t> & is
0844                    , LorentzVector<Coords> & v
0845         )
0846      {
0847         if( !is )  return is;
0848
0849         typename Coords::Scalar a, b, c, d;
0850
0851         if( detail::get_manip( is, detail::bitforbit ) )  {
0852            detail::set_manip( is, detail::bitforbit, '\00' );
0853            // TODO: call MF's bitwise-accurate functions on each of a, b, c
0854         }
0855         else  {
0856            detail::require_delim( is, detail::open  );  is >> a;
0857            detail::require_delim( is, detail::sep   );  is >> b;
0858            detail::require_delim( is, detail::sep   );  is >> c;
0859            detail::require_delim( is, detail::sep   );  is >> d;
0860            detail::require_delim( is, detail::close );
0861         }
0862
0863         if( is )
0864            v.SetCoordinates(a, b, c, d);
0865         return is;
0866
0867      }  // op>> <>()
0868
0869     // Structured bindings
0870     template <std::size_t I, class CoordSystem>
0871     typename CoordSystem::Scalar get(LorentzVector<CoordSystem> const& p)
0872     {
0873        static_assert(I < 4);
0874        if constexpr (I == 0) {
0875           return p.X();
0876        } else if constexpr (I == 1) {
0877           return p.Y();
0878        } else if constexpr (I == 2) {
0879           return p.Z();
0880        } else {
0881           return p.E();
0882        }
0883     }
0884
0885   } // end namespace Math
0886
0887 } // end namespace ROOT
0888
0889 // Structured bindings
0890 #include <tuple>
0891 namespace std {
0892    template <class CoordSystem>
0893    struct tuple_size<ROOT::Math::LorentzVector<CoordSystem>> : integral_constant<size_t, 4> {};
0894    template <size_t I, class CoordSystem>
0895    struct tuple_element<I, ROOT::Math::LorentzVector<CoordSystem>> {
0896       static_assert(I < 4);
0897       using type = typename CoordSystem::Scalar;
0898    };
0899 }
0900
0901 #include <sstream>
0902 namespace cling
0903 {
0904 template<typename CoordSystem>
0905 std::string printValue(const ROOT::Math::LorentzVector<CoordSystem> *v)
0906 {
0907    std::stringstream s;
0908    s << *v;
0909    return s.str();
0910 }
0911
0912 } // end namespace cling
0913
0914 #endif
0915
0916 //#include "Math/GenVector/LorentzVectorOperations.h"