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