Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:13

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 Transform3D
0012 //
0013 // Created by: Lorenzo Moneta  October 21 2005
0014 //
0015 //
0016 #ifndef ROOT_Math_GenVector_Transform3D
0017 #define ROOT_Math_GenVector_Transform3D  1
0018 
0019 
0020 
0021 #include "Math/GenVector/DisplacementVector3D.h"
0022 
0023 #include "Math/GenVector/PositionVector3D.h"
0024 
0025 #include "Math/GenVector/Rotation3D.h"
0026 
0027 #include "Math/GenVector/Translation3D.h"
0028 
0029 
0030 #include "Math/GenVector/AxisAnglefwd.h"
0031 #include "Math/GenVector/EulerAnglesfwd.h"
0032 #include "Math/GenVector/Quaternionfwd.h"
0033 #include "Math/GenVector/RotationZYXfwd.h"
0034 #include "Math/GenVector/RotationXfwd.h"
0035 #include "Math/GenVector/RotationYfwd.h"
0036 #include "Math/GenVector/RotationZfwd.h"
0037 
0038 #include <iostream>
0039 #include <type_traits>
0040 #include <cmath>
0041 
0042 //#include "Math/Vector3Dfwd.h"
0043 
0044 
0045 
0046 namespace ROOT {
0047 
0048 namespace Math {
0049 
0050 namespace Impl {
0051 
0052 //_________________________________________________________________________________________
0053 /**
0054     Basic 3D Transformation class describing  a rotation and then a translation
0055     The internal data are a 3D rotation data (represented as a 3x3 matrix) and a 3D vector data.
0056     They are represented and held in this class like a 3x4 matrix (a simple array of 12 numbers).
0057 
0058     The class can be constructed from any 3D rotation object
0059     (ROOT::Math::Rotation3D, ROOT::Math::AxisAngle, ROOT::Math::Quaternion, etc...) and/or
0060     a 3D Vector (ROOT::Math::DislacementVector3D or via ROOT::Math::Translation ) representing a Translation.
0061     The Transformation is defined by applying first the rotation and then the translation.
0062     A transformation defined by applying first a translation and then a rotation is equivalent to the
0063     transformation obtained applying first the rotation and then a translation equivalent to the rotated vector.
0064     The operator * can be used to obtain directly such transformations, in addition to combine various
0065     transformations.
0066     Keep in mind that the operator * (like in the case of rotations ) is not commutative.
0067     The operator * is used (in addition to operator() ) to apply a transformations on the vector
0068     (DisplacementVector3D and LorentzVector classes) and point (PositionVector3D)  classes.
0069     In the case of Vector objects the transformation only rotates them and does not translate them.
0070     Only Point objects are able to be both rotated and translated.
0071 
0072 
0073     @ingroup GenVector
0074 
0075     @sa Overview of the @ref GenVector "physics vector library"
0076 
0077 */
0078 
0079 template <typename T = double>
0080 class Transform3D {
0081 
0082 public:
0083    typedef T Scalar;
0084 
0085    typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
0086    typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag>     Point;
0087 
0088    enum ETransform3DMatrixIndex {
0089       kXX = 0, kXY = 1, kXZ = 2, kDX = 3,
0090       kYX = 4, kYY = 5, kYZ = 6, kDY = 7,
0091       kZX = 8, kZY = 9, kZZ =10, kDZ = 11
0092    };
0093 
0094 
0095 
0096    /**
0097        Default constructor (identy rotation) + zero translation
0098    */
0099    Transform3D()
0100    {
0101       SetIdentity();
0102    }
0103 
0104    /**
0105       Construct given a pair of pointers or iterators defining the
0106       beginning and end of an array of 12 Scalars
0107    */
0108    template<class IT>
0109    Transform3D(IT begin, IT end)
0110    {
0111       SetComponents(begin,end);
0112    }
0113 
0114    /**
0115       Construct from a rotation and then a translation described by a Vector
0116    */
0117    Transform3D( const Rotation3D & r, const Vector & v)
0118    {
0119       AssignFrom( r, v );
0120    }
0121    /**
0122       Construct from a rotation and then a translation described by a Translation3D class
0123    */
0124    Transform3D(const Rotation3D &r, const Translation3D<T> &t) { AssignFrom(r, t.Vect()); }
0125 
0126    /**
0127       Construct from a rotation (any rotation object)  and then a translation
0128       (represented by any DisplacementVector)
0129       The requirements on the rotation and vector objects are that they can be transformed in a
0130       Rotation3D class and in a Cartesian3D Vector
0131    */
0132    template <class ARotation, class CoordSystem, class Tag>
0133    Transform3D( const ARotation & r, const DisplacementVector3D<CoordSystem,Tag> & v)
0134    {
0135       AssignFrom( Rotation3D(r), Vector (v.X(),v.Y(),v.Z()) );
0136    }
0137 
0138    /**
0139       Construct from a rotation (any rotation object)  and then a translation
0140       represented by a Translation3D class
0141       The requirements on the rotation is that it can be transformed in a
0142       Rotation3D class
0143    */
0144    template <class ARotation>
0145    Transform3D(const ARotation &r, const Translation3D<T> &t)
0146    {
0147       AssignFrom( Rotation3D(r), t.Vect() );
0148    }
0149 
0150 
0151 #ifdef OLD_VERSION
0152    /**
0153       Construct from a translation and then a rotation (inverse assignment)
0154    */
0155    Transform3D( const Vector & v, const Rotation3D & r)
0156    {
0157       // is equivalent from having first the rotation and then the translation vector rotated
0158       AssignFrom( r, r(v) );
0159    }
0160 #endif
0161 
0162    /**
0163       Construct from a 3D Rotation only with zero translation
0164    */
0165    explicit constexpr Transform3D( const Rotation3D & r) {
0166       AssignFrom(r);
0167    }
0168 
0169    // convenience methods for constructing a Transform3D from all the 3D rotations classes
0170    // (cannot use templates for conflict with LA)
0171 
0172    explicit constexpr Transform3D( const AxisAngle & r) {
0173       AssignFrom(Rotation3D(r));
0174    }
0175    explicit constexpr Transform3D( const EulerAngles & r) {
0176       AssignFrom(Rotation3D(r));
0177    }
0178    explicit constexpr Transform3D( const Quaternion & r) {
0179       AssignFrom(Rotation3D(r));
0180    }
0181    explicit constexpr Transform3D( const RotationZYX & r) {
0182       AssignFrom(Rotation3D(r));
0183    }
0184 
0185    // Constructors from axial rotations
0186    // TO DO: implement direct methods for axial rotations without going through Rotation3D
0187    explicit constexpr Transform3D( const RotationX & r) {
0188       AssignFrom(Rotation3D(r));
0189    }
0190    explicit constexpr Transform3D( const RotationY & r) {
0191       AssignFrom(Rotation3D(r));
0192    }
0193    explicit constexpr Transform3D( const RotationZ & r) {
0194       AssignFrom(Rotation3D(r));
0195    }
0196 
0197    /**
0198       Construct from a translation only, represented by any DisplacementVector3D
0199       and with an identity rotation
0200    */
0201    template<class CoordSystem, class Tag>
0202    explicit constexpr Transform3D( const DisplacementVector3D<CoordSystem,Tag> & v) {
0203       AssignFrom(Vector(v.X(),v.Y(),v.Z()));
0204    }
0205    /**
0206       Construct from a translation only, represented by a Cartesian 3D Vector,
0207       and with an identity rotation
0208    */
0209    explicit constexpr Transform3D( const Vector & v) {
0210       AssignFrom(v);
0211    }
0212    /**
0213       Construct from a translation only, represented by a Translation3D class
0214       and with an identity rotation
0215    */
0216    explicit constexpr Transform3D(const Translation3D<T> &t) { AssignFrom(t.Vect()); }
0217 
0218    //#if !defined(__MAKECINT__) && !defined(G__DICTIONARY)  // this is ambiguous with double * , double *
0219 
0220 
0221 #ifdef OLD_VERSION
0222    /**
0223       Construct from a translation (using any type of DisplacementVector )
0224       and then a rotation (any rotation object).
0225       Requirement on the rotation and vector objects are that they can be transformed in a
0226       Rotation3D class and in a Vector
0227    */
0228    template <class ARotation, class CoordSystem, class Tag>
0229    Transform3D(const DisplacementVector3D<CoordSystem,Tag> & v , const ARotation & r)
0230    {
0231       // is equivalent from having first the rotation and then the translation vector rotated
0232       Rotation3D r3d(r);
0233       AssignFrom( r3d, r3d( Vector(v.X(),v.Y(),v.Z()) ) );
0234    }
0235 #endif
0236 
0237 public:
0238    /**
0239       Construct transformation from one coordinate system defined by three
0240       points (origin + two axis) to
0241       a new coordinate system defined by other three points (origin + axis)
0242       Scalar version.
0243       @param fr0  point defining origin of original reference system
0244       @param fr1  point defining first axis of original reference system
0245       @param fr2  point defining second axis of original reference system
0246       @param to0  point defining origin of transformed reference system
0247       @param to1  point defining first axis transformed reference system
0248       @param to2  point defining second axis transformed reference system
0249    */
0250    template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0251    Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0252                const Point &to2)
0253    {
0254       // takes impl. from CLHEP ( E.Chernyaev). To be checked
0255 
0256       Vector x1 = (fr1 - fr0).Unit();
0257       Vector y1 = (fr2 - fr0).Unit();
0258       Vector x2 = (to1 - to0).Unit();
0259       Vector y2 = (to2 - to0).Unit();
0260 
0261       //   C H E C K   A N G L E S
0262 
0263       const T cos1 = x1.Dot(y1);
0264       const T cos2 = x2.Dot(y2);
0265 
0266       if (std::fabs(T(1) - cos1) <= T(0.000001) || std::fabs(T(1) - cos2) <= T(0.000001)) {
0267          std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0268          SetIdentity();
0269       } else {
0270          if (std::fabs(cos1 - cos2) > T(0.000001)) {
0271             std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0272          }
0273 
0274          //   F I N D   R O T A T I O N   M A T R I X
0275 
0276          Vector z1 = (x1.Cross(y1)).Unit();
0277          y1        = z1.Cross(x1);
0278 
0279          Vector z2 = (x2.Cross(y2)).Unit();
0280          y2        = z2.Cross(x2);
0281 
0282          T x1x = x1.x();
0283          T x1y = x1.y();
0284          T x1z = x1.z();
0285          T y1x = y1.x();
0286          T y1y = y1.y();
0287          T y1z = y1.z();
0288          T z1x = z1.x();
0289          T z1y = z1.y();
0290          T z1z = z1.z();
0291 
0292          T x2x = x2.x();
0293          T x2y = x2.y();
0294          T x2z = x2.z();
0295          T y2x = y2.x();
0296          T y2y = y2.y();
0297          T y2z = y2.z();
0298          T z2x = z2.x();
0299          T z2y = z2.y();
0300          T z2z = z2.z();
0301 
0302          T detxx = (y1y * z1z - z1y * y1z);
0303          T detxy = -(y1x * z1z - z1x * y1z);
0304          T detxz = (y1x * z1y - z1x * y1y);
0305          T detyx = -(x1y * z1z - z1y * x1z);
0306          T detyy = (x1x * z1z - z1x * x1z);
0307          T detyz = -(x1x * z1y - z1x * x1y);
0308          T detzx = (x1y * y1z - y1y * x1z);
0309          T detzy = -(x1x * y1z - y1x * x1z);
0310          T detzz = (x1x * y1y - y1x * x1y);
0311 
0312          T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0313          T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0314          T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0315          T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0316          T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0317          T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0318          T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0319          T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0320          T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0321 
0322          //   S E T    T R A N S F O R M A T I O N
0323 
0324          T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0325          T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0326 
0327          SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0328                        dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0329       }
0330    }
0331 
0332    /**
0333       Construct transformation from one coordinate system defined by three
0334       points (origin + two axis) to
0335       a new coordinate system defined by other three points (origin + axis)
0336       Vectorised version.
0337       @param fr0  point defining origin of original reference system
0338       @param fr1  point defining first axis of original reference system
0339       @param fr2  point defining second axis of original reference system
0340       @param to0  point defining origin of transformed reference system
0341       @param to1  point defining first axis transformed reference system
0342       @param to2  point defining second axis transformed reference system
0343    */
0344    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0345    Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0346                const Point &to2)
0347    {
0348       // takes impl. from CLHEP ( E.Chernyaev). To be checked
0349 
0350       Vector x1 = (fr1 - fr0).Unit();
0351       Vector y1 = (fr2 - fr0).Unit();
0352       Vector x2 = (to1 - to0).Unit();
0353       Vector y2 = (to2 - to0).Unit();
0354 
0355       //   C H E C K   A N G L E S
0356 
0357       const T cos1 = x1.Dot(y1);
0358       const T cos2 = x2.Dot(y2);
0359 
0360       const auto m1 = (abs(T(1) - cos1) <= T(0.000001) || abs(T(1) - cos2) <= T(0.000001));
0361 
0362       const auto m2 = (abs(cos1 - cos2) > T(0.000001));
0363       if (any_of(m2)) {
0364          std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0365       }
0366 
0367       //   F I N D   R O T A T I O N   M A T R I X
0368 
0369       Vector z1 = (x1.Cross(y1)).Unit();
0370       y1        = z1.Cross(x1);
0371 
0372       Vector z2 = (x2.Cross(y2)).Unit();
0373       y2        = z2.Cross(x2);
0374 
0375       T x1x = x1.x();
0376       T x1y = x1.y();
0377       T x1z = x1.z();
0378       T y1x = y1.x();
0379       T y1y = y1.y();
0380       T y1z = y1.z();
0381       T z1x = z1.x();
0382       T z1y = z1.y();
0383       T z1z = z1.z();
0384 
0385       T x2x = x2.x();
0386       T x2y = x2.y();
0387       T x2z = x2.z();
0388       T y2x = y2.x();
0389       T y2y = y2.y();
0390       T y2z = y2.z();
0391       T z2x = z2.x();
0392       T z2y = z2.y();
0393       T z2z = z2.z();
0394 
0395       T detxx = (y1y * z1z - z1y * y1z);
0396       T detxy = -(y1x * z1z - z1x * y1z);
0397       T detxz = (y1x * z1y - z1x * y1y);
0398       T detyx = -(x1y * z1z - z1y * x1z);
0399       T detyy = (x1x * z1z - z1x * x1z);
0400       T detyz = -(x1x * z1y - z1x * x1y);
0401       T detzx = (x1y * y1z - y1y * x1z);
0402       T detzy = -(x1x * y1z - y1x * x1z);
0403       T detzz = (x1x * y1y - y1x * x1y);
0404 
0405       T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0406       T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0407       T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0408       T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0409       T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0410       T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0411       T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0412       T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0413       T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0414 
0415       //   S E T    T R A N S F O R M A T I O N
0416 
0417       T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0418       T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0419 
0420       SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0421                     dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0422 
0423       if (any_of(m1)) {
0424          std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0425          SetIdentity(m1);
0426       }
0427    }
0428 
0429    // use compiler generated copy ctor, copy assignment and destructor
0430 
0431    /**
0432       Construct from a linear algebra matrix of size at least 3x4,
0433       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0434       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0435       are described by the 4-th column
0436    */
0437    template<class ForeignMatrix>
0438    explicit constexpr Transform3D(const ForeignMatrix & m) {
0439       SetComponents(m);
0440    }
0441 
0442    /**
0443       Raw constructor from 12 Scalar components
0444    */
0445    Transform3D(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0446    {
0447       SetComponents (xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz);
0448    }
0449 
0450 
0451    /**
0452       Construct from a linear algebra matrix of size at least 3x4,
0453       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0454       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0455       are described by the 4-th column
0456    */
0457    template <class ForeignMatrix>
0458    Transform3D<T> &operator=(const ForeignMatrix &m)
0459    {
0460       SetComponents(m);
0461       return *this;
0462    }
0463 
0464 
0465    // ======== Components ==============
0466 
0467 
0468    /**
0469       Set the 12 matrix components given an iterator to the start of
0470       the desired data, and another to the end (12 past start).
0471    */
0472    template<class IT>
0473    void SetComponents(IT begin, IT end) {
0474       for (int i = 0; i <12; ++i) {
0475          fM[i] = *begin;
0476          ++begin;
0477       }
0478       (void)end;
0479       assert (end==begin);
0480    }
0481 
0482    /**
0483       Get the 12 matrix components into data specified by an iterator begin
0484       and another to the end of the desired data (12 past start).
0485    */
0486    template<class IT>
0487    void GetComponents(IT begin, IT end) const {
0488       for (int i = 0; i <12; ++i) {
0489          *begin = fM[i];
0490          ++begin;
0491       }
0492       (void)end;
0493       assert (end==begin);
0494    }
0495 
0496    /**
0497       Get the 12 matrix components into data specified by an iterator begin
0498    */
0499    template<class IT>
0500    void GetComponents(IT begin) const {
0501       std::copy(fM, fM + 12, begin);
0502    }
0503 
0504    /**
0505       Set components from a linear algebra matrix of size at least 3x4,
0506       which must support operator()(i,j) to obtain elements (0,0) thru (2,3).
0507       The 3x3 sub-block is assumed to be the rotation part and the translations vector
0508       are described by the 4-th column
0509    */
0510    template<class ForeignMatrix>
0511    void
0512    SetTransformMatrix (const ForeignMatrix & m) {
0513       fM[kXX]=m(0,0);  fM[kXY]=m(0,1);  fM[kXZ]=m(0,2); fM[kDX]=m(0,3);
0514       fM[kYX]=m(1,0);  fM[kYY]=m(1,1);  fM[kYZ]=m(1,2); fM[kDY]=m(1,3);
0515       fM[kZX]=m(2,0);  fM[kZY]=m(2,1);  fM[kZZ]=m(2,2); fM[kDZ]=m(2,3);
0516    }
0517 
0518    /**
0519       Get components into a linear algebra matrix of size at least 3x4,
0520       which must support operator()(i,j) for write access to elements
0521       (0,0) thru (2,3).
0522    */
0523    template<class ForeignMatrix>
0524    void
0525    GetTransformMatrix (ForeignMatrix & m) const {
0526       m(0,0)=fM[kXX];  m(0,1)=fM[kXY];  m(0,2)=fM[kXZ];  m(0,3)=fM[kDX];
0527       m(1,0)=fM[kYX];  m(1,1)=fM[kYY];  m(1,2)=fM[kYZ];  m(1,3)=fM[kDY];
0528       m(2,0)=fM[kZX];  m(2,1)=fM[kZY];  m(2,2)=fM[kZZ];  m(2,3)=fM[kDZ];
0529    }
0530 
0531 
0532    /**
0533       Set the components from 12 scalars
0534    */
0535    void SetComponents(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0536    {
0537       fM[kXX]=xx;  fM[kXY]=xy;  fM[kXZ]=xz;  fM[kDX]=dx;
0538       fM[kYX]=yx;  fM[kYY]=yy;  fM[kYZ]=yz;  fM[kDY]=dy;
0539       fM[kZX]=zx;  fM[kZY]=zy;  fM[kZZ]=zz;  fM[kDZ]=dz;
0540    }
0541 
0542    /**
0543       Get the components into 12 scalars
0544    */
0545    void GetComponents(T &xx, T &xy, T &xz, T &dx, T &yx, T &yy, T &yz, T &dy, T &zx, T &zy, T &zz, T &dz) const
0546    {
0547       xx=fM[kXX];  xy=fM[kXY];  xz=fM[kXZ];  dx=fM[kDX];
0548       yx=fM[kYX];  yy=fM[kYY];  yz=fM[kYZ];  dy=fM[kDY];
0549       zx=fM[kZX];  zy=fM[kZY];  zz=fM[kZZ];  dz=fM[kDZ];
0550    }
0551 
0552 
0553    /**
0554       Get the rotation and translation vector representing the 3D transformation
0555       in any rotation and any vector (the Translation class could also be used)
0556    */
0557    template<class AnyRotation, class V>
0558    void GetDecomposition(AnyRotation &r, V &v) const {
0559       GetRotation(r);
0560       GetTranslation(v);
0561    }
0562 
0563 
0564    /**
0565       Get the rotation and translation vector representing the 3D transformation
0566    */
0567    void GetDecomposition(Rotation3D &r, Vector &v) const {
0568       GetRotation(r);
0569       GetTranslation(v);
0570    }
0571 
0572    /**
0573       Get the 3D rotation representing the 3D transformation
0574    */
0575    Rotation3D Rotation() const {
0576       return Rotation3D( fM[kXX], fM[kXY], fM[kXZ],
0577                          fM[kYX], fM[kYY], fM[kYZ],
0578                          fM[kZX], fM[kZY], fM[kZZ] );
0579    }
0580 
0581    /**
0582       Get the rotation representing the 3D transformation
0583    */
0584    template <class AnyRotation>
0585    AnyRotation Rotation() const {
0586       return AnyRotation(Rotation3D(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]));
0587    }
0588 
0589    /**
0590       Get the  rotation (any type) representing the 3D transformation
0591    */
0592    template <class AnyRotation>
0593    void GetRotation(AnyRotation &r) const {
0594       r = Rotation();
0595    }
0596 
0597    /**
0598       Get the translation representing the 3D transformation in a Cartesian vector
0599    */
0600    Translation3D<T> Translation() const { return Translation3D<T>(fM[kDX], fM[kDY], fM[kDZ]); }
0601 
0602    /**
0603       Get the translation representing the 3D transformation in any vector
0604       which implements the SetXYZ method
0605    */
0606    template <class AnyVector>
0607    void GetTranslation(AnyVector &v) const {
0608       v.SetXYZ(fM[kDX], fM[kDY], fM[kDZ]);
0609    }
0610 
0611 
0612 
0613    // operations on points and vectors
0614 
0615    /**
0616       Transformation operation for Position Vector in Cartesian coordinate
0617       For a Position Vector first a rotation and then a translation is applied
0618    */
0619    Point operator() (const Point & p) const {
0620       return Point ( fM[kXX]*p.X() + fM[kXY]*p.Y() + fM[kXZ]*p.Z() + fM[kDX],
0621                      fM[kYX]*p.X() + fM[kYY]*p.Y() + fM[kYZ]*p.Z() + fM[kDY],
0622                      fM[kZX]*p.X() + fM[kZY]*p.Y() + fM[kZZ]*p.Z() + fM[kDZ] );
0623    }
0624 
0625 
0626    /**
0627       Transformation operation for Displacement Vectors in Cartesian coordinate
0628       For the Displacement Vectors only the rotation applies - no translations
0629    */
0630    Vector operator() (const Vector & v) const {
0631       return Vector( fM[kXX]*v.X() + fM[kXY]*v.Y() + fM[kXZ]*v.Z() ,
0632                      fM[kYX]*v.X() + fM[kYY]*v.Y() + fM[kYZ]*v.Z() ,
0633                      fM[kZX]*v.X() + fM[kZY]*v.Y() + fM[kZZ]*v.Z()  );
0634    }
0635 
0636 
0637    /**
0638       Transformation operation for Position Vector in any coordinate system
0639    */
0640    template <class CoordSystem>
0641    PositionVector3D<CoordSystem> operator()(const PositionVector3D<CoordSystem> &p) const
0642    {
0643       return PositionVector3D<CoordSystem>(operator()(Point(p)));
0644    }
0645    /**
0646       Transformation operation for Position Vector in any coordinate system
0647    */
0648    template <class CoordSystem>
0649    PositionVector3D<CoordSystem> operator*(const PositionVector3D<CoordSystem> &v) const
0650    {
0651       return operator()(v);
0652    }
0653 
0654    /**
0655       Transformation operation for Displacement Vector in any coordinate system
0656    */
0657    template<class CoordSystem >
0658    DisplacementVector3D<CoordSystem> operator() (const DisplacementVector3D <CoordSystem> & v) const {
0659       return DisplacementVector3D<CoordSystem>(operator()(Vector(v)));
0660    }
0661    /**
0662       Transformation operation for Displacement Vector in any coordinate system
0663    */
0664    template <class CoordSystem>
0665    DisplacementVector3D<CoordSystem> operator*(const DisplacementVector3D<CoordSystem> &v) const
0666    {
0667       return operator()(v);
0668    }
0669 
0670    /**
0671       Directly apply the inverse affine transformation on vectors.
0672       Avoids having to calculate the inverse as an intermediate result.
0673       This is possible since the inverse of a rotation is its transpose.
0674    */
0675    Vector ApplyInverse(const Vector &v) const
0676    {
0677       return Vector(fM[kXX] * v.X() + fM[kYX] * v.Y() + fM[kZX] * v.Z(),
0678                     fM[kXY] * v.X() + fM[kYY] * v.Y() + fM[kZY] * v.Z(),
0679                     fM[kXZ] * v.X() + fM[kYZ] * v.Y() + fM[kZZ] * v.Z());
0680    }
0681 
0682    /**
0683       Directly apply the inverse affine transformation on points
0684       (first inverse translation then inverse rotation).
0685       Avoids having to calculate the inverse as an intermediate result.
0686       This is possible since the inverse of a rotation is its transpose.
0687    */
0688    Point ApplyInverse(const Point &p) const
0689    {
0690       Point tmp(p.X() - fM[kDX], p.Y() - fM[kDY], p.Z() - fM[kDZ]);
0691       return Point(fM[kXX] * tmp.X() + fM[kYX] * tmp.Y() + fM[kZX] * tmp.Z(),
0692                    fM[kXY] * tmp.X() + fM[kYY] * tmp.Y() + fM[kZY] * tmp.Z(),
0693                    fM[kXZ] * tmp.X() + fM[kYZ] * tmp.Y() + fM[kZZ] * tmp.Z());
0694    }
0695 
0696    /**
0697       Directly apply the inverse affine transformation on an arbitrary
0698       coordinate-system point.
0699       Involves casting to Point(p) type.
0700    */
0701    template <class CoordSystem>
0702    PositionVector3D<CoordSystem> ApplyInverse(const PositionVector3D<CoordSystem> &p) const
0703    {
0704       return PositionVector3D<CoordSystem>(ApplyInverse(Point(p)));
0705    }
0706 
0707    /**
0708       Directly apply the inverse affine transformation on an arbitrary
0709       coordinate-system vector.
0710       Involves casting to Vector(p) type.
0711    */
0712    template <class CoordSystem>
0713    DisplacementVector3D<CoordSystem> ApplyInverse(const DisplacementVector3D<CoordSystem> &p) const
0714    {
0715       return DisplacementVector3D<CoordSystem>(ApplyInverse(Vector(p)));
0716    }
0717 
0718    /**
0719       Transformation operation for points between different coordinate system tags
0720    */
0721    template <class CoordSystem, class Tag1, class Tag2>
0722    void Transform(const PositionVector3D<CoordSystem, Tag1> &p1, PositionVector3D<CoordSystem, Tag2> &p2) const
0723    {
0724       const Point xyzNew = operator()(Point(p1.X(), p1.Y(), p1.Z()));
0725       p2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0726    }
0727 
0728 
0729    /**
0730       Transformation operation for Displacement Vector of different coordinate systems
0731    */
0732    template <class CoordSystem, class Tag1, class Tag2>
0733    void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1, DisplacementVector3D<CoordSystem, Tag2> &v2) const
0734    {
0735       const Vector xyzNew = operator()(Vector(v1.X(), v1.Y(), v1.Z()));
0736       v2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0737    }
0738 
0739    /**
0740       Transformation operation for a Lorentz Vector in any  coordinate system
0741    */
0742    template <class CoordSystem >
0743    LorentzVector<CoordSystem> operator() (const LorentzVector<CoordSystem> & q) const {
0744       const Vector xyzNew = operator()(Vector(q.Vect()));
0745       return LorentzVector<CoordSystem>(xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E());
0746    }
0747    /**
0748       Transformation operation for a Lorentz Vector in any  coordinate system
0749    */
0750    template <class CoordSystem>
0751    LorentzVector<CoordSystem> operator*(const LorentzVector<CoordSystem> &q) const
0752    {
0753       return operator()(q);
0754    }
0755 
0756    /**
0757       Transformation on a 3D plane
0758    */
0759    template <typename TYPE>
0760    Plane3D<TYPE> operator()(const Plane3D<TYPE> &plane) const
0761    {
0762       // transformations on a 3D plane
0763       const auto n = plane.Normal();
0764       // take a point on the plane. Use origin projection on the plane
0765       // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
0766       const auto d = plane.HesseDistance();
0767       Point p(-d * n.X(), -d * n.Y(), -d * n.Z());
0768       return Plane3D<TYPE>(operator()(n), operator()(p));
0769    }
0770 
0771    /// Multiplication operator for 3D plane
0772    template <typename TYPE>
0773    Plane3D<TYPE> operator*(const Plane3D<TYPE> &plane) const
0774    {
0775       return operator()(plane);
0776    }
0777 
0778    // skip transformation for arbitrary vectors - not really defined if point or displacement vectors
0779 
0780    /**
0781       multiply (combine) with another transformation in place
0782    */
0783    inline Transform3D<T> &operator*=(const Transform3D<T> &t);
0784 
0785    /**
0786       multiply (combine) two transformations
0787    */
0788    inline Transform3D<T> operator*(const Transform3D<T> &t) const;
0789 
0790    /**
0791        Invert the transformation in place (scalar)
0792    */
0793    template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0794    void Invert()
0795    {
0796       //
0797       // Name: Transform3D::inverse                     Date:    24.09.96
0798       // Author: E.Chernyaev (IHEP/Protvino)            Revised:
0799       //
0800       // Function: Find inverse affine transformation.
0801 
0802       T detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0803       T detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0804       T detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0805       T det   = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0806       if (det == T(0)) {
0807          std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0808          return;
0809       }
0810       det = T(1) / det;
0811       detxx *= det;
0812       detxy *= det;
0813       detxz *= det;
0814       T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0815       T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0816       T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0817       T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0818       T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0819       T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0820       SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0821                     detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0822                     -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0823    }
0824 
0825    /**
0826        Invert the transformation in place (vectorised)
0827    */
0828    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0829    void Invert()
0830    {
0831       //
0832       // Name: Transform3D::inverse                     Date:    24.09.96
0833       // Author: E.Chernyaev (IHEP/Protvino)            Revised:
0834       //
0835       // Function: Find inverse affine transformation.
0836 
0837       T          detxx    = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0838       T          detxy    = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0839       T          detxz    = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0840       T          det      = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0841       const auto detZmask = (det == T(0));
0842       if (any_of(detZmask)) {
0843          std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0844          det(detZmask) = T(1);
0845       }
0846       det = T(1) / det;
0847       detxx *= det;
0848       detxy *= det;
0849       detxz *= det;
0850       T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0851       T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0852       T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0853       T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0854       T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0855       T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0856       // Set det=0 cases to 0
0857       if (any_of(detZmask)) {
0858          detxx(detZmask) = T(0);
0859          detxy(detZmask) = T(0);
0860          detxz(detZmask) = T(0);
0861          detyx(detZmask) = T(0);
0862          detyy(detZmask) = T(0);
0863          detyz(detZmask) = T(0);
0864          detzx(detZmask) = T(0);
0865          detzy(detZmask) = T(0);
0866          detzz(detZmask) = T(0);
0867       }
0868       // set final components
0869       SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0870                     detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0871                     -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0872    }
0873 
0874    /**
0875       Return the inverse of the transformation.
0876    */
0877    Transform3D<T> Inverse() const
0878    {
0879       Transform3D<T> t(*this);
0880       t.Invert();
0881       return t;
0882    }
0883 
0884    /**
0885       Equality operator. Check equality for each element
0886       To do: use T tolerance
0887    */
0888    bool operator==(const Transform3D<T> &rhs) const
0889    {
0890       return (fM[0] == rhs.fM[0] && fM[1] == rhs.fM[1] && fM[2] == rhs.fM[2] && fM[3] == rhs.fM[3] &&
0891               fM[4] == rhs.fM[4] && fM[5] == rhs.fM[5] && fM[6] == rhs.fM[6] && fM[7] == rhs.fM[7] &&
0892               fM[8] == rhs.fM[8] && fM[9] == rhs.fM[9] && fM[10] == rhs.fM[10] && fM[11] == rhs.fM[11]);
0893    }
0894 
0895    /**
0896       Inequality operator. Check equality for each element
0897       To do: use T tolerance
0898    */
0899    bool operator!=(const Transform3D<T> &rhs) const { return !operator==(rhs); }
0900 
0901 protected:
0902 
0903    /**
0904       make transformation from first a rotation then a translation
0905    */
0906    void AssignFrom(const Rotation3D &r, const Vector &v)
0907    {
0908       // assignment  from rotation + translation
0909 
0910       T rotData[9];
0911       r.GetComponents(rotData, rotData + 9);
0912       // first raw
0913       for (int i = 0; i < 3; ++i) fM[i] = rotData[i];
0914       // second raw
0915       for (int i = 0; i < 3; ++i) fM[kYX + i] = rotData[3 + i];
0916       // third raw
0917       for (int i = 0; i < 3; ++i) fM[kZX + i] = rotData[6 + i];
0918 
0919       // translation data
0920       T vecData[3];
0921       v.GetCoordinates(vecData, vecData + 3);
0922       fM[kDX] = vecData[0];
0923       fM[kDY] = vecData[1];
0924       fM[kDZ] = vecData[2];
0925    }
0926 
0927    /**
0928       make transformation from only rotations (zero translation)
0929    */
0930    void AssignFrom(const Rotation3D &r)
0931    {
0932       // assign from only a rotation  (null translation)
0933       T rotData[9];
0934       r.GetComponents(rotData, rotData + 9);
0935       for (int i = 0; i < 3; ++i) {
0936          for (int j = 0; j < 3; ++j) fM[4 * i + j] = rotData[3 * i + j];
0937          // empty vector data
0938          fM[4 * i + 3] = T(0);
0939       }
0940    }
0941 
0942    /**
0943       make transformation from only translation (identity rotations)
0944    */
0945    void AssignFrom(const Vector &v)
0946    {
0947       // assign from a translation only (identity rotations)
0948       fM[kXX] = T(1);
0949       fM[kXY] = T(0);
0950       fM[kXZ] = T(0);
0951       fM[kDX] = v.X();
0952       fM[kYX] = T(0);
0953       fM[kYY] = T(1);
0954       fM[kYZ] = T(0);
0955       fM[kDY] = v.Y();
0956       fM[kZX] = T(0);
0957       fM[kZY] = T(0);
0958       fM[kZZ] = T(1);
0959       fM[kDZ] = v.Z();
0960    }
0961 
0962    /**
0963       Set identity transformation (identity rotation , zero translation)
0964    */
0965    void SetIdentity()
0966    {
0967       // set identity ( identity rotation and zero translation)
0968       fM[kXX] = T(1);
0969       fM[kXY] = T(0);
0970       fM[kXZ] = T(0);
0971       fM[kDX] = T(0);
0972       fM[kYX] = T(0);
0973       fM[kYY] = T(1);
0974       fM[kYZ] = T(0);
0975       fM[kDY] = T(0);
0976       fM[kZX] = T(0);
0977       fM[kZY] = T(0);
0978       fM[kZZ] = T(1);
0979       fM[kDZ] = T(0);
0980    }
0981 
0982    /**
0983       Set identity transformation (identity rotation , zero translation)
0984       vectorised version that sets using a mask
0985    */
0986    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0987    void SetIdentity(const typename SCALAR::mask_type m)
0988    {
0989       // set identity ( identity rotation and zero translation)
0990       fM[kXX](m) = T(1);
0991       fM[kXY](m) = T(0);
0992       fM[kXZ](m) = T(0);
0993       fM[kDX](m) = T(0);
0994       fM[kYX](m) = T(0);
0995       fM[kYY](m) = T(1);
0996       fM[kYZ](m) = T(0);
0997       fM[kDY](m) = T(0);
0998       fM[kZX](m) = T(0);
0999       fM[kZY](m) = T(0);
1000       fM[kZZ](m) = T(1);
1001       fM[kDZ](m) = T(0);
1002    }
1003 
1004 private:
1005    T fM[12]; // transformation elements (3x4 matrix)
1006 };
1007 
1008 
1009 
1010 
1011 // inline functions (combination of transformations)
1012 
1013 template <class T>
1014 inline Transform3D<T> &Transform3D<T>::operator*=(const Transform3D<T> &t)
1015 {
1016    // combination of transformations
1017 
1018    SetComponents(fM[kXX]*t.fM[kXX]+fM[kXY]*t.fM[kYX]+fM[kXZ]*t.fM[kZX],
1019                  fM[kXX]*t.fM[kXY]+fM[kXY]*t.fM[kYY]+fM[kXZ]*t.fM[kZY],
1020                  fM[kXX]*t.fM[kXZ]+fM[kXY]*t.fM[kYZ]+fM[kXZ]*t.fM[kZZ],
1021                  fM[kXX]*t.fM[kDX]+fM[kXY]*t.fM[kDY]+fM[kXZ]*t.fM[kDZ]+fM[kDX],
1022 
1023                  fM[kYX]*t.fM[kXX]+fM[kYY]*t.fM[kYX]+fM[kYZ]*t.fM[kZX],
1024                  fM[kYX]*t.fM[kXY]+fM[kYY]*t.fM[kYY]+fM[kYZ]*t.fM[kZY],
1025                  fM[kYX]*t.fM[kXZ]+fM[kYY]*t.fM[kYZ]+fM[kYZ]*t.fM[kZZ],
1026                  fM[kYX]*t.fM[kDX]+fM[kYY]*t.fM[kDY]+fM[kYZ]*t.fM[kDZ]+fM[kDY],
1027 
1028                  fM[kZX]*t.fM[kXX]+fM[kZY]*t.fM[kYX]+fM[kZZ]*t.fM[kZX],
1029                  fM[kZX]*t.fM[kXY]+fM[kZY]*t.fM[kYY]+fM[kZZ]*t.fM[kZY],
1030                  fM[kZX]*t.fM[kXZ]+fM[kZY]*t.fM[kYZ]+fM[kZZ]*t.fM[kZZ],
1031                  fM[kZX]*t.fM[kDX]+fM[kZY]*t.fM[kDY]+fM[kZZ]*t.fM[kDZ]+fM[kDZ]);
1032 
1033    return *this;
1034 }
1035 
1036 template <class T>
1037 inline Transform3D<T> Transform3D<T>::operator*(const Transform3D<T> &t) const
1038 {
1039    // combination of transformations
1040 
1041    return Transform3D<T>(fM[kXX] * t.fM[kXX] + fM[kXY] * t.fM[kYX] + fM[kXZ] * t.fM[kZX],
1042                          fM[kXX] * t.fM[kXY] + fM[kXY] * t.fM[kYY] + fM[kXZ] * t.fM[kZY],
1043                          fM[kXX] * t.fM[kXZ] + fM[kXY] * t.fM[kYZ] + fM[kXZ] * t.fM[kZZ],
1044                          fM[kXX] * t.fM[kDX] + fM[kXY] * t.fM[kDY] + fM[kXZ] * t.fM[kDZ] + fM[kDX],
1045 
1046                          fM[kYX] * t.fM[kXX] + fM[kYY] * t.fM[kYX] + fM[kYZ] * t.fM[kZX],
1047                          fM[kYX] * t.fM[kXY] + fM[kYY] * t.fM[kYY] + fM[kYZ] * t.fM[kZY],
1048                          fM[kYX] * t.fM[kXZ] + fM[kYY] * t.fM[kYZ] + fM[kYZ] * t.fM[kZZ],
1049                          fM[kYX] * t.fM[kDX] + fM[kYY] * t.fM[kDY] + fM[kYZ] * t.fM[kDZ] + fM[kDY],
1050 
1051                          fM[kZX] * t.fM[kXX] + fM[kZY] * t.fM[kYX] + fM[kZZ] * t.fM[kZX],
1052                          fM[kZX] * t.fM[kXY] + fM[kZY] * t.fM[kYY] + fM[kZZ] * t.fM[kZY],
1053                          fM[kZX] * t.fM[kXZ] + fM[kZY] * t.fM[kYZ] + fM[kZZ] * t.fM[kZZ],
1054                          fM[kZX] * t.fM[kDX] + fM[kZY] * t.fM[kDY] + fM[kZZ] * t.fM[kDZ] + fM[kDZ]);
1055 }
1056 
1057 
1058 
1059 
1060 //--- global functions resulting in Transform3D -------
1061 
1062 
1063 // ------ combination of a  translation (first)  and a rotation ------
1064 
1065 
1066 /**
1067    combine a translation and a rotation to give a transform3d
1068    First the translation then the rotation
1069  */
1070 template <class T>
1071 inline Transform3D<T> operator*(const Rotation3D &r, const Translation3D<T> &t)
1072 {
1073    return Transform3D<T>(r, r(t.Vect()));
1074 }
1075 template <class T>
1076 inline Transform3D<T> operator*(const RotationX &r, const Translation3D<T> &t)
1077 {
1078    Rotation3D r3(r);
1079    return Transform3D<T>(r3, r3(t.Vect()));
1080 }
1081 template <class T>
1082 inline Transform3D<T> operator*(const RotationY &r, const Translation3D<T> &t)
1083 {
1084    Rotation3D r3(r);
1085    return Transform3D<T>(r3, r3(t.Vect()));
1086 }
1087 template <class T>
1088 inline Transform3D<T> operator*(const RotationZ &r, const Translation3D<T> &t)
1089 {
1090    Rotation3D r3(r);
1091    return Transform3D<T>(r3, r3(t.Vect()));
1092 }
1093 template <class T>
1094 inline Transform3D<T> operator*(const RotationZYX &r, const Translation3D<T> &t)
1095 {
1096    Rotation3D r3(r);
1097    return Transform3D<T>(r3, r3(t.Vect()));
1098 }
1099 template <class T>
1100 inline Transform3D<T> operator*(const AxisAngle &r, const Translation3D<T> &t)
1101 {
1102    Rotation3D r3(r);
1103    return Transform3D<T>(r3, r3(t.Vect()));
1104 }
1105 template <class T>
1106 inline Transform3D<T> operator*(const EulerAngles &r, const Translation3D<T> &t)
1107 {
1108    Rotation3D r3(r);
1109    return Transform3D<T>(r3, r3(t.Vect()));
1110 }
1111 template <class T>
1112 inline Transform3D<T> operator*(const Quaternion &r, const Translation3D<T> &t)
1113 {
1114    Rotation3D r3(r);
1115    return Transform3D<T>(r3, r3(t.Vect()));
1116 }
1117 
1118 // ------ combination of a  rotation (first)  and then a translation ------
1119 
1120 /**
1121    combine a rotation and a translation to give a transform3d
1122    First a rotation then the translation
1123  */
1124 template <class T>
1125 inline Transform3D<T> operator*(const Translation3D<T> &t, const Rotation3D &r)
1126 {
1127    return Transform3D<T>(r, t.Vect());
1128 }
1129 template <class T>
1130 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationX &r)
1131 {
1132    return Transform3D<T>(Rotation3D(r), t.Vect());
1133 }
1134 template <class T>
1135 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationY &r)
1136 {
1137    return Transform3D<T>(Rotation3D(r), t.Vect());
1138 }
1139 template <class T>
1140 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZ &r)
1141 {
1142    return Transform3D<T>(Rotation3D(r), t.Vect());
1143 }
1144 template <class T>
1145 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZYX &r)
1146 {
1147    return Transform3D<T>(Rotation3D(r), t.Vect());
1148 }
1149 template <class T>
1150 inline Transform3D<T> operator*(const Translation3D<T> &t, const EulerAngles &r)
1151 {
1152    return Transform3D<T>(Rotation3D(r), t.Vect());
1153 }
1154 template <class T>
1155 inline Transform3D<T> operator*(const Translation3D<T> &t, const Quaternion &r)
1156 {
1157    return Transform3D<T>(Rotation3D(r), t.Vect());
1158 }
1159 template <class T>
1160 inline Transform3D<T> operator*(const Translation3D<T> &t, const AxisAngle &r)
1161 {
1162    return Transform3D<T>(Rotation3D(r), t.Vect());
1163 }
1164 
1165 // ------ combination of a Transform3D and a pure translation------
1166 
1167 /**
1168    combine a transformation and a translation to give a transform3d
1169    First the translation then the transform3D
1170  */
1171 template <class T>
1172 inline Transform3D<T> operator*(const Transform3D<T> &t, const Translation3D<T> &d)
1173 {
1174    Rotation3D r = t.Rotation();
1175    return Transform3D<T>(r, r(d.Vect()) + t.Translation().Vect());
1176 }
1177 
1178 /**
1179    combine a translation and a transformation to give a transform3d
1180    First the transformation then the translation
1181  */
1182 template <class T>
1183 inline Transform3D<T> operator*(const Translation3D<T> &d, const Transform3D<T> &t)
1184 {
1185    return Transform3D<T>(t.Rotation(), t.Translation().Vect() + d.Vect());
1186 }
1187 
1188 // ------ combination of a Transform3D and any rotation------
1189 
1190 
1191 /**
1192    combine a transformation and a rotation to give a transform3d
1193    First the rotation then the transform3D
1194  */
1195 template <class T>
1196 inline Transform3D<T> operator*(const Transform3D<T> &t, const Rotation3D &r)
1197 {
1198    return Transform3D<T>(t.Rotation() * r, t.Translation());
1199 }
1200 template <class T>
1201 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationX &r)
1202 {
1203    return Transform3D<T>(t.Rotation() * r, t.Translation());
1204 }
1205 template <class T>
1206 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationY &r)
1207 {
1208    return Transform3D<T>(t.Rotation() * r, t.Translation());
1209 }
1210 template <class T>
1211 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZ &r)
1212 {
1213    return Transform3D<T>(t.Rotation() * r, t.Translation());
1214 }
1215 template <class T>
1216 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZYX &r)
1217 {
1218    return Transform3D<T>(t.Rotation() * r, t.Translation());
1219 }
1220 template <class T>
1221 inline Transform3D<T> operator*(const Transform3D<T> &t, const EulerAngles &r)
1222 {
1223    return Transform3D<T>(t.Rotation() * r, t.Translation());
1224 }
1225 template <class T>
1226 inline Transform3D<T> operator*(const Transform3D<T> &t, const AxisAngle &r)
1227 {
1228    return Transform3D<T>(t.Rotation() * r, t.Translation());
1229 }
1230 template <class T>
1231 inline Transform3D<T> operator*(const Transform3D<T> &t, const Quaternion &r)
1232 {
1233    return Transform3D<T>(t.Rotation() * r, t.Translation());
1234 }
1235 
1236 
1237 
1238 /**
1239    combine a rotation and a transformation to give a transform3d
1240    First the transformation then the rotation
1241  */
1242 template <class T>
1243 inline Transform3D<T> operator*(const Rotation3D &r, const Transform3D<T> &t)
1244 {
1245    return Transform3D<T>(r * t.Rotation(), r * t.Translation().Vect());
1246 }
1247 template <class T>
1248 inline Transform3D<T> operator*(const RotationX &r, const Transform3D<T> &t)
1249 {
1250    Rotation3D r3d(r);
1251    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1252 }
1253 template <class T>
1254 inline Transform3D<T> operator*(const RotationY &r, const Transform3D<T> &t)
1255 {
1256    Rotation3D r3d(r);
1257    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1258 }
1259 template <class T>
1260 inline Transform3D<T> operator*(const RotationZ &r, const Transform3D<T> &t)
1261 {
1262    Rotation3D r3d(r);
1263    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1264 }
1265 template <class T>
1266 inline Transform3D<T> operator*(const RotationZYX &r, const Transform3D<T> &t)
1267 {
1268    Rotation3D r3d(r);
1269    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1270 }
1271 template <class T>
1272 inline Transform3D<T> operator*(const EulerAngles &r, const Transform3D<T> &t)
1273 {
1274    Rotation3D r3d(r);
1275    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1276 }
1277 template <class T>
1278 inline Transform3D<T> operator*(const AxisAngle &r, const Transform3D<T> &t)
1279 {
1280    Rotation3D r3d(r);
1281    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1282 }
1283 template <class T>
1284 inline Transform3D<T> operator*(const Quaternion &r, const Transform3D<T> &t)
1285 {
1286    Rotation3D r3d(r);
1287    return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1288 }
1289 
1290 
1291 //---I/O functions
1292 // TODO - I/O should be put in the manipulator form
1293 
1294 /**
1295    print the 12 components of the Transform3D
1296  */
1297 template <class T>
1298 std::ostream &operator<<(std::ostream &os, const Transform3D<T> &t)
1299 {
1300    // TODO - this will need changing for machine-readable issues
1301    //        and even the human readable form needs formatting improvements
1302 
1303    T m[12];
1304    t.GetComponents(m, m + 12);
1305    os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "  " << m[3];
1306    os << "\n" << m[4] << "  " << m[5] << "  " << m[6] << "  " << m[7];
1307    os << "\n" << m[8] << "  " << m[9] << "  " << m[10] << "  " << m[11] << "\n";
1308    return os;
1309 }
1310 
1311 } // end namespace Impl
1312 
1313 // typedefs for double and float versions
1314 typedef Impl::Transform3D<double> Transform3D;
1315 typedef Impl::Transform3D<float>  Transform3DF;
1316 
1317 } // end namespace Math
1318 
1319 } // end namespace ROOT
1320 
1321 
1322 #endif /* ROOT_Math_GenVector_Transform3D */