Math/GenVector/EulerAngles.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 EulerAngles
0012 //
0013 // Created by: Lorenzo Moneta  at Tue May 10 17:55:10 2005
0014 //
0015 // Last update: Tue May 10 17:55:10 2005
0016 //
0017 #ifndef ROOT_Math_GenVector_EulerAngles
0018 #define ROOT_Math_GenVector_EulerAngles  1
0019
0020 #include "Math/GenVector/Rotation3D.h"
0021 #include "Math/GenVector/DisplacementVector3D.h"
0022 #include "Math/GenVector/PositionVector3D.h"
0023 #include "Math/GenVector/LorentzVector.h"
0024 #include "Math/GenVector/3DConversions.h"
0025 #include <algorithm>
0026 #include <cassert>
0027
0028 namespace ROOT {
0029 namespace Math {
0030
0031
0032 //__________________________________________________________________________________________
0033    /**
0034       EulerAngles class describing rotation as three angles (Euler Angles).
0035       The Euler angles definition matches that of Classical Mechanics (Goldstein).
0036       It is also the same convention defined in
0037       <A HREF="http://mathworld.wolfram.com/EulerAngles.html">mathworld</A>
0038       and used in Mathematica and CLHEP. Note that the ROOT class TRotation defines
0039       a slightly different convention.
0040
0041       @ingroup GenVector
0042
0043       @sa Overview of the @ref GenVector "physics vector library"
0044    */
0045 class EulerAngles {
0046
0047 public:
0048
0049    typedef double Scalar;
0050
0051   /**
0052      Default constructor
0053   */
0054    constexpr EulerAngles() : fPhi(0.0), fTheta(0.0), fPsi(0.0) { }
0055
0056    /**
0057       Constructor from phi, theta and psi
0058    */
0059    EulerAngles( Scalar phi, Scalar theta, Scalar psi ) :
0060       fPhi(phi), fTheta(theta), fPsi(psi)
0061    {Rectify();} // Added 27 Jan. 06   JMM
0062
0063    /**
0064       Construct given a pair of pointers or iterators defining the
0065       beginning and end of an array of three Scalars, to be treated as
0066       the angles phi, theta and psi.
0067    */
0068    template<class IT>
0069    constexpr EulerAngles(IT begin, IT end) { SetComponents(begin,end); }
0070
0071    // The compiler-generated copy ctor, copy assignment, and dtor are OK.
0072
0073    /**
0074       Re-adjust components place angles in canonical ranges
0075    */
0076    void Rectify();
0077
0078
0079    // ======== Construction and assignment from any other rotation ==================
0080
0081    /**
0082       Create from any other supported rotation (see gv_detail::convert )
0083     */
0084    template <class OtherRotation>
0085    explicit constexpr EulerAngles(const OtherRotation & r) {gv_detail::convert(r,*this);}
0086
0087    /**
0088       Assign from any other rotation (see gv_detail::convert )
0089    */
0090    template <class OtherRotation>
0091    EulerAngles &  operator=( OtherRotation const  & r ) {
0092       gv_detail::convert(r,*this);
0093       return *this;
0094    }
0095
0096 #ifdef OLD
0097    explicit EulerAngles(const Rotation3D & r) {gv_detail::convert(r,*this);}
0098
0099    /**
0100       Construct from a rotation matrix
0101    */
0102    explicit EulerAngles(const Rotation3D & r) {gv_detail::convert(r,*this);}
0103
0104    /**
0105       Construct from a rotation represented by a Quaternion
0106    */
0107    explicit EulerAngles(const Quaternion & q) {gv_detail::convert(q,*this);}
0108
0109    /**
0110       Construct from an AxisAngle
0111    */
0112    explicit EulerAngles(const AxisAngle & a ) { gv_detail::convert(a, *this); }
0113
0114    /**
0115       Construct from an axial rotation
0116    */
0117    explicit EulerAngles( RotationZ const & r ) { gv_detail::convert(r, *this); }
0118    explicit EulerAngles( RotationY const & r ) { gv_detail::convert(r, *this); }
0119    explicit EulerAngles( RotationX const & r ) { gv_detail::convert(r, *this); }
0120
0121
0122    /**
0123       Assign from an AxisAngle
0124    */
0125    EulerAngles &
0126    operator=( AxisAngle const & a ) { return operator=(EulerAngles(a)); }
0127
0128    /**
0129       Assign from a Quaternion
0130    */
0131    EulerAngles &
0132    operator=( Quaternion const  & q ) {return operator=(EulerAngles(q)); }
0133
0134    /**
0135       Assign from an axial rotation
0136    */
0137    EulerAngles &
0138    operator=( RotationZ const & r ) { return operator=(EulerAngles(r)); }
0139    EulerAngles &
0140    operator=( RotationY const & r ) { return operator=(EulerAngles(r)); }
0141    EulerAngles &
0142    operator=( RotationX const & r ) { return operator=(EulerAngles(r)); }
0143
0144 #endif
0145
0146    // ======== Components ==============
0147
0148    /**
0149       Set the three Euler angles given a pair of pointers or iterators
0150       defining the beginning and end of an array of three Scalars.
0151    */
0152    template<class IT>
0153    void SetComponents(IT begin, IT end) {
0154       fPhi   = *begin++;
0155       fTheta = *begin++;
0156       fPsi   = *begin++;
0157       (void)end;
0158       assert(begin == end);
0159       Rectify(); // Added 27 Jan. 06   JMM
0160    }
0161
0162    /**
0163       Get the axis and then the angle into data specified by an iterator begin
0164       and another to the end of the desired data (4 past start).
0165    */
0166    template<class IT>
0167    void GetComponents(IT begin, IT end) const {
0168       *begin++ = fPhi;
0169       *begin++ = fTheta;
0170       *begin++ = fPsi;
0171       (void)end;
0172       assert(begin == end);
0173    }
0174
0175    /**
0176       Get the axis and then the angle into data specified by an iterator begin
0177    */
0178    template<class IT>
0179    void GetComponents(IT begin) const {
0180       *begin++ = fPhi;
0181       *begin++ = fTheta;
0182       *begin   = fPsi;
0183    }
0184
0185    /**
0186       Set the components phi, theta, psi based on three Scalars.
0187    */
0188    void SetComponents(Scalar phi, Scalar theta, Scalar psi) {
0189       fPhi=phi; fTheta=theta; fPsi=psi;
0190       Rectify(); // Added 27 Jan. 06   JMM
0191    }
0192
0193    /**
0194       Get the components phi, theta, psi into three Scalars.
0195    */
0196    void GetComponents(Scalar & phi, Scalar & theta, Scalar & psi) const {
0197       phi=fPhi; theta=fTheta; psi=fPsi;
0198    }
0199
0200    /**
0201       Set Phi Euler angle // JMM 30 Jan. 2006
0202    */
0203    void SetPhi(Scalar phi) { fPhi=phi; Rectify(); }
0204
0205    /**
0206       Return Phi Euler angle
0207    */
0208    Scalar Phi() const { return fPhi; }
0209
0210    /**
0211       Set Theta Euler angle // JMM 30 Jan. 2006
0212    */
0213    void SetTheta(Scalar theta) { fTheta=theta; Rectify(); }
0214
0215    /**
0216       Return Theta Euler angle
0217    */
0218    Scalar Theta() const { return fTheta; }
0219
0220    /**
0221       Set Psi Euler angle // JMM 30 Jan. 2006
0222    */
0223    void SetPsi(Scalar psi) { fPsi=psi; Rectify(); }
0224
0225    /**
0226       Return Psi Euler angle
0227    */
0228    Scalar Psi() const { return fPsi; }
0229
0230    // =========== operations ==============
0231
0232
0233    /**
0234       Rotation operation on a displacement vector in any coordinate system and tag
0235    */
0236    template <class CoordSystem, class U>
0237    DisplacementVector3D<CoordSystem,U>
0238    operator() (const DisplacementVector3D<CoordSystem,U> & v) const {
0239       return Rotation3D(*this) ( v );
0240    }
0241
0242    /**
0243       Rotation operation on a position vector in any coordinate system
0244    */
0245    template <class CoordSystem, class U>
0246    PositionVector3D<CoordSystem, U>
0247    operator() (const PositionVector3D<CoordSystem,U> & v) const {
0248       DisplacementVector3D< Cartesian3D<double>,U > xyz(v);
0249       DisplacementVector3D< Cartesian3D<double>,U > rxyz = operator()(xyz);
0250       return PositionVector3D<CoordSystem,U> ( rxyz );
0251    }
0252
0253    /**
0254       Rotation operation on a Lorentz vector in any 4D coordinate system
0255    */
0256    template <class CoordSystem>
0257    LorentzVector<CoordSystem>
0258    operator() (const LorentzVector<CoordSystem> & v) const {
0259       DisplacementVector3D< Cartesian3D<double> > xyz(v.Vect());
0260       xyz = operator()(xyz);
0261       LorentzVector< PxPyPzE4D<double> > xyzt (xyz.X(), xyz.Y(), xyz.Z(), v.E());
0262       return LorentzVector<CoordSystem> ( xyzt );
0263    }
0264
0265    /**
0266       Rotation operation on an arbitrary vector v.
0267       Preconditions:  v must implement methods x(), y(), and z()
0268       and the arbitrary vector type must have a constructor taking (x,y,z)
0269    */
0270    template <class ForeignVector>
0271    ForeignVector
0272    operator() (const  ForeignVector & v) const {
0273       DisplacementVector3D< Cartesian3D<double> > xyz(v);
0274       DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);
0275       return ForeignVector ( rxyz.X(), rxyz.Y(), rxyz.Z() );
0276    }
0277
0278    /**
0279       Overload operator * for rotation on a vector
0280    */
0281    template <class AVector>
0282    inline
0283    AVector operator* (const AVector & v) const
0284    {
0285       return operator()(v);
0286    }
0287
0288    /**
0289       Invert a rotation in place
0290    */
0291    // theta stays the same and negative rotation in Theta is done via a rotation
0292    // of + PI in phi and Psi
0293    void Invert() {
0294       Scalar tmp = -fPhi;
0295       fPhi = -fPsi + Pi();
0296       fPsi=tmp + Pi();
0297    }
0298
0299    /**
0300       Return inverse of a rotation
0301    */
0302    EulerAngles Inverse() const { return EulerAngles(-fPsi + Pi(), fTheta, -fPhi + Pi()); }
0303
0304    // ========= Multi-Rotation Operations ===============
0305
0306    /**
0307       Multiply (combine) two rotations
0308    */
0309    EulerAngles operator * (const Rotation3D  & r) const;
0310    EulerAngles operator * (const AxisAngle   & a) const;
0311    EulerAngles operator * (const EulerAngles & e) const;
0312    EulerAngles operator * (const Quaternion  & q) const;
0313    EulerAngles operator * (const RotationX  & rx) const;
0314    EulerAngles operator * (const RotationY  & ry) const;
0315    EulerAngles operator * (const RotationZ  & rz) const;
0316
0317    /**
0318       Post-Multiply (on right) by another rotation :  T = T*R
0319    */
0320    template <class R>
0321    EulerAngles & operator *= (const R & r) { return *this = (*this)*r; }
0322
0323    /**
0324       Distance between two rotations
0325    */
0326    template <class R>
0327    Scalar Distance ( const R & r ) const {return gv_detail::dist(*this,r);}
0328
0329    /**
0330       Equality/inequality operators
0331    */
0332    bool operator == (const EulerAngles & rhs) const {
0333       if( fPhi   != rhs.fPhi   ) return false;
0334       if( fTheta != rhs.fTheta ) return false;
0335       if( fPsi   != rhs.fPsi   ) return false;
0336       return true;
0337    }
0338    bool operator != (const EulerAngles & rhs) const {
0339       return ! operator==(rhs);
0340    }
0341
0342 private:
0343
0344    double fPhi;     // Z rotation angle (first)  defined in [-PI,PI]
0345    double fTheta;   // X rotation angle (second) defined only [0,PI]
0346    double fPsi;     // Z rotation angle (third)  defined in [-PI,PI]
0347
0348    static double Pi() { return M_PI; }
0349
0350 };  // EulerAngles
0351
0352 /**
0353    Distance between two rotations
0354  */
0355 template <class R>
0356 inline
0357 typename EulerAngles::Scalar
0358 Distance ( const EulerAngles& r1, const R & r2) {return gv_detail::dist(r1,r2);}
0359
0360 /**
0361    Multiplication of an axial rotation by an AxisAngle
0362  */
0363 EulerAngles operator* (RotationX const & r1, EulerAngles const & r2);
0364 EulerAngles operator* (RotationY const & r1, EulerAngles const & r2);
0365 EulerAngles operator* (RotationZ const & r1, EulerAngles const & r2);
0366
0367 /**
0368    Stream Output and Input
0369  */
0370   // TODO - I/O should be put in the manipulator form
0371
0372 std::ostream & operator<< (std::ostream & os, const EulerAngles & e);
0373
0374 } // namespace Math
0375 } // namespace ROOT
0376
0377
0378 #endif /* ROOT_Math_GenVector_EulerAngles  */