Math/GenVector/VectorUtil.h

0001 // @(#)root/mathcore:$Id: 9ef2a4a7bd1b62c1293920c2af2f64791c75bdd8 $
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003
0004
0005 /**********************************************************************
0006  *                                                                    *
0007  * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0008  *                                                                    *
0009  *                                                                    *
0010  **********************************************************************/
0011
0012 // Header file for Vector Utility functions
0013 //
0014 // Created by: moneta  at Tue May 31 21:10:29 2005
0015 //
0016 // Last update: Tue May 31 21:10:29 2005
0017 //
0018 #ifndef ROOT_Math_GenVector_VectorUtil
0019 #define ROOT_Math_GenVector_VectorUtil  1
0020
0021 #include "TMath.h"
0022
0023 #include "Math/GenVector/Boost.h"
0024
0025 namespace ROOT {
0026
0027    namespace Math {
0028
0029
0030       // utility functions for vector classes
0031
0032
0033
0034       /**
0035        Global Helper functions for generic Vector classes. Any Vector classes implementing some defined member functions,
0036        like  Phi() or Eta() or mag() can use these functions.
0037        The functions returning a scalar value, returns always double precision number even if the vector are
0038        based on another precision type
0039
0040        @ingroup GenVector
0041
0042        @see GenVector
0043        */
0044
0045
0046       namespace VectorUtil {
0047
0048       /// \addtogroup GenVector
0049       /// @{
0050
0051          // methods for 3D vectors
0052
0053          /**
0054           Find aximutal Angle difference between two generic vectors ( v2.Phi() - v1.Phi() )
0055           The only requirements on the Vector classes is that they implement the Phi() method
0056           \param v1  Vector of any type implementing the Phi() operator
0057           \param v2  Vector of any type implementing the Phi() operator
0058           \return  Phi difference
0059           \f[ \Delta \phi = \phi_2 - \phi_1 \f]
0060           */
0061          template <class Vector1, class Vector2>
0062          inline typename Vector1::Scalar DeltaPhi( const Vector1 & v1, const Vector2 & v2) {
0063             typename Vector1::Scalar dphi = v2.Phi() - v1.Phi();
0064             if (dphi > TMath::Pi()) {
0065                dphi -= 2.0 * TMath::Pi();
0066             } else if (dphi <= -TMath::Pi()) {
0067                dphi += 2.0 * TMath::Pi();
0068             }
0069             return dphi;
0070          }
0071
0072
0073
0074          /**
0075           Find square of the difference in pseudorapidity (Eta) and Phi between two generic vectors
0076           The only requirements on the Vector classes is that they implement the Phi() and Eta() method
0077           \param v1  Vector 1
0078           \param v2  Vector 2
0079           \return   Angle between the two vectors
0080           \f[ \Delta R2 = ( \Delta \phi )^2 + ( \Delta \eta )^2  \f]
0081           */
0082          template <class Vector1, class Vector2>
0083          inline typename Vector1::Scalar DeltaR2( const Vector1 & v1, const Vector2 & v2) {
0084             typename Vector1::Scalar dphi = DeltaPhi(v1,v2);
0085             typename Vector1::Scalar deta = v2.Eta() - v1.Eta();
0086             return dphi*dphi + deta*deta;
0087          }
0088
0089      /**
0090       Find square of the difference in true rapidity (y) and Phi between two generic vectors
0091       The only requirements on the Vector classes is that they implement the Phi() and Rapidity() method
0092       \param v1  Vector 1
0093       \param v2  Vector 2
0094       \return   Angle between the two vectors
0095       \f[ \Delta R2 = ( \Delta \phi )^2 + ( \Delta \y )^2  \f]
0096       */
0097      template <class Vector1, class Vector2>
0098      inline typename Vector1::Scalar DeltaR2RapidityPhi( const Vector1 & v1, const Vector2 & v2) {
0099         typename Vector1::Scalar dphi = DeltaPhi(v1,v2);
0100         typename Vector1::Scalar drap = v2.Rapidity() - v1.Rapidity();
0101         return dphi*dphi + drap*drap;
0102      }
0103
0104          /**
0105           Find difference in pseudorapidity (Eta) and Phi between two generic vectors
0106           The only requirements on the Vector classes is that they implement the Phi() and Eta() method
0107           \param v1  Vector 1
0108           \param v2  Vector 2
0109           \return   Angle between the two vectors
0110           \f[ \Delta R = \sqrt{  ( \Delta \phi )^2 + ( \Delta \eta )^2 } \f]
0111           */
0112          template <class Vector1, class Vector2>
0113          inline typename Vector1::Scalar DeltaR( const Vector1 & v1, const Vector2 & v2) {
0114             using std::sqrt;
0115             return sqrt( DeltaR2(v1,v2) );
0116          }
0117
0118     /**
0119           Find difference in Rapidity (y) and Phi between two generic vectors
0120           The only requirements on the Vector classes is that they implement the Phi() and Rapidity() method
0121           \param v1  Vector 1
0122           \param v2  Vector 2
0123           \return   Angle between the two vectors
0124           \f[ \Delta R = \sqrt{  ( \Delta \phi )^2 + ( \Delta y )^2 } \f],
0125           */
0126          template <class Vector1, class Vector2>
0127          inline typename Vector1::Scalar DeltaRapidityPhi( const Vector1 & v1, const Vector2 & v2) {
0128             using std::sqrt;
0129             return sqrt( DeltaR2RapidityPhi(v1,v2) );
0130          }
0131
0132          /**
0133           Find CosTheta Angle between two generic 3D vectors
0134           pre-requisite: vectors implement the X(), Y() and Z()
0135           \param v1  Vector v1
0136           \param v2  Vector v2
0137           \return   cosine of Angle between the two vectors
0138           \f[ \cos \theta = \frac { \vec{v1} \cdot \vec{v2} }{ | \vec{v1} | | \vec{v2} | } \f]
0139           */
0140          // this cannot be made all generic since Mag2() for 2, 3 or 4 D is different
0141          // need to have a specialization for polar Coordinates ??
0142          template <class Vector1, class Vector2>
0143          double CosTheta( const Vector1 &  v1, const Vector2  & v2) {
0144             double arg;
0145             double v1_r2 = v1.X()*v1.X() + v1.Y()*v1.Y() + v1.Z()*v1.Z();
0146             double v2_r2 = v2.X()*v2.X() + v2.Y()*v2.Y() + v2.Z()*v2.Z();
0147             double ptot2 = v1_r2*v2_r2;
0148             if(ptot2 <= 0) {
0149                arg = 0.0;
0150             }else{
0151                double pdot = v1.X()*v2.X() + v1.Y()*v2.Y() + v1.Z()*v2.Z();
0152                using std::sqrt;
0153                arg = pdot/sqrt(ptot2);
0154                if(arg >  1.0) arg =  1.0;
0155                if(arg < -1.0) arg = -1.0;
0156             }
0157             return arg;
0158          }
0159
0160
0161          /**
0162           Find Angle between two vectors.
0163           Use the CosTheta() function
0164           \param v1  Vector v1
0165           \param v2  Vector v2
0166           \return   Angle between the two vectors
0167           \f[ \theta = \cos ^{-1} \frac { \vec{v1} \cdot \vec{v2} }{ | \vec{v1} | | \vec{v2} | } \f]
0168           */
0169          template <class Vector1, class Vector2>
0170          inline double Angle( const  Vector1 & v1, const Vector2 & v2) {
0171             using std::acos;
0172             return acos( CosTheta(v1, v2) );
0173          }
0174
0175          /**
0176           Find the projection of v along the given direction u.
0177           \param v  Vector v for which the propjection is to be found
0178           \param u  Vector specifying the direction
0179           \return   Vector projection (same type of v)
0180           \f[ \vec{proj} = \frac{ \vec{v}  \cdot \vec{u} }{|\vec{u}|}\vec{u} \f]
0181           Precondition is that Vector1 implements Dot function and Vector2 implements X(),Y() and Z()
0182           */
0183          template <class Vector1, class Vector2>
0184          Vector1 ProjVector( const  Vector1 & v, const Vector2 & u) {
0185             double magU2 = u.X()*u.X() + u.Y()*u.Y() + u.Z()*u.Z();
0186             if (magU2 == 0) return Vector1(0,0,0);
0187             double d = v.Dot(u)/magU2;
0188             return Vector1( u.X() * d, u.Y() * d, u.Z() * d);
0189          }
0190
0191          /**
0192           Find the vector component of v perpendicular to the given direction of u
0193           \param v  Vector v for which the perpendicular component is to be found
0194           \param u  Vector specifying the direction
0195           \return   Vector component of v which is perpendicular to u
0196           \f[ \vec{perp} = \vec{v} -  \frac{ \vec{v}  \cdot \vec{u} }{|\vec{u}|}\vec{u} \f]
0197           Precondition is that Vector1 implements Dot function and Vector2 implements X(),Y() and Z()
0198           */
0199          template <class Vector1, class Vector2>
0200          inline Vector1 PerpVector( const  Vector1 & v, const Vector2 & u) {
0201             return v - ProjVector(v,u);
0202          }
0203
0204          /**
0205           Find the magnitude square of the vector component of v perpendicular to the given direction of u
0206           \param v  Vector v for which the perpendicular component is to be found
0207           \param u  Vector specifying the direction
0208           \return   square value of the component of v which is perpendicular to u
0209           \f[ perp = | \vec{v} -  \frac{ \vec{v}  \cdot \vec{u} }{|\vec{u}|}\vec{u} |^2 \f]
0210           Precondition is that Vector1 implements Dot function and Vector2 implements X(),Y() and Z()
0211           */
0212          template <class Vector1, class Vector2>
0213          inline double Perp2( const  Vector1 & v, const Vector2 & u) {
0214             double magU2 = u.X()*u.X() + u.Y()*u.Y() + u.Z()*u.Z();
0215             double prjvu = v.Dot(u);
0216             double magV2 = v.Dot(v);
0217             return magU2 > 0.0 ? magV2-prjvu*prjvu/magU2 : magV2;
0218          }
0219
0220          /**
0221           Find the magnitude of the vector component of v perpendicular to the given direction of u
0222           \param v  Vector v for which the perpendicular component is to be found
0223           \param u  Vector specifying the direction
0224           \return   value of the component of v which is perpendicular to u
0225           \f[ perp = | \vec{v} -  \frac{ \vec{v}  \cdot \vec{u} }{|\vec{u}|}\vec{u} | \f]
0226           Precondition is that Vector1 implements Dot function and Vector2 implements X(),Y() and Z()
0227           */
0228          template <class Vector1, class Vector2>
0229          inline double Perp( const  Vector1 & v, const Vector2 & u) {
0230             using std::sqrt;
0231             return sqrt(Perp2(v,u) );
0232          }
0233
0234
0235
0236          // Lorentz Vector functions
0237
0238
0239          /**
0240           return the invariant mass of two LorentzVector
0241           The only requirement on the LorentzVector is that they need to implement the
0242           X() , Y(), Z() and E() methods.
0243           \param v1 LorenzVector 1
0244           \param v2 LorenzVector 2
0245           \return invariant mass M
0246           \f[ M_{12} = \sqrt{ (\vec{v1} + \vec{v2} ) \cdot (\vec{v1} + \vec{v2} ) } \f]
0247           */
0248          template <class Vector1, class Vector2>
0249          inline typename Vector1::Scalar InvariantMass( const Vector1 & v1, const Vector2 & v2) {
0250             typedef typename  Vector1::Scalar Scalar;
0251             Scalar ee = (v1.E() + v2.E() );
0252             Scalar xx = (v1.X() + v2.X() );
0253             Scalar yy = (v1.Y() + v2.Y() );
0254             Scalar zz = (v1.Z() + v2.Z() );
0255             Scalar mm2 = ee*ee - xx*xx - yy*yy - zz*zz;
0256             using std::sqrt;
0257             return mm2 < 0.0 ? -sqrt(-mm2) : sqrt(mm2);
0258             //  PxPyPzE4D<double> q(xx,yy,zz,ee);
0259             //  return q.M();
0260             //return ( v1 + v2).mag();
0261          }
0262
0263          /// @brief Returns the square of what InvariantMass(const Vector1&, const Vector2&) would return.
0264          /// This is mostly useful for speed optimisations, because the expensive sqrt operation is skipped.
0265          template <class Vector1, class Vector2>
0266          inline typename Vector1::Scalar InvariantMass2( const Vector1 & v1, const Vector2 & v2) {
0267             typedef typename  Vector1::Scalar Scalar;
0268             Scalar ee = (v1.E() + v2.E() );
0269             Scalar xx = (v1.X() + v2.X() );
0270             Scalar yy = (v1.Y() + v2.Y() );
0271             Scalar zz = (v1.Z() + v2.Z() );
0272             Scalar mm2 = ee*ee - xx*xx - yy*yy - zz*zz;
0273             return mm2 ; // < 0.0 ? -std::sqrt(-mm2) : std::sqrt(mm2);
0274                          //  PxPyPzE4D<double> q(xx,yy,zz,ee);
0275                          //  return q.M();
0276                          //return ( v1 + v2).mag();
0277          }
0278
0279          // rotation and transformations
0280
0281          /**
0282           rotation along X axis for a generic vector by an Angle alpha
0283           returning a new vector.
0284           The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
0285           and SetXYZ methods.
0286           */
0287          template <class Vector>
0288          Vector RotateX(const Vector & v, double alpha) {
0289             if (std::fmod(alpha, 2 * TMath::Pi()) == 0.)
0290                return v;
0291             using std::sin;
0292             double sina = sin(alpha);
0293             using std::cos;
0294             double cosa = cos(alpha);
0295             double y2 = v.Y() * cosa - v.Z()*sina;
0296             double z2 = v.Z() * cosa + v.Y() * sina;
0297             Vector vrot;
0298             vrot.SetXYZ(v.X(), y2, z2);
0299             return vrot;
0300          }
0301
0302          /**
0303           rotation along Y axis for a generic vector by an Angle alpha
0304           returning a new vector.
0305           The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
0306           and SetXYZ methods.
0307           */
0308          template <class Vector>
0309          Vector RotateY(const Vector & v, double alpha) {
0310             if (std::fmod(alpha, 2 * TMath::Pi()) == 0.)
0311                return v;
0312             using std::sin;
0313             double sina = sin(alpha);
0314             using std::cos;
0315             double cosa = cos(alpha);
0316             double x2 = v.X() * cosa + v.Z() * sina;
0317             double z2 = v.Z() * cosa - v.X() * sina;
0318             Vector vrot;
0319             vrot.SetXYZ(x2, v.Y(), z2);
0320             return vrot;
0321          }
0322
0323          /**
0324           rotation along Z axis for a generic vector by an Angle alpha
0325           returning a new vector.
0326           The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
0327           and SetXYZ methods.
0328           */
0329          template <class Vector>
0330          Vector RotateZ(const Vector & v, double alpha) {
0331             if (std::fmod(alpha, 2 * TMath::Pi()) == 0.)
0332                return v;
0333             using std::sin;
0334             double sina = sin(alpha);
0335             using std::cos;
0336             double cosa = cos(alpha);
0337             double x2 = v.X() * cosa - v.Y() * sina;
0338             double y2 = v.Y() * cosa + v.X() * sina;
0339             Vector vrot;
0340             vrot.SetXYZ(x2, y2, v.Z());
0341             return vrot;
0342          }
0343
0344          /**
0345           rotation along a custom axis for a generic vector by an Angle alpha (in rad)
0346           returning a new vector.
0347           The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
0348           and SetXYZ methods.
0349           */
0350          template <class Vector>
0351          Vector Rotate(const Vector &v, double alpha, const Vector &axis)
0352          {
0353             if (std::fmod(alpha, 2 * TMath::Pi()) == 0.)
0354                return v;
0355             const double ll = std::sqrt(axis.X() * axis.X() + axis.Y() * axis.Y() + axis.Z() * axis.Z());
0356             if (ll == 0.)
0357                GenVector_Throw("Axis Vector has zero magnitude");
0358             const double sa = std::sin(alpha);
0359             const double ca = std::cos(alpha);
0360             const double dx = axis.X() / ll;
0361             const double dy = axis.Y() / ll;
0362             const double dz = axis.Z() / ll;
0363             // clang-format off
0364             const double rot00 = (1 - ca) * dx * dx + ca     , rot01 = (1 - ca) * dx * dy - sa * dz, rot02 = (1 - ca) * dx * dz + sa * dy,
0365                          rot10 = (1 - ca) * dy * dx + sa * dz, rot11 = (1 - ca) * dy * dy + ca     , rot12 = (1 - ca) * dy * dz - sa * dx,
0366                          rot20 = (1 - ca) * dz * dx - sa * dy, rot21 = (1 - ca) * dz * dy + sa * dx, rot22 = (1 - ca) * dz * dz + ca     ;
0367             // clang-format on
0368             const double xX = v.X();
0369             const double yY = v.Y();
0370             const double zZ = v.Z();
0371             const double x2 = rot00 * xX + rot01 * yY + rot02 * zZ;
0372             const double y2 = rot10 * xX + rot11 * yY + rot12 * zZ;
0373             const double z2 = rot20 * xX + rot21 * yY + rot22 * zZ;
0374             Vector vrot;
0375             vrot.SetXYZ(x2, y2, z2);
0376             return vrot;
0377          }
0378
0379          /**
0380           rotation on a generic vector using a generic rotation class.
0381           The only requirement on the vector is that implements the
0382           X(), Y(), Z() and SetXYZ methods.
0383           The requirement on the rotation matrix is that need to implement the
0384           (i,j) operator returning the matrix element with R(0,0) = xx element
0385           */
0386          template<class Vector, class RotationMatrix>
0387          Vector Rotate(const Vector &v, const RotationMatrix & rot) {
0388             double xX = v.X();
0389             double yY = v.Y();
0390             double zZ = v.Z();
0391             double x2 =  rot(0,0)*xX + rot(0,1)*yY + rot(0,2)*zZ;
0392             double y2 =  rot(1,0)*xX + rot(1,1)*yY + rot(1,2)*zZ;
0393             double z2 =  rot(2,0)*xX + rot(2,1)*yY + rot(2,2)*zZ;
0394             Vector vrot;
0395             vrot.SetXYZ(x2,y2,z2);
0396             return vrot;
0397          }
0398
0399          /**
0400           Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost
0401           The only requirement on the vector is that implements the
0402           X(), Y(), Z(), T() and SetXYZT methods.
0403           The requirement on the boost vector is that needs to implement the
0404           X(), Y() , Z()  retorning the vector elements describing the boost
0405           The beta of the boost must be <= 1 or a nul Lorentz Vector will be returned
0406           */
0407          template <class LVector, class BoostVector>
0408          LVector boost(const LVector & v, const BoostVector & b) {
0409             double bx = b.X();
0410             double by = b.Y();
0411             double bz = b.Z();
0412             double b2 = bx*bx + by*by + bz*bz;
0413             if (b2 >= 1) {
0414                GenVector_Throw("Beta Vector supplied to set Boost represents speed >= c");
0415                return LVector();
0416             }
0417             using std::sqrt;
0418             double gamma = 1.0 / sqrt(1.0 - b2);
0419             double bp = bx*v.X() + by*v.Y() + bz*v.Z();
0420             double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
0421             double x2 = v.X() + gamma2*bp*bx + gamma*bx*v.T();
0422             double y2 = v.Y() + gamma2*bp*by + gamma*by*v.T();
0423             double z2 = v.Z() + gamma2*bp*bz + gamma*bz*v.T();
0424             double t2 = gamma*(v.T() + bp);
0425             LVector lv;
0426             lv.SetXYZT(x2,y2,z2,t2);
0427             return lv;
0428          }
0429
0430
0431          /**
0432           Boost a generic Lorentz Vector class along the X direction with a factor beta
0433           The only requirement on the vector is that implements the
0434           X(), Y(), Z(), T()  and SetXYZT methods.
0435           The beta of the boost must be <= 1 or a nul Lorentz Vector will be returned
0436           */
0437          template <class LVector, class T>
0438          LVector boostX(const LVector & v, T beta) {
0439             if (beta >= 1) {
0440                GenVector_Throw("Beta Vector supplied to set Boost represents speed >= c");
0441                return LVector();
0442             }
0443             using std::sqrt;
0444             T gamma = 1.0/ sqrt(1.0 - beta*beta);
0445             typename LVector::Scalar x2 = gamma * v.X() + gamma * beta * v.T();
0446             typename LVector::Scalar t2 = gamma * beta * v.X() + gamma * v.T();
0447
0448             LVector lv;
0449             lv.SetXYZT(x2,v.Y(),v.Z(),t2);
0450             return lv;
0451          }
0452
0453          /**
0454           Boost a generic Lorentz Vector class along the Y direction with a factor beta
0455           The only requirement on the vector is that implements the
0456           X(), Y(), Z(), T()  methods and be constructed from x,y,z,t values
0457           The beta of the boost must be <= 1 or a nul Lorentz Vector will be returned
0458           */
0459          template <class LVector>
0460          LVector boostY(const LVector & v, double beta) {
0461             if (beta >= 1) {
0462                GenVector_Throw("Beta Vector supplied to set Boost represents speed >= c");
0463                return LVector();
0464             }
0465             using std::sqrt;
0466             double gamma = 1.0/ sqrt(1.0 - beta*beta);
0467             double y2 = gamma * v.Y() + gamma * beta * v.T();
0468             double t2 = gamma * beta * v.Y() + gamma * v.T();
0469             LVector lv;
0470             lv.SetXYZT(v.X(),y2,v.Z(),t2);
0471             return lv;
0472          }
0473
0474          /**
0475           Boost a generic Lorentz Vector class along the Z direction with a factor beta
0476           The only requirement on the vector is that implements the
0477           X(), Y(), Z(), T()  methods and be constructed from x,y,z,t values
0478           The beta of the boost must be <= 1 or a nul Lorentz Vector will be returned
0479           */
0480          template <class LVector>
0481          LVector boostZ(const LVector & v, double beta) {
0482             if (beta >= 1) {
0483                GenVector_Throw("Beta Vector supplied to set Boost represents speed >= c");
0484                return LVector();
0485             }
0486             using std::sqrt;
0487             double gamma = 1.0/ sqrt(1.0 - beta*beta);
0488             double z2 = gamma * v.Z() + gamma * beta * v.T();
0489             double t2 = gamma * beta * v.Z() + gamma * v.T();
0490             LVector lv;
0491             lv.SetXYZT(v.X(),v.Y(),z2,t2);
0492             return lv;
0493          }
0494
0495          // MATRIX VECTOR MULTIPLICATION
0496          // cannot define an operator * otherwise conflicts with rotations
0497          // operations like Rotation3D * vector use Mult
0498
0499          /**
0500           Multiplications of a generic matrices with a  DisplacementVector3D of any coordinate system.
0501           Assume that the matrix implements the operator( i,j) and that it has at least         3 columns and 3 rows. There is no check on the matrix size !!
0502           */
0503          template<class Matrix, class CoordSystem, class U>
0504          inline
0505          DisplacementVector3D<CoordSystem,U> Mult (const Matrix & m, const DisplacementVector3D<CoordSystem,U> & v) {
0506             DisplacementVector3D<CoordSystem,U> vret;
0507             vret.SetXYZ( m(0,0) * v.x() + m(0,1) * v.y() + m(0,2) * v.z() ,
0508                         m(1,0) * v.x() + m(1,1) * v.y() + m(1,2) * v.z() ,
0509                         m(2,0) * v.x() + m(2,1) * v.y() + m(2,2) * v.z() );
0510             return vret;
0511          }
0512
0513
0514          /**
0515           Multiplications of a generic matrices with a generic PositionVector
0516           Assume that the matrix implements the operator( i,j) and that it has at least         3 columns and 3 rows. There is no check on the matrix size !!
0517           */
0518          template<class Matrix, class CoordSystem, class U>
0519          inline
0520          PositionVector3D<CoordSystem,U> Mult (const Matrix & m, const PositionVector3D<CoordSystem,U> & p) {
0521             DisplacementVector3D<CoordSystem,U> pret;
0522             pret.SetXYZ( m(0,0) * p.x() + m(0,1) * p.y() + m(0,2) * p.z() ,
0523                         m(1,0) * p.x() + m(1,1) * p.y() + m(1,2) * p.z() ,
0524                         m(2,0) * p.x() + m(2,1) * p.y() + m(2,2) * p.z() );
0525             return pret;
0526          }
0527
0528
0529          /**
0530           Multiplications of a generic matrices with a  LorentzVector described
0531           in any coordinate system.
0532           Assume that the matrix implements the operator( i,j) and that it has at least         4 columns and 4 rows. There is no check on the matrix size !!
0533           */
0534          // this will not be ambiguous with operator*(Scalar, LorentzVector) since that one     // Scalar is passed by value
0535          template<class CoordSystem, class Matrix>
0536          inline
0537          LorentzVector<CoordSystem> Mult (const Matrix & m, const LorentzVector<CoordSystem> & v) {
0538             LorentzVector<CoordSystem> vret;
0539             vret.SetXYZT( m(0,0)*v.x() + m(0,1)*v.y() + m(0,2)*v.z() + m(0,3)* v.t() ,
0540                          m(1,0)*v.x() + m(1,1)*v.y() + m(1,2)*v.z() + m(1,3)* v.t() ,
0541                          m(2,0)*v.x() + m(2,1)*v.y() + m(2,2)*v.z() + m(2,3)* v.t() ,
0542                          m(3,0)*v.x() + m(3,1)*v.y() + m(3,2)*v.z() + m(3,3)* v.t() );
0543             return vret;
0544          }
0545
0546
0547
0548          // non-template utility functions for all objects
0549
0550
0551          /**
0552           Return a phi angle in the interval (0,2*PI]
0553           */
0554          double Phi_0_2pi(double phi);
0555          /**
0556           Returns phi angle in the interval (-PI,PI]
0557           */
0558          double  Phi_mpi_pi(double phi);
0559
0560          // Close of doxygen group:
0561          /// @}
0562
0563       }  // end namespace Vector Util
0564
0565
0566
0567    } // end namespace Math
0568
0569 } // end namespace ROOT
0570
0571
0572 #endif /* ROOT_Math_GenVector_VectorUtil  */