Math/GenVector/RotationZYX.h

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