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