Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-01-01 09:59:21 

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 
0711     private:
0712 
0713        CoordSystem  fCoordinates;    // internal coordinate system
0714        static constexpr unsigned int fDimension = CoordinateType::Dimension;
0715 
0716     };  // LorentzVector<>
0717 
0718 
0719 
0720   // global methods
0721 
0722     /**
0723        Scale of a LorentzVector with a scalar quantity a
0724        \param a  scalar quantity of type a
0725        \param v  LorentzVector based on any coordinate system
0726        \return a new mathcoreLorentzVector q = v * a same type as v
0727      */
0728     template <class CoordSystem>
0729     inline LorentzVector<CoordSystem>
0730     operator*(const typename LorentzVector<CoordSystem>::Scalar &a, const LorentzVector<CoordSystem> &v)
0731     {
0732        LorentzVector<CoordSystem> tmp(v);
0733        tmp *= a;
0734        return tmp;
0735     }
0736 
0737     /**
0738        pair (p+ p-) acoplanarity `alpha = 1 - |phi+ - phi-|/pi`.
0739        \param pp p+, LorentzVector based on any coordinate system
0740        \param pm p-, LorentzVector based on any coordinate system
0741        \return a scalar
0742        \ingroup GenVector
0743        \see http://doi.org/10.1103/PhysRevLett.121.212301
0744      */
0745     template <class CoordSystem>
0746     typename LorentzVector<CoordSystem>::Scalar
0747     Acoplanarity(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0748     {
0749        auto dphi = pp.Phi() - pm.Phi();
0750        // convert dphi angle to the interval (-PI,PI]
0751        if (dphi > TMath::Pi() || dphi <= -TMath::Pi()) {
0752           if (dphi > 0) {
0753              int n = static_cast<int>(dphi / TMath::TwoPi() + 0.5);
0754              dphi -= TMath::TwoPi() * n;
0755           } else {
0756              int n = static_cast<int>(0.5 - dphi / TMath::TwoPi());
0757              dphi += TMath::TwoPi() * n;
0758           }
0759        }
0760        return 1 - std::abs(dphi) / TMath::Pi();
0761     }
0762 
0763     /**
0764        pair (p+ p-) vectorial asymmetry `Av = |Pt+ - Pt-|/|Pt+ + Pt-|`.
0765        In an experimental setting, it reflects a convolution of the experimental resolutions
0766        of particle energy and azimuthal angle measurement.
0767        \param pp p+, LorentzVector based on any coordinate system
0768        \param pm p-, LorentzVector based on any coordinate system
0769        \return a scalar. Returns -1 if both momenta are exactly mirrored.
0770        \ingroup GenVector
0771        \see http://doi.org/10.1103/PhysRevLett.121.212301, https://doi.org/10.1103/PhysRevD.99.093013
0772      */
0773     template <class CoordSystem>
0774     typename LorentzVector<CoordSystem>::Scalar
0775     AsymmetryVectorial(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0776     {
0777        ROOT::Math::XYVector vp(pp.Px(), pp.Py());
0778        ROOT::Math::XYVector vm(pm.Px(), pm.Py());
0779        auto denom = (vp + vm).R();
0780        if (denom == 0.)
0781           return -1;
0782        return (vp - vm).R() / denom;
0783     }
0784 
0785     /**
0786        pair (p+ p-) scalar asymmetry `As = ||Pt+| - |Pt-|/||Pt+| + |Pt-||`.
0787        Measures the relative difference in transverse momentum of the pair, e.g. two photons,
0788        and would be ideally zero for two back-to-back photons.
0789        \param pp p+, LorentzVector based on any coordinate system
0790        \param pm p-, LorentzVector based on any coordinate system
0791        \return a scalar. Returns 0 if both transverse momenta are zero
0792        \ingroup GenVector
0793        \see http://doi.org/10.1103/PhysRevLett.121.212301, https://doi.org/10.1103/PhysRevD.99.093013
0794      */
0795     template< class CoordSystem >
0796     typename LorentzVector<CoordSystem>::Scalar AsymmetryScalar(LorentzVector<CoordSystem> const &pp, LorentzVector<CoordSystem> const &pm)
0797     {
0798        auto ptp = pp.Pt();
0799        auto ptm = pm.Pt();
0800        if (ptp == 0 && ptm == 0)
0801           return 0;
0802        return std::abs(ptp - ptm) / (ptp + ptm);
0803     }
0804 
0805     // ------------- I/O to/from streams -------------
0806 
0807     template< class char_t, class traits_t, class Coords >
0808     inline
0809     std::basic_ostream<char_t,traits_t> &
0810     operator << ( std::basic_ostream<char_t,traits_t> & os
0811                   , LorentzVector<Coords> const & v
0812        )
0813     {
0814        if( !os )  return os;
0815 
0816        typename Coords::Scalar a, b, c, d;
0817        v.GetCoordinates(a, b, c, d);
0818 
0819        if( detail::get_manip( os, detail::bitforbit ) )  {
0820         detail::set_manip( os, detail::bitforbit, '\00' );
0821         // TODO: call MF's bitwise-accurate functions on each of a, b, c, d
0822        }
0823        else  {
0824           os << detail::get_manip( os, detail::open  ) << a
0825              << detail::get_manip( os, detail::sep   ) << b
0826              << detail::get_manip( os, detail::sep   ) << c
0827              << detail::get_manip( os, detail::sep   ) << d
0828              << detail::get_manip( os, detail::close );
0829        }
0830 
0831        return os;
0832 
0833     }  // op<< <>()
0834 
0835 
0836      template< class char_t, class traits_t, class Coords >
0837      inline
0838      std::basic_istream<char_t,traits_t> &
0839      operator >> ( std::basic_istream<char_t,traits_t> & is
0840                    , LorentzVector<Coords> & v
0841         )
0842      {
0843         if( !is )  return is;
0844 
0845         typename Coords::Scalar a, b, c, d;
0846 
0847         if( detail::get_manip( is, detail::bitforbit ) )  {
0848            detail::set_manip( is, detail::bitforbit, '\00' );
0849            // TODO: call MF's bitwise-accurate functions on each of a, b, c
0850         }
0851         else  {
0852            detail::require_delim( is, detail::open  );  is >> a;
0853            detail::require_delim( is, detail::sep   );  is >> b;
0854            detail::require_delim( is, detail::sep   );  is >> c;
0855            detail::require_delim( is, detail::sep   );  is >> d;
0856            detail::require_delim( is, detail::close );
0857         }
0858 
0859         if( is )
0860            v.SetCoordinates(a, b, c, d);
0861         return is;
0862 
0863      }  // op>> <>()
0864 
0865     // Structured bindings
0866     template <std::size_t I, class CoordSystem>
0867     typename CoordSystem::Scalar get(LorentzVector<CoordSystem> const& p)
0868     {
0869        static_assert(I < 4);
0870        if constexpr (I == 0) {
0871           return p.X();
0872        } else if constexpr (I == 1) {
0873           return p.Y();
0874        } else if constexpr (I == 2) {
0875           return p.Z();
0876        } else {
0877           return p.E();
0878        }
0879     }
0880 
0881   } // end namespace Math
0882 
0883 } // end namespace ROOT
0884 
0885 // Structured bindings
0886 #include <tuple>
0887 namespace std {
0888    template <class CoordSystem>
0889    struct tuple_size<ROOT::Math::LorentzVector<CoordSystem>> : integral_constant<size_t, 4> {};
0890    template <size_t I, class CoordSystem>
0891    struct tuple_element<I, ROOT::Math::LorentzVector<CoordSystem>> {
0892       static_assert(I < 4);
0893       using type = typename CoordSystem::Scalar;
0894    };
0895 }
0896 
0897 #include <sstream>
0898 namespace cling
0899 {
0900 template<typename CoordSystem>
0901 std::string printValue(const ROOT::Math::LorentzVector<CoordSystem> *v)
0902 {
0903    std::stringstream s;
0904    s << *v;
0905    return s.str();
0906 }
0907 
0908 } // end namespace cling
0909 
0910 #endif
0911 
0912 //#include "Math/GenVector/LorentzVectorOperations.h"
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!