TGeoMatrix.cxx

Go to the documentation of this file.

// @(#)root/geom:$Id$
// Author: Andrei Gheata   25/10/01

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

/** \class TGeoMatrix
\ingroup Geometry_classes

Geometrical transformation package.

  All geometrical transformations handled by the modeller are provided as a
built-in package. This was designed to minimize memory requirements and
optimize performance of point/vector master-to-local and local-to-master
computation. We need to have in mind that a transformation in TGeo has 2
major use-cases. The first one is for defining the placement of a volume
with respect to its container reference frame. This frame will be called
'master' and the frame of the positioned volume - 'local'. If T is a
transformation used for positioning volume daughters, then:

~~~ {.cpp}
         MASTER = T * LOCAL
~~~

  Therefore a local-to-master conversion will be performed by using T, while
a master-to-local by using its inverse. The second use case is the computation
of the global transformation of a given object in the geometry. Since the
geometry is built as 'volumes-inside-volumes', this global transformation
represent the pile-up of all local transformations in the corresponding
branch. The conversion from the global reference frame and the given object
is also called master-to-local, but it is handled by the manager class.
  A general homogenous transformation is defined as a 4x4 matrix embedding
a rotation, a translation and a scale. The advantage of this description
is that each basic transformation can be represented as a homogenous matrix,
composition being performed as simple matrix multiplication.

  Rotation:                      Inverse rotation:

~~~ {.cpp}
        r11  r12  r13   0              r11  r21  r31   0
        r21  r22  r23   0              r12  r22  r32   0
        r31  r32  r33   0              r13  r23  r33   0
         0    0    0    1               0    0    0    1
~~~

  Translation:                   Inverse translation:

~~~ {.cpp}
         1    0    0    tx               1    0    0   -tx
         0    1    0    ty               0    1    0   -ty
         0    0    1    tz               0    0    1   -tz
         0    0    0    1                0    0    0   1
~~~

  Scale:                         Inverse scale:

~~~ {.cpp}
         sx   0    0    0              1/sx  0    0    0
         0    sy   0    0               0   1/sy  0    0
         0    0    sz   0               0    0   1/sz  0
         0    0    0    1               0    0    0    1
~~~

 where:
       - `rij` are the 3x3 rotation matrix components,
       - `tx`, `ty`, `tz` are the translation components
       - `sx`, `sy`, `sz` are arbitrary scale constants on each axis,

  The disadvantage in using this approach is that computation for 4x4 matrices
is expensive. Even combining two translation would become a multiplication
of their corresponding matrices, which is quite an undesired effect. On the
other hand, it is not a good idea to store a translation as a block of 16
numbers. We have therefore chosen to implement each basic transformation type
as a class deriving from the same basic abstract class and handling its specific
data and point/vector transformation algorithms.

\image html geom_transf.jpg

### The base class TGeoMatrix defines abstract metods for:

#### translation, rotation and scale getters. Every derived class stores only
  its specific data, e.g. a translation stores an array of 3 doubles and a
  rotation an array of 9. However, asking which is the rotation array of a
  TGeoTranslation through the base TGeoMatrix interface is a legal operation.
  The answer in this case is a pointer to a global constant array representing
  an identity rotation.

~~~ {.cpp}
     Double_t *TGeoMatrix::GetTranslation()
     Double_t *TGeoMatrix::GetRotation()
     Double_t *TGeoMatrix::GetScale()
~~~

#### MasterToLocal() and LocalToMaster() point and vector transformations :

~~~ {.cpp}
     void      TGeoMatrix::MasterToLocal(const Double_t *master, Double_t *local)
     void      TGeoMatrix::LocalToMaster(const Double_t *local, Double_t *master)
     void      TGeoMatrix::MasterToLocalVect(const Double_t *master, Double_t *local)
     void      TGeoMatrix::LocalToMasterVect(const Double_t *local, Double_t *master)
~~~

  These allow correct conversion also for reflections.

#### Transformation type getters :

~~~ {.cpp}
     Bool_t    TGeoMatrix::IsIdentity()
     Bool_t    TGeoMatrix::IsTranslation()
     Bool_t    TGeoMatrix::IsRotation()
     Bool_t    TGeoMatrix::IsScale()
     Bool_t    TGeoMatrix::IsCombi() (translation + rotation)
     Bool_t    TGeoMatrix::IsGeneral() (translation + rotation + scale)
~~~

  Combinations of basic transformations are represented by specific classes
deriving from TGeoMatrix. In order to define a matrix as a combination of several
others, a special class TGeoHMatrix is provided. Here is an example of matrix
creation :

### Matrix creation example:

~~~ {.cpp}
  root[0] TGeoRotation r1,r2;
          r1.SetAngles(90,0,30);        // rotation defined by Euler angles
          r2.SetAngles(90,90,90,180,0,0); // rotation defined by GEANT3 angles
          TGeoTranslation t1(-10,10,0);
          TGeoTranslation t2(10,-10,5);
          TGeoCombiTrans c1(t1,r1);
          TGeoCombiTrans c2(t2,r2);
          TGeoHMatrix h = c1 * c2; // composition is done via TGeoHMatrix class
  root[7] TGeoHMatrix *ph = new TGeoHMatrix(hm); // this is the one we want to
                                               // use for positioning a volume
  root[8] ph->Print();
          ...
          pVolume->AddNode(pVolDaughter,id,ph) // now ph is owned by the manager
~~~

### Rule for matrix creation:
 Unless explicitly used for positioning nodes (TGeoVolume::AddNode()) all
matrices deletion have to be managed by users. Matrices passed to geometry
have to be created by using new() operator and their deletion is done by
TGeoManager class.

### Available geometrical transformations

#### TGeoTranslation
Represent a (dx,dy,dz) translation. Data members:
   Double_t fTranslation[3]. Translations can be added/subtracted.

~~~ {.cpp}
   TGeoTranslation t1;
   t1->SetTranslation(-5,10,4);
   TGeoTranslation *t2 = new TGeoTranslation(4,3,10);
   t2->Subtract(&t1);
~~~

#### Rotations
 Represent a pure rotation. Data members: Double_t fRotationMatrix[3*3].
   Rotations can be defined either by Euler angles, either, by GEANT3 angles :

~~~ {.cpp}
   TGeoRotation *r1 = new TGeoRotation();
   r1->SetAngles(phi, theta, psi); // all angles in degrees
~~~

   This represent the composition of : first a rotation about Z axis with
   angle phi, then a rotation with theta about the rotated X axis, and
   finally a rotation with psi about the new Z axis.

~~~ {.cpp}
   r1->SetAngles(th1,phi1, th2,phi2, th3,phi3)
~~~

   This is a rotation defined in GEANT3 style. Theta and phi are the spherical
   angles of each axis of the rotated coordinate system with respect to the
   initial one. This construction allows definition of malformed rotations,
   e.g. not orthogonal. A check is performed and an error message is issued
   in this case.

   Specific utilities : determinant, inverse.

#### Scale transformations
     Represent a scale shrinking/enlargement. Data
     members :Double_t fScale[3]. Not fully implemented yet.

#### Combined transformations
Represent a rotation followed by a translation.
Data members: Double_t fTranslation[3], TGeoRotation *fRotation.

~~~ {.cpp}
   TGeoRotation *rot = new TGeoRotation("rot",10,20,30);
   TGeoTranslation trans;
   ...
   TGeoCombiTrans *c1 = new TGeoCombiTrans(trans, rot);
   TGeoCombiTrans *c2 = new TGeoCombiTrans("somename",10,20,30,rot)
~~~


#### TGeoGenTrans
Combined transformations including a scale. Not implemented.

#### TGeoIdentity
A generic singleton matrix representing a identity transformation
   NOTE: identified by the global variable gGeoIdentity.
*/

#include "Riostream.h"
#include "TObjArray.h"

#include "TGeoManager.h"
#include "TGeoMatrix.h"
#include "TMath.h"

TGeoIdentity *gGeoIdentity = 0;
const Int_t kN3 = 3*sizeof(Double_t);
const Int_t kN9 = 9*sizeof(Double_t);

// statics and globals

ClassImp(TGeoMatrix);

////////////////////////////////////////////////////////////////////////////////
/// dummy constructor

TGeoMatrix::TGeoMatrix()
{
  ResetBit(kGeoMatrixBits);
}

////////////////////////////////////////////////////////////////////////////////
/// copy constructor

TGeoMatrix::TGeoMatrix(const TGeoMatrix &other)
           :TNamed(other)
{
   ResetBit(kGeoRegistered);
}

////////////////////////////////////////////////////////////////////////////////
/// Constructor

TGeoMatrix::TGeoMatrix(const char *name)
           :TNamed(name, "")
{
}

////////////////////////////////////////////////////////////////////////////////
/// Destructor

TGeoMatrix::~TGeoMatrix()
{
   if (IsRegistered() && gGeoManager) {
      if (!gGeoManager->IsCleaning()) {
         gGeoManager->GetListOfMatrices()->Remove(this);
         Warning("dtor", "Registered matrix %s was removed", GetName());
      }
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Returns true if no rotation or the rotation is about Z axis

Bool_t TGeoMatrix::IsRotAboutZ() const
{
   if (IsIdentity()) return kTRUE;
   const Double_t *rot = GetRotationMatrix();
   if (TMath::Abs(rot[6])>1E-9) return kFALSE;
   if (TMath::Abs(rot[7])>1E-9) return kFALSE;
   if ((1.-TMath::Abs(rot[8]))>1E-9) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Get total size in bytes of this

Int_t TGeoMatrix::GetByteCount() const
{
   Int_t count = 4+28+strlen(GetName())+strlen(GetTitle()); // fId + TNamed
   if (IsTranslation()) count += 12;
   if (IsScale()) count += 12;
   if (IsCombi() || IsGeneral()) count += 4 + 36;
   return count;
}

////////////////////////////////////////////////////////////////////////////////
/// Provide a pointer name containing uid.

char *TGeoMatrix::GetPointerName() const
{
   static TString name;
   name = TString::Format("pMatrix%d", GetUniqueID());
   return (char*)name.Data();
}

////////////////////////////////////////////////////////////////////////////////
/// The homogenous matrix associated with the transformation is used for
/// piling up's and visualization. A homogenous matrix is a 4*4 array
/// containing the translation, the rotation and the scale components
/// ~~~ {.cpp}
///          | R00*sx  R01    R02    dx |
///          | R10    R11*sy  R12    dy |
///          | R20     R21   R22*sz  dz |
///          |  0       0      0      1 |
/// ~~~
///   where Rij is the rotation matrix, (sx, sy, sz) is the scale
/// transformation and (dx, dy, dz) is the translation.

void TGeoMatrix::GetHomogenousMatrix(Double_t *hmat) const
{
   Double_t *hmatrix = hmat;
   const Double_t *mat = GetRotationMatrix();
   for (Int_t i=0; i<3; i++) {
      memcpy(hmatrix, mat, kN3);
      mat     += 3;
      hmatrix += 3;
      *hmatrix = 0.0;
      hmatrix++;
   }
   memcpy(hmatrix, GetTranslation(), kN3);
   hmatrix = hmat;
   if (IsScale()) {
      for (Int_t i=0; i<3; i++) {
         *hmatrix *= GetScale()[i];
         hmatrix  += 5;
      }
   }
   hmatrix[15] = 1.;
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoMatrix::LocalToMaster(const Double_t *local, Double_t *master) const
{
   if (IsIdentity()) {
      memcpy(master, local, kN3);
      return;
   }
   Int_t i;
   const Double_t *tr = GetTranslation();
   if (!IsRotation()) {
      for (i=0; i<3; i++) master[i] = tr[i] + local[i];
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   for (i=0; i<3; i++) {
      master[i] = tr[i]
                + local[0]*rot[3*i]
                + local[1]*rot[3*i+1]
                + local[2]*rot[3*i+2];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoMatrix::LocalToMasterVect(const Double_t *local, Double_t *master) const
{
   if (!IsRotation()) {
      memcpy(master, local, kN3);
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   for (Int_t i=0; i<3; i++) {
      master[i] = local[0]*rot[3*i]
                + local[1]*rot[3*i+1]
                + local[2]*rot[3*i+2];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoMatrix::LocalToMasterBomb(const Double_t *local, Double_t *master) const
{
   if (IsIdentity()) {
      memcpy(master, local, kN3);
      return;
   }
   Int_t i;
   const Double_t *tr = GetTranslation();
   Double_t bombtr[3] = {0.,0.,0.};
   gGeoManager->BombTranslation(tr, &bombtr[0]);
   if (!IsRotation()) {
      for (i=0; i<3; i++) master[i] = bombtr[i] + local[i];
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   for (i=0; i<3; i++) {
      master[i] = bombtr[i]
                + local[0]*rot[3*i]
                + local[1]*rot[3*i+1]
                + local[2]*rot[3*i+2];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoMatrix::MasterToLocal(const Double_t *master, Double_t *local) const
{
   if (IsIdentity()) {
      memcpy(local, master, kN3);
      return;
   }
   const Double_t *tr  = GetTranslation();
   Double_t mt0  = master[0]-tr[0];
   Double_t mt1  = master[1]-tr[1];
   Double_t mt2  = master[2]-tr[2];
   if (!IsRotation()) {
      local[0] = mt0;
      local[1] = mt1;
      local[2] = mt2;
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   local[0] = mt0*rot[0] + mt1*rot[3] + mt2*rot[6];
   local[1] = mt0*rot[1] + mt1*rot[4] + mt2*rot[7];
   local[2] = mt0*rot[2] + mt1*rot[5] + mt2*rot[8];
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoMatrix::MasterToLocalVect(const Double_t *master, Double_t *local) const
{
   if (!IsRotation()) {
      memcpy(local, master, kN3);
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   for (Int_t i=0; i<3; i++) {
      local[i] = master[0]*rot[i]
               + master[1]*rot[i+3]
               + master[2]*rot[i+6];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoMatrix::MasterToLocalBomb(const Double_t *master, Double_t *local) const
{
   if (IsIdentity()) {
      memcpy(local, master, kN3);
      return;
   }
   const Double_t *tr = GetTranslation();
   Double_t bombtr[3] = {0.,0.,0.};
   Int_t i;
   gGeoManager->UnbombTranslation(tr, &bombtr[0]);
   if (!IsRotation()) {
      for (i=0; i<3; i++) local[i] = master[i]-bombtr[i];
      return;
   }
   const Double_t *rot = GetRotationMatrix();
   for (i=0; i<3; i++) {
      local[i] = (master[0]-bombtr[0])*rot[i]
               + (master[1]-bombtr[1])*rot[i+3]
               + (master[2]-bombtr[2])*rot[i+6];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Normalize a vector.

void TGeoMatrix::Normalize(Double_t *vect)
{
   Double_t normfactor = vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2];
   if (normfactor <= 1E-10) return;
   normfactor = 1./TMath::Sqrt(normfactor);
   vect[0] *= normfactor;
   vect[1] *= normfactor;
   vect[2] *= normfactor;
}

////////////////////////////////////////////////////////////////////////////////
/// print the matrix in 4x4 format

void TGeoMatrix::Print(Option_t *) const
{
   const Double_t *rot = GetRotationMatrix();
   const Double_t *tr  = GetTranslation();
   printf("matrix %s - tr=%d  rot=%d  refl=%d  scl=%d shr=%d reg=%d own=%d\n", GetName(),(Int_t)IsTranslation(),
          (Int_t)IsRotation(), (Int_t)IsReflection(), (Int_t)IsScale(), (Int_t)IsShared(), (Int_t)IsRegistered(),
          (Int_t)IsOwned());
   printf("%10.6f%12.6f%12.6f    Tx = %10.6f\n", rot[0], rot[1], rot[2], tr[0]);
   printf("%10.6f%12.6f%12.6f    Ty = %10.6f\n", rot[3], rot[4], rot[5], tr[1]);
   printf("%10.6f%12.6f%12.6f    Tz = %10.6f\n", rot[6], rot[7], rot[8], tr[2]);
   if (IsScale()) {
      const Double_t *scl  = GetScale();
      printf("Sx=%10.6fSy=%12.6fSz=%12.6f\n", scl[0], scl[1], scl[2]);
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to YZ.

void TGeoMatrix::ReflectX(Bool_t, Bool_t)
{
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to ZX.

void TGeoMatrix::ReflectY(Bool_t, Bool_t)
{
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to XY.

void TGeoMatrix::ReflectZ(Bool_t, Bool_t)
{
}

////////////////////////////////////////////////////////////////////////////////
/// Register the matrix in the current manager, which will become the owner.

void TGeoMatrix::RegisterYourself()
{
   if (!gGeoManager) {
      Warning("RegisterYourself", "cannot register without geometry");
      return;
   }
   if (!IsRegistered()) {
      gGeoManager->RegisterMatrix(this);
      SetBit(kGeoRegistered);
   }
}

////////////////////////////////////////////////////////////////////////////////
/// If no name was supplied in the ctor, the type of transformation is checked.
/// A letter will be prepended to the name :
///  - t - translation
///  - r - rotation
///  - s - scale
///  - c - combi (translation + rotation)
///  - g - general (tr+rot+scale)
/// The index of the transformation in gGeoManager list of transformations will
/// be appended.

void TGeoMatrix::SetDefaultName()
{
   if (!gGeoManager) return;
   if (strlen(GetName())) return;
   char type = 'n';
   if (IsTranslation()) type = 't';
   if (IsRotation()) type = 'r';
   if (IsScale()) type = 's';
   if (IsCombi()) type = 'c';
   if (IsGeneral()) type = 'g';
   TObjArray *matrices = gGeoManager->GetListOfMatrices();
   Int_t index = 0;
   if (matrices) index =matrices->GetEntriesFast() - 1;
   TString name = TString::Format("%c%d", type, index);
   SetName(name);
}

/** \class TGeoTranslation
\ingroup Geometry_classes

Class describing translations. A translation is
basically an array of 3 doubles matching the positions 12, 13
and 14 in the homogenous matrix description.
*/

ClassImp(TGeoTranslation);

////////////////////////////////////////////////////////////////////////////////
/// Default constructor

TGeoTranslation::TGeoTranslation()
{
   for (Int_t i=0; i<3; i++) fTranslation[i] = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// Copy ctor.

TGeoTranslation::TGeoTranslation(const TGeoTranslation &other)
                :TGeoMatrix(other)
{
   SetTranslation(other);
}

////////////////////////////////////////////////////////////////////////////////
/// Ctor. based on a general matrix

TGeoTranslation::TGeoTranslation(const TGeoMatrix &other)
                :TGeoMatrix(other)
{
   ResetBit(kGeoRotation);
   ResetBit(kGeoScale);
   SetTranslation(other);
}

////////////////////////////////////////////////////////////////////////////////
/// Default constructor defining the translation

TGeoTranslation::TGeoTranslation(Double_t dx, Double_t dy, Double_t dz)
                :TGeoMatrix("")
{
   if (dx || dy || dz) SetBit(kGeoTranslation);
   SetTranslation(dx, dy, dz);
}

////////////////////////////////////////////////////////////////////////////////
/// Default constructor defining the translation

TGeoTranslation::TGeoTranslation(const char *name, Double_t dx, Double_t dy, Double_t dz)
                :TGeoMatrix(name)
{
   if (dx || dy || dz) SetBit(kGeoTranslation);
   SetTranslation(dx, dy, dz);
}

////////////////////////////////////////////////////////////////////////////////
/// Assignment from a general matrix

TGeoTranslation& TGeoTranslation::operator = (const TGeoMatrix &matrix)
{
   if (&matrix == this) return *this;
   Bool_t registered = TestBit(kGeoRegistered);
   TNamed::operator=(matrix);
   SetTranslation(matrix);
   SetBit(kGeoRegistered,registered);
   ResetBit(kGeoRotation);
   ResetBit(kGeoScale);
   return *this;
}

////////////////////////////////////////////////////////////////////////////////
/// Translation composition

TGeoTranslation &TGeoTranslation::operator*=(const TGeoTranslation &right)
{
   const Double_t *tr = right.GetTranslation();
   fTranslation[0] += tr[0];
   fTranslation[1] += tr[1];
   fTranslation[2] += tr[2];
   if (!IsTranslation()) SetBit(kGeoTranslation, right.IsTranslation());
   return *this;
}

TGeoTranslation TGeoTranslation::operator*(const TGeoTranslation &right) const
{
   TGeoTranslation t = *this;
   t *= right;
   return t;
}

TGeoHMatrix TGeoTranslation::operator*(const TGeoMatrix &right) const
{
   TGeoHMatrix t = *this;
   t *= right;
   return t;
}

////////////////////////////////////////////////////////////////////////////////
/// Is-equal operator

Bool_t TGeoTranslation::operator ==(const TGeoTranslation &other) const
{
   if (&other == this) return kTRUE;
   const Double_t *tr = GetTranslation();
   const Double_t *otr = other.GetTranslation();
   for (auto i=0; i<3; i++)
      if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoTranslation::Inverse() const
{
   TGeoHMatrix h;
   h = *this;
   Double_t tr[3];
   tr[0] = -fTranslation[0];
   tr[1] = -fTranslation[1];
   tr[2] = -fTranslation[2];
   h.SetTranslation(tr);
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// Adding a translation to this one

void TGeoTranslation::Add(const TGeoTranslation *other)
{
   const Double_t *trans = other->GetTranslation();
   for (Int_t i=0; i<3; i++)
      fTranslation[i] += trans[i];
}

////////////////////////////////////////////////////////////////////////////////
/// Make a clone of this matrix.

TGeoMatrix *TGeoTranslation::MakeClone() const
{
   TGeoMatrix *matrix = new TGeoTranslation(*this);
   return matrix;
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about X axis of the master frame with angle expressed in degrees.

void TGeoTranslation::RotateX(Double_t /*angle*/)
{
   Warning("RotateX", "Not implemented. Use TGeoCombiTrans instead");
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Y axis of the master frame with angle expressed in degrees.

void TGeoTranslation::RotateY(Double_t /*angle*/)
{
   Warning("RotateY", "Not implemented. Use TGeoCombiTrans instead");
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Z axis of the master frame with angle expressed in degrees.

void TGeoTranslation::RotateZ(Double_t /*angle*/)
{
   Warning("RotateZ", "Not implemented. Use TGeoCombiTrans instead");
}

////////////////////////////////////////////////////////////////////////////////
/// Subtracting a translation from this one

void TGeoTranslation::Subtract(const TGeoTranslation *other)
{
   const Double_t *trans = other->GetTranslation();
   for (Int_t i=0; i<3; i++)
      fTranslation[i] -= trans[i];
}

////////////////////////////////////////////////////////////////////////////////
/// Set translation components

void TGeoTranslation::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
{
   fTranslation[0] = dx;
   fTranslation[1] = dy;
   fTranslation[2] = dz;
   if (dx || dy || dz) SetBit(kGeoTranslation);
   else                ResetBit(kGeoTranslation);
}

////////////////////////////////////////////////////////////////////////////////
/// Set translation components

void TGeoTranslation::SetTranslation(const TGeoMatrix &other)
{
   SetBit(kGeoTranslation, other.IsTranslation());
   const Double_t *transl = other.GetTranslation();
   memcpy(fTranslation, transl, kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoTranslation::LocalToMaster(const Double_t *local, Double_t *master) const
{
   const Double_t *tr = GetTranslation();
   for (Int_t i=0; i<3; i++)
      master[i] = tr[i] + local[i];
}

////////////////////////////////////////////////////////////////////////////////
/// convert a vector to MARS

void TGeoTranslation::LocalToMasterVect(const Double_t *local, Double_t *master) const
{
   memcpy(master, local, kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoTranslation::LocalToMasterBomb(const Double_t *local, Double_t *master) const
{
   const Double_t *tr = GetTranslation();
   Double_t bombtr[3] = {0.,0.,0.};
   gGeoManager->BombTranslation(tr, &bombtr[0]);
   for (Int_t i=0; i<3; i++)
      master[i] = bombtr[i] + local[i];
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoTranslation::MasterToLocal(const Double_t *master, Double_t *local) const
{
   const Double_t *tr = GetTranslation();
   for (Int_t i=0; i<3; i++)
      local[i] =  master[i]-tr[i];
}

////////////////////////////////////////////////////////////////////////////////
/// convert a vector from MARS to local

void TGeoTranslation::MasterToLocalVect(const Double_t *master, Double_t *local) const
{
   memcpy(local, master, kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoTranslation::MasterToLocalBomb(const Double_t *master, Double_t *local) const
{
   const Double_t *tr = GetTranslation();
   Double_t bombtr[3] = {0.,0.,0.};
   gGeoManager->UnbombTranslation(tr, &bombtr[0]);
   for (Int_t i=0; i<3; i++)
      local[i] =  master[i]-bombtr[i];
}

////////////////////////////////////////////////////////////////////////////////
/// Save a primitive as a C++ statement(s) on output stream "out".

void TGeoTranslation::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
{
   if (TestBit(kGeoSavePrimitive)) return;
   out << "   // Translation: " << GetName() << std::endl;
   out << "   dx = " << fTranslation[0] << ";" << std::endl;
   out << "   dy = " << fTranslation[1] << ";" << std::endl;
   out << "   dz = " << fTranslation[2] << ";" << std::endl;
   out << "   TGeoTranslation *" << GetPointerName() << " = new TGeoTranslation(\"" << GetName() << "\",dx,dy,dz);" << std::endl;
   TObject::SetBit(kGeoSavePrimitive);
}

/** \class TGeoRotation
\ingroup Geometry_classes
Class describing rotations. A rotation is a 3*3 array
Column vectors has to be orthogonal unit vectors.
*/

ClassImp(TGeoRotation);

////////////////////////////////////////////////////////////////////////////////
/// Default constructor.

TGeoRotation::TGeoRotation()
{
   for (Int_t i=0; i<9; i++) {
      if (i%4) fRotationMatrix[i] = 0;
      else fRotationMatrix[i] = 1.0;
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Copy ctor.

TGeoRotation::TGeoRotation(const TGeoRotation &other)
             :TGeoMatrix(other)
{
   SetRotation(other);
}

////////////////////////////////////////////////////////////////////////////////
/// Copy ctor.

TGeoRotation::TGeoRotation(const TGeoMatrix &other)
             :TGeoMatrix(other)
{
   ResetBit(kGeoTranslation);
   ResetBit(kGeoScale);
   SetRotation(other);
}

////////////////////////////////////////////////////////////////////////////////
/// Named rotation constructor

TGeoRotation::TGeoRotation(const char *name)
             :TGeoMatrix(name)
{
   for (Int_t i=0; i<9; i++) {
      if (i%4) fRotationMatrix[i] = 0;
      else fRotationMatrix[i] = 1.0;
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Default rotation constructor with Euler angles. Phi is the rotation angle about
/// Z axis  and is done first, theta is the rotation about new X and is done
/// second, psi is the rotation angle about new Z and is done third. All angles are in
/// degrees.

TGeoRotation::TGeoRotation(const char *name, Double_t phi, Double_t theta, Double_t psi)
             :TGeoMatrix(name)
{
   SetAngles(phi, theta, psi);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotation constructor a la GEANT3. Angles theta(i), phi(i) are the polar and azimuthal
/// angles of the (i) axis of the rotated system with respect to the initial non-rotated
/// system.
///   Example : the identity matrix (no rotation) is composed by
///      theta1=90, phi1=0, theta2=90, phi2=90, theta3=0, phi3=0
///   SetBit(kGeoRotation);

TGeoRotation::TGeoRotation(const char *name, Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
                           Double_t theta3, Double_t phi3)
             :TGeoMatrix(name)
{
   SetAngles(theta1, phi1, theta2, phi2, theta3, phi3);
}

////////////////////////////////////////////////////////////////////////////////
/// Assignment from a general matrix

TGeoRotation& TGeoRotation::operator = (const TGeoMatrix &other)
{
   if (&other == this) return *this;
   Bool_t registered = TestBit(kGeoRegistered);
   TNamed::operator=(other);
   SetRotation(other);
   SetBit(kGeoRegistered,registered);
   ResetBit(kGeoTranslation);
   ResetBit(kGeoScale);
   return *this;
}

////////////////////////////////////////////////////////////////////////////////
/// Composition

TGeoRotation &TGeoRotation::operator*=(const TGeoRotation &right)
{
   if (!right.IsRotation()) return *this;
   MultiplyBy(&right, true);
   return *this;
}

TGeoRotation TGeoRotation::operator*(const TGeoRotation &right) const
{
   TGeoRotation r = *this;
   r *= right;
   return r;
}

TGeoHMatrix TGeoRotation::operator*(const TGeoMatrix &right) const
{
   TGeoHMatrix t = *this;
   t *= right;
   return t;
}

////////////////////////////////////////////////////////////////////////////////
/// Is-equal operator

Bool_t TGeoRotation::operator ==(const TGeoRotation &other) const
{
   if (&other == this) return kTRUE;
   const Double_t *rot = GetRotationMatrix();
   const Double_t *orot = other.GetRotationMatrix();
   for (auto i=0; i<9; i++)
      if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoRotation::Inverse() const
{
   TGeoHMatrix h;
   h = *this;
   Double_t newrot[9];
   newrot[0] = fRotationMatrix[0];
   newrot[1] = fRotationMatrix[3];
   newrot[2] = fRotationMatrix[6];
   newrot[3] = fRotationMatrix[1];
   newrot[4] = fRotationMatrix[4];
   newrot[5] = fRotationMatrix[7];
   newrot[6] = fRotationMatrix[2];
   newrot[7] = fRotationMatrix[5];
   newrot[8] = fRotationMatrix[8];
   h.SetRotation(newrot);
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// Perform orthogonality test for rotation.

Bool_t TGeoRotation::IsValid() const
{
   const Double_t *r = fRotationMatrix;
   Double_t cij;
   for (Int_t i=0; i<2; i++) {
      for (Int_t j=i+1; j<3; j++) {
         // check columns
         cij = TMath::Abs(r[i]*r[j]+r[i+3]*r[j+3]+r[i+6]*r[j+6]);
         if (cij>1E-4) return kFALSE;
         // check rows
         cij = TMath::Abs(r[3*i]*r[3*j]+r[3*i+1]*r[3*j+1]+r[3*i+2]*r[3*j+2]);
         if (cij>1E-4) return kFALSE;
      }
   }
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// reset data members

void TGeoRotation::Clear(Option_t *)
{
   memcpy(fRotationMatrix,kIdentityMatrix,kN9);
   ResetBit(kGeoRotation);
}

////////////////////////////////////////////////////////////////////////////////
/// Perform a rotation about Z having the sine/cosine of the rotation angle.

void TGeoRotation::FastRotZ(const Double_t *sincos)
{
   fRotationMatrix[0] = sincos[1];
   fRotationMatrix[1] = -sincos[0];
   fRotationMatrix[3] = sincos[0];
   fRotationMatrix[4] = sincos[1];
   SetBit(kGeoRotation);
}

////////////////////////////////////////////////////////////////////////////////
/// Returns rotation angle about Z axis in degrees. If the rotation is a pure
/// rotation about Z, fixX parameter does not matter, otherwise its meaning is:
///    - fixX = true  : result is the phi angle of the projection of the rotated X axis in the un-rotated XY
///    - fixX = false : result is the phi angle of the projection of the rotated Y axis - 90 degrees

Double_t TGeoRotation::GetPhiRotation(Bool_t fixX) const
{
   Double_t phi;
   if (fixX) phi = 180.*TMath::ATan2(-fRotationMatrix[1],fRotationMatrix[4])/TMath::Pi();
   else      phi = 180.*TMath::ATan2(fRotationMatrix[3], fRotationMatrix[0])/TMath::Pi();
   return phi;
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse

void TGeoRotation::LocalToMaster(const Double_t *local, Double_t *master) const
{
   const Double_t *rot = GetRotationMatrix();
   for (Int_t i=0; i<3; i++) {
      master[i] = local[0]*rot[3*i]
                + local[1]*rot[3*i+1]
                + local[2]*rot[3*i+2];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// convert a point by multiplying its column vector (x, y, z, 1) to matrix

void TGeoRotation::MasterToLocal(const Double_t *master, Double_t *local) const
{
   const Double_t *rot = GetRotationMatrix();
   for (Int_t i=0; i<3; i++) {
      local[i] = master[0]*rot[i]
               + master[1]*rot[i+3]
               + master[2]*rot[i+6];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Make a clone of this matrix.

TGeoMatrix *TGeoRotation::MakeClone() const
{
   TGeoMatrix *matrix = new TGeoRotation(*this);
   return matrix;
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about X axis of the master frame with angle expressed in degrees.

void TGeoRotation::RotateX(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = fRotationMatrix[0];
   v[1] = fRotationMatrix[1];
   v[2] = fRotationMatrix[2];
   v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
   v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
   v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
   v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
   v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
   v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];

   memcpy(fRotationMatrix, v, kN9);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Y axis of the master frame with angle expressed in degrees.

void TGeoRotation::RotateY(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
   v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
   v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
   v[3] = fRotationMatrix[3];
   v[4] = fRotationMatrix[4];
   v[5] = fRotationMatrix[5];
   v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
   v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
   v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];

   memcpy(fRotationMatrix, v, kN9);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Z axis of the master frame with angle expressed in degrees.

void TGeoRotation::RotateZ(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
   v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
   v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
   v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
   v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
   v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
   v[6] = fRotationMatrix[6];
   v[7] = fRotationMatrix[7];
   v[8] = fRotationMatrix[8];

   memcpy(&fRotationMatrix[0],v,kN9);
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to YZ.

void TGeoRotation::ReflectX(Bool_t leftside, Bool_t)
{
   if (leftside) {
      fRotationMatrix[0]=-fRotationMatrix[0];
      fRotationMatrix[1]=-fRotationMatrix[1];
      fRotationMatrix[2]=-fRotationMatrix[2];
   } else {
      fRotationMatrix[0]=-fRotationMatrix[0];
      fRotationMatrix[3]=-fRotationMatrix[3];
      fRotationMatrix[6]=-fRotationMatrix[6];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to ZX.

void TGeoRotation::ReflectY(Bool_t leftside, Bool_t)
{
   if (leftside) {
      fRotationMatrix[3]=-fRotationMatrix[3];
      fRotationMatrix[4]=-fRotationMatrix[4];
      fRotationMatrix[5]=-fRotationMatrix[5];
   } else {
      fRotationMatrix[1]=-fRotationMatrix[1];
      fRotationMatrix[4]=-fRotationMatrix[4];
      fRotationMatrix[7]=-fRotationMatrix[7];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to XY.

void TGeoRotation::ReflectZ(Bool_t leftside, Bool_t)
{
   if (leftside) {
      fRotationMatrix[6]=-fRotationMatrix[6];
      fRotationMatrix[7]=-fRotationMatrix[7];
      fRotationMatrix[8]=-fRotationMatrix[8];
   } else {
      fRotationMatrix[2]=-fRotationMatrix[2];
      fRotationMatrix[5]=-fRotationMatrix[5];
      fRotationMatrix[8]=-fRotationMatrix[8];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Save a primitive as a C++ statement(s) on output stream "out".

void TGeoRotation::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
{
   if (TestBit(kGeoSavePrimitive)) return;
   out << "   // Rotation: " << GetName() << std::endl;
   Double_t th1,ph1,th2,ph2,th3,ph3;
   GetAngles(th1,ph1,th2,ph2,th3,ph3);
   out << "   thx = " << th1 << ";    phx = " << ph1 << ";" << std::endl;
   out << "   thy = " << th2 << ";    phy = " << ph2 << ";" << std::endl;
   out << "   thz = " << th3 << ";    phz = " << ph3 << ";" << std::endl;
   out << "   TGeoRotation *" << GetPointerName() << " = new TGeoRotation(\"" << GetName() << "\",thx,phx,thy,phy,thz,phz);" << std::endl;
   TObject::SetBit(kGeoSavePrimitive);
}

////////////////////////////////////////////////////////////////////////////////
/// Copy rotation elements from other rotation matrix.

void TGeoRotation::SetRotation(const TGeoMatrix &other)
{
   SetBit(kGeoRotation, other.IsRotation());
   SetMatrix(other.GetRotationMatrix());
}

////////////////////////////////////////////////////////////////////////////////
/// Set matrix elements according to Euler angles. Phi is the rotation angle about
/// Z axis  and is done first, theta is the rotation about new X and is done
/// second, psi is the rotation angle about new Z and is done third. All angles are in
/// degrees.

void TGeoRotation::SetAngles(Double_t phi, Double_t theta, Double_t psi)
{
   Double_t degrad = TMath::Pi()/180.;
   Double_t sinphi = TMath::Sin(degrad*phi);
   Double_t cosphi = TMath::Cos(degrad*phi);
   Double_t sinthe = TMath::Sin(degrad*theta);
   Double_t costhe = TMath::Cos(degrad*theta);
   Double_t sinpsi = TMath::Sin(degrad*psi);
   Double_t cospsi = TMath::Cos(degrad*psi);

   fRotationMatrix[0] =  cospsi*cosphi - costhe*sinphi*sinpsi;
   fRotationMatrix[1] = -sinpsi*cosphi - costhe*sinphi*cospsi;
   fRotationMatrix[2] =  sinthe*sinphi;
   fRotationMatrix[3] =  cospsi*sinphi + costhe*cosphi*sinpsi;
   fRotationMatrix[4] = -sinpsi*sinphi + costhe*cosphi*cospsi;
   fRotationMatrix[5] = -sinthe*cosphi;
   fRotationMatrix[6] =  sinpsi*sinthe;
   fRotationMatrix[7] =  cospsi*sinthe;
   fRotationMatrix[8] =  costhe;

   if (!IsValid()) Error("SetAngles", "invalid rotation (Euler angles : phi=%f theta=%f psi=%f)",phi,theta,psi);
   CheckMatrix();
}

////////////////////////////////////////////////////////////////////////////////
/// Set matrix elements in the GEANT3 way

void TGeoRotation::SetAngles(Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
                             Double_t theta3, Double_t phi3)
{
   Double_t degrad = TMath::Pi()/180.;
   fRotationMatrix[0] = TMath::Cos(degrad*phi1)*TMath::Sin(degrad*theta1);
   fRotationMatrix[3] = TMath::Sin(degrad*phi1)*TMath::Sin(degrad*theta1);
   fRotationMatrix[6] = TMath::Cos(degrad*theta1);
   fRotationMatrix[1] = TMath::Cos(degrad*phi2)*TMath::Sin(degrad*theta2);
   fRotationMatrix[4] = TMath::Sin(degrad*phi2)*TMath::Sin(degrad*theta2);
   fRotationMatrix[7] = TMath::Cos(degrad*theta2);
   fRotationMatrix[2] = TMath::Cos(degrad*phi3)*TMath::Sin(degrad*theta3);
   fRotationMatrix[5] = TMath::Sin(degrad*phi3)*TMath::Sin(degrad*theta3);
   fRotationMatrix[8] = TMath::Cos(degrad*theta3);
   // do the trick to eliminate most of the floating point errors
   for (Int_t i=0; i<9; i++) {
      if (TMath::Abs(fRotationMatrix[i])<1E-15) fRotationMatrix[i] = 0;
      if (TMath::Abs(fRotationMatrix[i]-1)<1E-15) fRotationMatrix[i] = 1;
      if (TMath::Abs(fRotationMatrix[i]+1)<1E-15) fRotationMatrix[i] = -1;
   }
   if (!IsValid()) Error("SetAngles", "invalid rotation (G3 angles, th1=%f phi1=%f, th2=%f ph2=%f, th3=%f phi3=%f)",
                          theta1,phi1,theta2,phi2,theta3,phi3);
   CheckMatrix();
}

////////////////////////////////////////////////////////////////////////////////
/// Retrieve rotation angles

void TGeoRotation::GetAngles(Double_t &theta1, Double_t &phi1, Double_t &theta2, Double_t &phi2,
                             Double_t &theta3, Double_t &phi3) const
{
   Double_t raddeg = 180./TMath::Pi();
   theta1 = raddeg*TMath::ACos(fRotationMatrix[6]);
   theta2 = raddeg*TMath::ACos(fRotationMatrix[7]);
   theta3 = raddeg*TMath::ACos(fRotationMatrix[8]);
   if (TMath::Abs(fRotationMatrix[0])<1E-6 && TMath::Abs(fRotationMatrix[3])<1E-6) phi1=0.;
   else phi1 = raddeg*TMath::ATan2(fRotationMatrix[3],fRotationMatrix[0]);
   if (phi1<0) phi1+=360.;
   if (TMath::Abs(fRotationMatrix[1])<1E-6 && TMath::Abs(fRotationMatrix[4])<1E-6) phi2=0.;
   else phi2 = raddeg*TMath::ATan2(fRotationMatrix[4],fRotationMatrix[1]);
   if (phi2<0) phi2+=360.;
   if (TMath::Abs(fRotationMatrix[2])<1E-6 && TMath::Abs(fRotationMatrix[5])<1E-6) phi3=0.;
   else phi3 = raddeg*TMath::ATan2(fRotationMatrix[5],fRotationMatrix[2]);
   if (phi3<0) phi3+=360.;
}

////////////////////////////////////////////////////////////////////////////////
/// Retrieve Euler angles.

void TGeoRotation::GetAngles(Double_t &phi, Double_t &theta, Double_t &psi) const
{
   const Double_t *m = fRotationMatrix;
   // Check if theta is 0 or 180.
   if (TMath::Abs(1.-TMath::Abs(m[8]))<1.e-9) {
      theta = TMath::ACos(m[8])*TMath::RadToDeg();
      phi = TMath::ATan2(-m[8]*m[1],m[0])*TMath::RadToDeg();
      psi = 0.; // convention, phi+psi matters
      return;
   }
   // sin(theta) != 0
   phi = TMath::ATan2(m[2],-m[5]);
   Double_t sphi = TMath::Sin(phi);
   if (TMath::Abs(sphi)<1.e-9) theta = -TMath::ASin(m[5]/TMath::Cos(phi))*TMath::RadToDeg();
   else theta = TMath::ASin(m[2]/sphi)*TMath::RadToDeg();
   phi *= TMath::RadToDeg();
   psi = TMath::ATan2(m[6],m[7])*TMath::RadToDeg();
}

////////////////////////////////////////////////////////////////////////////////
/// computes determinant of the rotation matrix

Double_t TGeoRotation::Determinant() const
{
   Double_t
   det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
         fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
         fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
         fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
         fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
         fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
   return det;
}

////////////////////////////////////////////////////////////////////////////////
/// performes an orthogonality check and finds if the matrix is a reflection
///   Warning("CheckMatrix", "orthogonality check not performed yet");

void TGeoRotation::CheckMatrix()
{
   if (Determinant() < 0) SetBit(kGeoReflection);
   Double_t dd = fRotationMatrix[0] + fRotationMatrix[4] + fRotationMatrix[8] - 3.;
   if (TMath::Abs(dd) < 1.E-12) ResetBit(kGeoRotation);
   else                         SetBit(kGeoRotation);
}

////////////////////////////////////////////////////////////////////////////////
/// Get the inverse rotation matrix (which is simply the transpose)

void TGeoRotation::GetInverse(Double_t *invmat) const
{
   if (!invmat) {
      Error("GetInverse", "no place to store the inverse matrix");
      return;
   }
   for (Int_t i=0; i<3; i++) {
      for (Int_t j=0; j<3; j++) {
         invmat[3*i+j] = fRotationMatrix[3*j+i];
      }
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply this rotation with the one specified by ROT.
/// -   after=TRUE (default): THIS*ROT
/// -   after=FALSE         : ROT*THIS

void TGeoRotation::MultiplyBy(const TGeoRotation *rot, Bool_t after)
{
   const Double_t *matleft, *matright;
   SetBit(kGeoRotation);
   Double_t  newmat[9] = {0};
   if (after) {
      matleft  = &fRotationMatrix[0];
      matright = rot->GetRotationMatrix();
   } else {
      matleft  = rot->GetRotationMatrix();
      matright = &fRotationMatrix[0];
   }
   for (Int_t i=0; i<3; i++) {
      for (Int_t j=0; j<3; j++) {
         for (Int_t k=0; k<3; k++) {
            newmat[3*i+j] += matleft[3*i+k] * matright[3*k+j];
         }
      }
   }
   memcpy(&fRotationMatrix[0], &newmat[0], kN9);
}

/** \class TGeoScale
\ingroup Geometry_classes
Class describing scale transformations. A scale is an
array of 3 doubles (sx, sy, sz) multiplying elements 0, 5 and 10
of the homogenous matrix. A scale is normalized : sx*sy*sz = 1
*/

ClassImp(TGeoScale);

////////////////////////////////////////////////////////////////////////////////
/// default constructor

TGeoScale::TGeoScale()
{
   SetBit(kGeoScale);
   for (Int_t i=0; i<3; i++) fScale[i] = 1.;
}

////////////////////////////////////////////////////////////////////////////////
/// Copy constructor

TGeoScale::TGeoScale(const TGeoScale &other)
          :TGeoMatrix(other)
{
   SetScale(other);
}

////////////////////////////////////////////////////////////////////////////////
/// Ctor. based on a general matrix

TGeoScale::TGeoScale(const TGeoMatrix &other)
                :TGeoMatrix(other)
{
   ResetBit(kGeoTranslation);
   ResetBit(kGeoRotation);
   SetScale(other);
}

////////////////////////////////////////////////////////////////////////////////
/// default constructor

TGeoScale::TGeoScale(Double_t sx, Double_t sy, Double_t sz)
          :TGeoMatrix("")
{
   SetBit(kGeoScale);
   SetScale(sx, sy, sz);
}

////////////////////////////////////////////////////////////////////////////////
/// default constructor

TGeoScale::TGeoScale(const char *name, Double_t sx, Double_t sy, Double_t sz)
          :TGeoMatrix(name)
{
   SetBit(kGeoScale);
   SetScale(sx, sy, sz);
}

////////////////////////////////////////////////////////////////////////////////
/// destructor

TGeoScale::~TGeoScale()
{
}

////////////////////////////////////////////////////////////////////////////////
/// Assignment from a general matrix

TGeoScale& TGeoScale::operator = (const TGeoMatrix &matrix)
{
   if (&matrix == this) return *this;
   Bool_t registered = TestBit(kGeoRegistered);
   TNamed::operator=(matrix);
   SetScale(matrix);
   SetBit(kGeoRegistered,registered);
   ResetBit(kGeoTranslation);
   ResetBit(kGeoRotation);
   return *this;
}

////////////////////////////////////////////////////////////////////////////////
/// Scale composition

TGeoScale &TGeoScale::operator*=(const TGeoScale &right)
{
   const Double_t *scl = right.GetScale();
   fScale[0] *= scl[0];
   fScale[1] *= scl[1];
   fScale[2] *= scl[2];
   SetBit(kGeoReflection, fScale[0] * fScale[1] * fScale[2] < 0);
   if (!IsScale()) SetBit(kGeoScale, right.IsScale());
   return *this;
}

TGeoScale TGeoScale::operator*(const TGeoScale &right) const
{
   TGeoScale s = *this;
   s *= right;
   return s;
}

TGeoHMatrix TGeoScale::operator*(const TGeoMatrix &right) const
{
   TGeoHMatrix t = *this;
   t *= right;
   return t;
}

////////////////////////////////////////////////////////////////////////////////
/// Is-equal operator

Bool_t TGeoScale::operator ==(const TGeoScale &other) const
{
   if (&other == this) return kTRUE;
   const Double_t *scl = GetScale();
   const Double_t *oscl = other.GetScale();
   for (auto i=0; i<3; i++)
      if (TMath::Abs(scl[i]-oscl[i])>1.E-10) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoScale::Inverse() const
{
   TGeoHMatrix h;
   h = *this;
   Double_t scale[3];
   scale[0] = 1./fScale[0];
   scale[1] = 1./fScale[1];
   scale[2] = 1./fScale[2];
   h.SetScale(scale);
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// scale setter

void TGeoScale::SetScale(Double_t sx, Double_t sy, Double_t sz)
{
   if (TMath::Abs(sx*sy*sz) < 1.E-10) {
      Error("SetScale", "Invalid scale %f, %f, %f for transformation %s",sx,sy,sx,GetName());
      return;
   }
   fScale[0] = sx;
   fScale[1] = sy;
   fScale[2] = sz;
   if (sx*sy*sz<0) SetBit(kGeoReflection);
   else            SetBit(kGeoReflection, kFALSE);
}

////////////////////////////////////////////////////////////////////////////////
/// Set scale from other transformation

void TGeoScale::SetScale(const TGeoMatrix &other)
{
   SetBit(kGeoScale, other.IsScale());
   SetBit(kGeoReflection, other.IsReflection());
   memcpy(fScale, other.GetScale(), kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Convert a local point to the master frame.

void TGeoScale::LocalToMaster(const Double_t *local, Double_t *master) const
{
   master[0] = local[0]*fScale[0];
   master[1] = local[1]*fScale[1];
   master[2] = local[2]*fScale[2];
}

////////////////////////////////////////////////////////////////////////////////
/// Convert the local distance along unit vector DIR to master frame. If DIR
/// is not specified perform a conversion such as the returned distance is the
/// the minimum for all possible directions.

Double_t TGeoScale::LocalToMaster(Double_t dist, const Double_t *dir) const
{
   Double_t scale;
   if (!dir) {
      scale = TMath::Abs(fScale[0]);
      if (TMath::Abs(fScale[1])<scale) scale = TMath::Abs(fScale[1]);
      if (TMath::Abs(fScale[2])<scale) scale = TMath::Abs(fScale[2]);
   } else {
      scale = fScale[0]*fScale[0]*dir[0]*dir[0] +
              fScale[1]*fScale[1]*dir[1]*dir[1] +
              fScale[2]*fScale[2]*dir[2]*dir[2];
      scale = TMath::Sqrt(scale);
   }
   return scale*dist;
}

////////////////////////////////////////////////////////////////////////////////
/// Make a clone of this matrix.

TGeoMatrix *TGeoScale::MakeClone() const
{
   TGeoMatrix *matrix = new TGeoScale(*this);
   return matrix;
}

////////////////////////////////////////////////////////////////////////////////
/// Convert a global point to local frame.

void TGeoScale::MasterToLocal(const Double_t *master, Double_t *local) const
{
   local[0] = master[0]/fScale[0];
   local[1] = master[1]/fScale[1];
   local[2] = master[2]/fScale[2];
}

////////////////////////////////////////////////////////////////////////////////
/// Convert the distance along unit vector DIR to local frame. If DIR
/// is not specified perform a conversion such as the returned distance is the
/// the minimum for all possible directions.

Double_t TGeoScale::MasterToLocal(Double_t dist, const Double_t *dir) const
{
   Double_t scale;
   if (!dir) {
      scale = TMath::Abs(fScale[0]);
      if (TMath::Abs(fScale[1])>scale) scale = TMath::Abs(fScale[1]);
      if (TMath::Abs(fScale[2])>scale) scale = TMath::Abs(fScale[2]);
      scale = 1./scale;
   } else {
      scale = (dir[0]*dir[0])/(fScale[0]*fScale[0]) +
              (dir[1]*dir[1])/(fScale[1]*fScale[1]) +
              (dir[2]*dir[2])/(fScale[2]*fScale[2]);
      scale = TMath::Sqrt(scale);
   }
   return scale*dist;
}

/** \class TGeoCombiTrans
\ingroup Geometry_classes
Class describing rotation + translation. Most frequently used in the description
of TGeoNode 's
*/

ClassImp(TGeoCombiTrans);

////////////////////////////////////////////////////////////////////////////////
/// dummy ctor

TGeoCombiTrans::TGeoCombiTrans()
{
   for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// Copy ctor from generic matrix.

TGeoCombiTrans::TGeoCombiTrans(const TGeoMatrix &other)
               :TGeoMatrix(other)
{
   ResetBit(kGeoScale);
   if (other.IsTranslation()) {
      SetBit(kGeoTranslation);
      memcpy(fTranslation,other.GetTranslation(),kN3);
   } else {
      for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   }
   if (other.IsRotation()) {
      SetBit(kGeoRotation);
      SetBit(kGeoMatrixOwned);
      fRotation = new TGeoRotation(other);
   }
   else fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// Constructor from a translation and a rotation.

TGeoCombiTrans::TGeoCombiTrans(const TGeoTranslation &tr, const TGeoRotation &rot)
{
   if (tr.IsTranslation()) {
      SetBit(kGeoTranslation);
      const Double_t *trans = tr.GetTranslation();
      memcpy(fTranslation, trans, kN3);
   } else {
      for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   }
   if (rot.IsRotation()) {
      SetBit(kGeoRotation);
      SetBit(kGeoMatrixOwned);
      fRotation = new TGeoRotation(rot);
      SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
   }
   else fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// Named ctor.

TGeoCombiTrans::TGeoCombiTrans(const char *name)
               :TGeoMatrix(name)
{
   for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// Constructor from a translation specified by X,Y,Z and a pointer to a rotation. The rotation will not be owned by this.

TGeoCombiTrans::TGeoCombiTrans(Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
               :TGeoMatrix("")
{
   SetTranslation(dx, dy, dz);
   fRotation = 0;
   SetRotation(rot);
}

////////////////////////////////////////////////////////////////////////////////
/// Named ctor

TGeoCombiTrans::TGeoCombiTrans(const char *name, Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
               :TGeoMatrix(name)
{
   SetTranslation(dx, dy, dz);
   fRotation = 0;
   SetRotation(rot);
}

////////////////////////////////////////////////////////////////////////////////
/// Assignment operator with generic matrix.

TGeoCombiTrans &TGeoCombiTrans::operator=(const TGeoMatrix &matrix)
{
   if (&matrix == this) return *this;
   Bool_t registered = TestBit(kGeoRegistered);
   Clear();
   TNamed::operator=(matrix);

   if (matrix.IsTranslation()) {
      memcpy(fTranslation,matrix.GetTranslation(),kN3);
   }
   if (matrix.IsRotation()) {
      if (!fRotation) {
         fRotation = new TGeoRotation();
         SetBit(kGeoMatrixOwned);
      } else {
         if (!TestBit(kGeoMatrixOwned)) {
            fRotation = new TGeoRotation();
            SetBit(kGeoMatrixOwned);
         }
      }
      fRotation->SetMatrix(matrix.GetRotationMatrix());
      fRotation->SetBit(kGeoReflection, matrix.TestBit(kGeoReflection));
      fRotation->SetBit(kGeoRotation);
   } else {
      if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
      ResetBit(kGeoMatrixOwned);
      fRotation = 0;
   }
   SetBit(kGeoRegistered,registered);
   ResetBit(kGeoScale);
   return *this;
}

////////////////////////////////////////////////////////////////////////////////
/// Is-equal operator

Bool_t TGeoCombiTrans::operator==(const TGeoMatrix &other) const
{
   if (&other == this) return kTRUE;
   const Double_t *tr = GetTranslation();
   const Double_t *otr = other.GetTranslation();
   for (auto i=0; i<3; i++) if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
   const Double_t *rot = GetRotationMatrix();
   const Double_t *orot = other.GetRotationMatrix();
   for (auto i=0; i<9; i++) if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Composition

TGeoCombiTrans &TGeoCombiTrans::operator*=(const TGeoMatrix &right)
{
   Multiply(&right);
   return *this;
}

TGeoCombiTrans TGeoCombiTrans::operator*(const TGeoMatrix &right) const
{
   TGeoHMatrix h = *this;
   h *= right;
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// destructor

TGeoCombiTrans::~TGeoCombiTrans()
{
   if (fRotation) {
      if(TestBit(TGeoMatrix::kGeoMatrixOwned) && !fRotation->IsRegistered()) delete fRotation;
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Reset translation/rotation to identity

void TGeoCombiTrans::Clear(Option_t *)
{
   if (IsTranslation()) {
      ResetBit(kGeoTranslation);
      memset(fTranslation, 0, kN3);
   }
   if (fRotation) {
      if (TestBit(kGeoMatrixOwned)) delete fRotation;
      fRotation = 0;
   }
   ResetBit(kGeoRotation);
   ResetBit(kGeoTranslation);
   ResetBit(kGeoMatrixOwned);
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoCombiTrans::Inverse() const
{
   TGeoHMatrix h;
   h = *this;
   Bool_t is_tr = IsTranslation();
   Bool_t is_rot = IsRotation();
   Double_t tr[3];
   Double_t newrot[9];
   const Double_t *rot = GetRotationMatrix();
   tr[0] = -fTranslation[0]*rot[0] - fTranslation[1]*rot[3] - fTranslation[2]*rot[6];
   tr[1] = -fTranslation[0]*rot[1] - fTranslation[1]*rot[4] - fTranslation[2]*rot[7];
   tr[2] = -fTranslation[0]*rot[2] - fTranslation[1]*rot[5] - fTranslation[2]*rot[8];
   h.SetTranslation(tr);
   newrot[0] = rot[0];
   newrot[1] = rot[3];
   newrot[2] = rot[6];
   newrot[3] = rot[1];
   newrot[4] = rot[4];
   newrot[5] = rot[7];
   newrot[6] = rot[2];
   newrot[7] = rot[5];
   newrot[8] = rot[8];
   h.SetRotation(newrot);
   h.SetBit(kGeoTranslation,is_tr);
   h.SetBit(kGeoRotation,is_rot);
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// Make a clone of this matrix.

TGeoMatrix *TGeoCombiTrans::MakeClone() const
{
   TGeoMatrix *matrix = new TGeoCombiTrans(*this);
   return matrix;
}

////////////////////////////////////////////////////////////////////////////////
/// multiply to the right with an other transformation
/// if right is identity matrix, just return

void TGeoCombiTrans::Multiply(const TGeoMatrix *right)
{
   if (right->IsIdentity()) return;
   TGeoHMatrix h = *this;
   h.Multiply(right);
   operator=(h);
}

////////////////////////////////////////////////////////////////////////////////
/// Register the matrix in the current manager, which will become the owner.

void TGeoCombiTrans::RegisterYourself()
{
   TGeoMatrix::RegisterYourself();
   if (fRotation && fRotation->IsRotation()) fRotation->RegisterYourself();
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about X axis with angle expressed in degrees.

void TGeoCombiTrans::RotateX(Double_t angle)
{
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   const Double_t *rot = fRotation->GetRotationMatrix();
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = rot[0];
   v[1] = rot[1];
   v[2] = rot[2];
   v[3] = c*rot[3]-s*rot[6];
   v[4] = c*rot[4]-s*rot[7];
   v[5] = c*rot[5]-s*rot[8];
   v[6] = s*rot[3]+c*rot[6];
   v[7] = s*rot[4]+c*rot[7];
   v[8] = s*rot[5]+c*rot[8];
   fRotation->SetMatrix(v);
   fRotation->SetBit(kGeoRotation);
   if (!IsTranslation()) return;
   v[0] = fTranslation[0];
   v[1] = c*fTranslation[1]-s*fTranslation[2];
   v[2] = s*fTranslation[1]+c*fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Y axis with angle expressed in degrees.

void TGeoCombiTrans::RotateY(Double_t angle)
{
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   const Double_t *rot = fRotation->GetRotationMatrix();
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*rot[0]+s*rot[6];
   v[1] = c*rot[1]+s*rot[7];
   v[2] = c*rot[2]+s*rot[8];
   v[3] = rot[3];
   v[4] = rot[4];
   v[5] = rot[5];
   v[6] = -s*rot[0]+c*rot[6];
   v[7] = -s*rot[1]+c*rot[7];
   v[8] = -s*rot[2]+c*rot[8];
   fRotation->SetMatrix(v);
   fRotation->SetBit(kGeoRotation);
   if (!IsTranslation()) return;
   v[0] = c*fTranslation[0]+s*fTranslation[2];
   v[1] = fTranslation[1];
   v[2] = -s*fTranslation[0]+c*fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Z axis with angle expressed in degrees.

void TGeoCombiTrans::RotateZ(Double_t angle)
{
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   const Double_t *rot = fRotation->GetRotationMatrix();
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*rot[0]-s*rot[3];
   v[1] = c*rot[1]-s*rot[4];
   v[2] = c*rot[2]-s*rot[5];
   v[3] = s*rot[0]+c*rot[3];
   v[4] = s*rot[1]+c*rot[4];
   v[5] = s*rot[2]+c*rot[5];
   v[6] = rot[6];
   v[7] = rot[7];
   v[8] = rot[8];
   fRotation->SetMatrix(v);
   fRotation->SetBit(kGeoRotation);
   if (!IsTranslation()) return;
   v[0] = c*fTranslation[0]-s*fTranslation[1];
   v[1] = s*fTranslation[0]+c*fTranslation[1];
   v[2] = fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to YZ.

void TGeoCombiTrans::ReflectX(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   fRotation->ReflectX(leftside);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to ZX.

void TGeoCombiTrans::ReflectY(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   fRotation->ReflectY(leftside);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to XY.

void TGeoCombiTrans::ReflectZ(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
   if (!fRotation || !TestBit(kGeoMatrixOwned)) {
      if (fRotation) fRotation = new TGeoRotation(*fRotation);
      else fRotation = new TGeoRotation();
      SetBit(kGeoMatrixOwned);
   }
   SetBit(kGeoRotation);
   fRotation->ReflectZ(leftside);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Save a primitive as a C++ statement(s) on output stream "out".

void TGeoCombiTrans::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
{
   if (TestBit(kGeoSavePrimitive)) return;
   out << "   // Combi transformation: " << GetName() << std::endl;
   out << "   dx = " << fTranslation[0] << ";" << std::endl;
   out << "   dy = " << fTranslation[1] << ";" << std::endl;
   out << "   dz = " << fTranslation[2] << ";" << std::endl;
   if (fRotation && fRotation->IsRotation()) {
      fRotation->SavePrimitive(out,option);
      out << "   " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\", dx,dy,dz,";
      out << fRotation->GetPointerName() << ");" << std::endl;
   } else {
      out << "   " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\");" << std::endl;
      out << "   " << GetPointerName() << "->SetTranslation(dx,dy,dz);" << std::endl;
   }
   TObject::SetBit(kGeoSavePrimitive);
}

////////////////////////////////////////////////////////////////////////////////
/// Assign a foreign rotation to the combi. The rotation is NOT owned by this.

void TGeoCombiTrans::SetRotation(const TGeoRotation *rot)
{
   if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
   fRotation = 0;
   ResetBit(TGeoMatrix::kGeoMatrixOwned);
   ResetBit(kGeoRotation);
   ResetBit(kGeoReflection);
   if (!rot) return;
   if (!rot->IsRotation()) return;

   SetBit(kGeoRotation);
   SetBit(kGeoReflection, rot->TestBit(kGeoReflection));
   TGeoRotation *rr = (TGeoRotation*)rot;
   fRotation = rr;
}

////////////////////////////////////////////////////////////////////////////////
/// Copy the rotation from another one.

void TGeoCombiTrans::SetRotation(const TGeoRotation &rot)
{
   if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
   fRotation = 0;
   if (!rot.IsRotation()) {
      ResetBit(kGeoRotation);
      ResetBit(kGeoReflection);
      ResetBit(TGeoMatrix::kGeoMatrixOwned);
      return;
   }

   SetBit(kGeoRotation);
   SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
   fRotation = new TGeoRotation(rot);
   SetBit(kGeoMatrixOwned);
}

////////////////////////////////////////////////////////////////////////////////
/// copy the translation component

void TGeoCombiTrans::SetTranslation(const TGeoTranslation &tr)
{
   if (tr.IsTranslation()) {
      SetBit(kGeoTranslation);
      const Double_t *trans = tr.GetTranslation();
      memcpy(fTranslation, trans, kN3);
   } else {
      if (!IsTranslation()) return;
      memset(fTranslation, 0, kN3);
      ResetBit(kGeoTranslation);
   }
}

////////////////////////////////////////////////////////////////////////////////
/// set the translation component

void TGeoCombiTrans::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
{
   fTranslation[0] = dx;
   fTranslation[1] = dy;
   fTranslation[2] = dz;
   if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
   else ResetBit(kGeoTranslation);
}

////////////////////////////////////////////////////////////////////////////////
/// set the translation component

void TGeoCombiTrans::SetTranslation(Double_t *vect)
{
   fTranslation[0] = vect[0];
   fTranslation[1] = vect[1];
   fTranslation[2] = vect[2];
   if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
   else ResetBit(kGeoTranslation);
}

////////////////////////////////////////////////////////////////////////////////
/// get the rotation array

const Double_t *TGeoCombiTrans::GetRotationMatrix() const
{
   if (!fRotation) return kIdentityMatrix;
   return fRotation->GetRotationMatrix();
}

/** \class TGeoGenTrans
\ingroup Geometry_classes
Most general transformation, holding a translation, a rotation and a scale
*/

ClassImp(TGeoGenTrans);

////////////////////////////////////////////////////////////////////////////////
/// dummy ctor

TGeoGenTrans::TGeoGenTrans()
{
   SetBit(kGeoGenTrans);
   for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
   fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// constructor

TGeoGenTrans::TGeoGenTrans(const char *name)
             :TGeoCombiTrans(name)
{
   SetBit(kGeoGenTrans);
   for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
   for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
   fRotation = 0;
}

////////////////////////////////////////////////////////////////////////////////
/// constructor

TGeoGenTrans::TGeoGenTrans(Double_t dx, Double_t dy, Double_t dz,
                           Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
             :TGeoCombiTrans("")
{
   SetBit(kGeoGenTrans);
   SetTranslation(dx, dy, dz);
   SetScale(sx, sy, sz);
   SetRotation(rot);
}

////////////////////////////////////////////////////////////////////////////////
/// constructor

TGeoGenTrans::TGeoGenTrans(const char *name, Double_t dx, Double_t dy, Double_t dz,
                           Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
             :TGeoCombiTrans(name)
{
   SetBit(kGeoGenTrans);
   SetTranslation(dx, dy, dz);
   SetScale(sx, sy, sz);
   SetRotation(rot);
}

////////////////////////////////////////////////////////////////////////////////
/// destructor

TGeoGenTrans::~TGeoGenTrans()
{
}

////////////////////////////////////////////////////////////////////////////////
/// clear the fields of this transformation

void TGeoGenTrans::Clear(Option_t *)
{
   memset(&fTranslation[0], 0, kN3);
   memset(&fScale[0], 0, kN3);
   if (fRotation) fRotation->Clear();
}

////////////////////////////////////////////////////////////////////////////////
/// set the scale

void TGeoGenTrans::SetScale(Double_t sx, Double_t sy, Double_t sz)
{
   if (sx<1.E-5 || sy<1.E-5 || sz<1.E-5) {
      Error("ctor", "Invalid scale");
      return;
   }
   fScale[0] = sx;
   fScale[1] = sy;
   fScale[2] = sz;
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoGenTrans::Inverse() const
{
   TGeoHMatrix h = *this;
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// A scale transformation should be normalized by sx*sy*sz factor

Bool_t TGeoGenTrans::Normalize()
{
   Double_t normfactor = fScale[0]*fScale[1]*fScale[2];
   if (normfactor <= 1E-5) return kFALSE;
   for (Int_t i=0; i<3; i++)
      fScale[i] /= normfactor;
   return kTRUE;
}

/** \class TGeoIdentity
\ingroup Geometry_classes
An identity transformation. It holds no data member
and returns pointers to static null translation and identity
transformations for rotation and scale
*/

ClassImp(TGeoIdentity);

////////////////////////////////////////////////////////////////////////////////
/// dummy ctor

TGeoIdentity::TGeoIdentity()
{
   if (!gGeoIdentity) gGeoIdentity = this;
   RegisterYourself();
}

////////////////////////////////////////////////////////////////////////////////
/// constructor

TGeoIdentity::TGeoIdentity(const char *name)
             :TGeoMatrix(name)
{
   if (!gGeoIdentity) gGeoIdentity = this;
   RegisterYourself();
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoIdentity::Inverse() const
{
   TGeoHMatrix h = *gGeoIdentity;
   return h;
}

/** \class TGeoHMatrix
\ingroup Geometry_classes

Matrix class used for computing global transformations
Should NOT be used for node definition. An instance of this class
is generally used to pile-up local transformations starting from
the top level physical node, down to the current node.
*/

ClassImp(TGeoHMatrix);

////////////////////////////////////////////////////////////////////////////////
/// dummy ctor

TGeoHMatrix::TGeoHMatrix()
{
   memset(&fTranslation[0], 0, kN3);
   memcpy(fRotationMatrix,kIdentityMatrix,kN9);
   memcpy(fScale,kUnitScale,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// constructor

TGeoHMatrix::TGeoHMatrix(const char* name)
            :TGeoMatrix(name)
{
   memset(&fTranslation[0], 0, kN3);
   memcpy(fRotationMatrix,kIdentityMatrix,kN9);
   memcpy(fScale,kUnitScale,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// assignment

TGeoHMatrix::TGeoHMatrix(const TGeoMatrix &matrix)
            :TGeoMatrix(matrix)
{
   memset(&fTranslation[0], 0, kN3);
   memcpy(fRotationMatrix,kIdentityMatrix,kN9);
   memcpy(fScale,kUnitScale,kN3);
   if (matrix.IsIdentity()) return;
   if (matrix.IsTranslation())
      SetTranslation(matrix.GetTranslation());
   if (matrix.IsRotation())
      memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
   if (matrix.IsScale())
      memcpy(fScale,matrix.GetScale(),kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// destructor

TGeoHMatrix::~TGeoHMatrix()
{
}

////////////////////////////////////////////////////////////////////////////////
/// assignment

TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix *matrix)
{
   return TGeoHMatrix::operator=(*matrix);
}

////////////////////////////////////////////////////////////////////////////////
/// assignment

TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix &matrix)
{
   if (&matrix == this) return *this;
   Clear();
   Bool_t registered = TestBit(kGeoRegistered);
   TNamed::operator=(matrix);
   if (matrix.IsIdentity()) return *this;
   if (matrix.IsTranslation())
      memcpy(fTranslation,matrix.GetTranslation(),kN3);
   if (matrix.IsRotation())
      memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
   if (matrix.IsScale())
      memcpy(fScale,matrix.GetScale(),kN3);
   SetBit(kGeoRegistered,registered);
   return *this;
}

////////////////////////////////////////////////////////////////////////////////
/// Composition

TGeoHMatrix &TGeoHMatrix::operator*=(const TGeoMatrix &right)
{
   Multiply(&right);
   return *this;
}

TGeoHMatrix TGeoHMatrix::operator*(const TGeoMatrix &right) const
{
   TGeoHMatrix h = *this;
   h *= right;
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// Is-equal operator

Bool_t TGeoHMatrix::operator==(const TGeoMatrix &other) const
{
   if (&other == this) return kTRUE;
   const Double_t *tr = GetTranslation();
   const Double_t *otr = other.GetTranslation();
   for (auto i=0; i<3; i++) if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
   const Double_t *rot = GetRotationMatrix();
   const Double_t *orot = other.GetRotationMatrix();
   for (auto i=0; i<9; i++) if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
   const Double_t *scl = GetScale();
   const Double_t *oscl = other.GetScale();
   for (auto i=0; i<3; i++) if (TMath::Abs(scl[i]-oscl[i])>1.E-10) return kFALSE;
   return kTRUE;
}

////////////////////////////////////////////////////////////////////////////////
/// Fast copy method.

void TGeoHMatrix::CopyFrom(const TGeoMatrix *other)
{
   SetBit(kGeoTranslation, other->IsTranslation());
   SetBit(kGeoRotation, other->IsRotation());
   SetBit(kGeoReflection, other->IsReflection());
   memcpy(fTranslation,other->GetTranslation(),kN3);
   memcpy(fRotationMatrix,other->GetRotationMatrix(),kN9);
}

////////////////////////////////////////////////////////////////////////////////
/// clear the data for this matrix

void TGeoHMatrix::Clear(Option_t *)
{
   SetBit(kGeoReflection, kFALSE);
   if (IsIdentity()) return;
   ResetBit(kGeoTranslation);
   ResetBit(kGeoRotation);
   ResetBit(kGeoScale);
   memcpy(fTranslation,kNullVector,kN3);
   memcpy(fRotationMatrix,kIdentityMatrix,kN9);
   memcpy(fScale,kUnitScale,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Make a clone of this matrix.

TGeoMatrix *TGeoHMatrix::MakeClone() const
{
   TGeoMatrix *matrix = new TGeoHMatrix(*this);
   return matrix;
}

////////////////////////////////////////////////////////////////////////////////
/// Perform a rotation about Z having the sine/cosine of the rotation angle.

void TGeoHMatrix::FastRotZ(const Double_t *sincos)
{
   fRotationMatrix[0] = sincos[1];
   fRotationMatrix[1] = -sincos[0];
   fRotationMatrix[3] = sincos[0];
   fRotationMatrix[4] = sincos[1];
   SetBit(kGeoRotation);
}

////////////////////////////////////////////////////////////////////////////////
/// Return a temporary inverse of this.

TGeoHMatrix TGeoHMatrix::Inverse() const
{
   TGeoHMatrix h;
   h = *this;
   if (IsTranslation()) {
      Double_t tr[3];
      tr[0] = -fTranslation[0]*fRotationMatrix[0] - fTranslation[1]*fRotationMatrix[3] - fTranslation[2]*fRotationMatrix[6];
      tr[1] = -fTranslation[0]*fRotationMatrix[1] - fTranslation[1]*fRotationMatrix[4] - fTranslation[2]*fRotationMatrix[7];
      tr[2] = -fTranslation[0]*fRotationMatrix[2] - fTranslation[1]*fRotationMatrix[5] - fTranslation[2]*fRotationMatrix[8];
      h.SetTranslation(tr);
   }
   if (IsRotation()) {
      Double_t newrot[9];
      newrot[0] = fRotationMatrix[0];
      newrot[1] = fRotationMatrix[3];
      newrot[2] = fRotationMatrix[6];
      newrot[3] = fRotationMatrix[1];
      newrot[4] = fRotationMatrix[4];
      newrot[5] = fRotationMatrix[7];
      newrot[6] = fRotationMatrix[2];
      newrot[7] = fRotationMatrix[5];
      newrot[8] = fRotationMatrix[8];
      h.SetRotation(newrot);
   }
   if (IsScale()) {
      Double_t sc[3];
      sc[0] = 1./fScale[0];
      sc[1] = 1./fScale[1];
      sc[2] = 1./fScale[2];
      h.SetScale(sc);
   }
   return h;
}

////////////////////////////////////////////////////////////////////////////////
/// computes determinant of the rotation matrix

Double_t TGeoHMatrix::Determinant() const
{
   Double_t
   det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
         fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
         fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
         fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
         fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
         fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
   return det;
}

////////////////////////////////////////////////////////////////////////////////
/// multiply to the right with an other transformation
/// if right is identity matrix, just return

void TGeoHMatrix::Multiply(const TGeoMatrix *right)
{
   if (right->IsIdentity()) return;
   const Double_t *r_tra = right->GetTranslation();
   const Double_t *r_rot = right->GetRotationMatrix();
   const Double_t *r_scl = right->GetScale();
   if (IsIdentity()) {
      if (right->IsRotation()) {
         SetBit(kGeoRotation);
         memcpy(fRotationMatrix,r_rot,kN9);
         if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
      }
      if (right->IsScale()) {
         SetBit(kGeoScale);
         memcpy(fScale,r_scl,kN3);
      }
      if (right->IsTranslation()) {
         SetBit(kGeoTranslation);
         memcpy(fTranslation,r_tra,kN3);
      }
      return;
   }
   Int_t i, j;
   Double_t new_rot[9];

   if (right->IsRotation())    {
      SetBit(kGeoRotation);
      if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
   }
   if (right->IsScale())       SetBit(kGeoScale);
   if (right->IsTranslation()) SetBit(kGeoTranslation);

   // new translation
   if (IsTranslation()) {
      for (i=0; i<3; i++) {
         fTranslation[i] += fRotationMatrix[3*i]*r_tra[0]
                         + fRotationMatrix[3*i+1]*r_tra[1]
                         + fRotationMatrix[3*i+2]*r_tra[2];
      }
   }
   if (IsRotation()) {
      // new rotation
      for (i=0; i<3; i++) {
         for (j=0; j<3; j++) {
               new_rot[3*i+j] = fRotationMatrix[3*i]*r_rot[j] +
                                fRotationMatrix[3*i+1]*r_rot[3+j] +
                                fRotationMatrix[3*i+2]*r_rot[6+j];
         }
      }
      memcpy(fRotationMatrix,new_rot,kN9);
   }
   // new scale
   if (IsScale()) {
      for (i=0; i<3; i++) fScale[i] *= r_scl[i];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// multiply to the left with an other transformation
/// if right is identity matrix, just return

void TGeoHMatrix::MultiplyLeft(const TGeoMatrix *left)
{
   if (left == gGeoIdentity) return;
   const Double_t *l_tra = left->GetTranslation();
   const Double_t *l_rot = left->GetRotationMatrix();
   const Double_t *l_scl = left->GetScale();
   if (IsIdentity()) {
      if (left->IsRotation()) {
         if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
         SetBit(kGeoRotation);
         memcpy(fRotationMatrix,l_rot,kN9);
      }
      if (left->IsScale()) {
         SetBit(kGeoScale);
         memcpy(fScale,l_scl,kN3);
      }
      if (left->IsTranslation()) {
         SetBit(kGeoTranslation);
         memcpy(fTranslation,l_tra,kN3);
      }
      return;
   }
   Int_t i, j;
   Double_t new_tra[3];
   Double_t new_rot[9];

   if (left->IsRotation()) {
      SetBit(kGeoRotation);
      if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
   }
   if (left->IsScale())       SetBit(kGeoScale);
   if (left->IsTranslation()) SetBit(kGeoTranslation);

   // new translation
   if (IsTranslation()) {
      for (i=0; i<3; i++) {
         new_tra[i] = l_tra[i]
                      + l_rot[3*i]*  fTranslation[0]
                      + l_rot[3*i+1]*fTranslation[1]
                      + l_rot[3*i+2]*fTranslation[2];
      }
      memcpy(fTranslation,new_tra,kN3);
   }
   if (IsRotation()) {
      // new rotation
      for (i=0; i<3; i++) {
         for (j=0; j<3; j++) {
               new_rot[3*i+j] = l_rot[3*i]*fRotationMatrix[j] +
                                l_rot[3*i+1]*fRotationMatrix[3+j] +
                                l_rot[3*i+2]*fRotationMatrix[6+j];
         }
      }
      memcpy(fRotationMatrix,new_rot,kN9);
   }
   // new scale
   if (IsScale()) {
      for (i=0; i<3; i++) fScale[i] *= l_scl[i];
   }
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about X axis with angle expressed in degrees.

void TGeoHMatrix::RotateX(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = fRotationMatrix[0];
   v[1] = fRotationMatrix[1];
   v[2] = fRotationMatrix[2];
   v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
   v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
   v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
   v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
   v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
   v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];
   memcpy(fRotationMatrix, v, kN9);

   v[0] = fTranslation[0];
   v[1] = c*fTranslation[1]-s*fTranslation[2];
   v[2] = s*fTranslation[1]+c*fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Y axis with angle expressed in degrees.

void TGeoHMatrix::RotateY(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
   v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
   v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
   v[3] = fRotationMatrix[3];
   v[4] = fRotationMatrix[4];
   v[5] = fRotationMatrix[5];
   v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
   v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
   v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];
   memcpy(fRotationMatrix, v, kN9);

   v[0] = c*fTranslation[0]+s*fTranslation[2];
   v[1] = fTranslation[1];
   v[2] = -s*fTranslation[0]+c*fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Rotate about Z axis with angle expressed in degrees.

void TGeoHMatrix::RotateZ(Double_t angle)
{
   SetBit(kGeoRotation);
   Double_t phi = angle*TMath::DegToRad();
   Double_t c = TMath::Cos(phi);
   Double_t s = TMath::Sin(phi);
   Double_t v[9];
   v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
   v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
   v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
   v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
   v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
   v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
   v[6] = fRotationMatrix[6];
   v[7] = fRotationMatrix[7];
   v[8] = fRotationMatrix[8];
   memcpy(&fRotationMatrix[0],v,kN9);

   v[0] = c*fTranslation[0]-s*fTranslation[1];
   v[1] = s*fTranslation[0]+c*fTranslation[1];
   v[2] = fTranslation[2];
   memcpy(fTranslation,v,kN3);
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to YZ.

void TGeoHMatrix::ReflectX(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
   if (leftside) {
      fRotationMatrix[0]=-fRotationMatrix[0];
      fRotationMatrix[1]=-fRotationMatrix[1];
      fRotationMatrix[2]=-fRotationMatrix[2];
   } else {
      fRotationMatrix[0]=-fRotationMatrix[0];
      fRotationMatrix[3]=-fRotationMatrix[3];
      fRotationMatrix[6]=-fRotationMatrix[6];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to ZX.

void TGeoHMatrix::ReflectY(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
   if (leftside) {
      fRotationMatrix[3]=-fRotationMatrix[3];
      fRotationMatrix[4]=-fRotationMatrix[4];
      fRotationMatrix[5]=-fRotationMatrix[5];
   } else {
      fRotationMatrix[1]=-fRotationMatrix[1];
      fRotationMatrix[4]=-fRotationMatrix[4];
      fRotationMatrix[7]=-fRotationMatrix[7];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Multiply by a reflection respect to XY.

void TGeoHMatrix::ReflectZ(Bool_t leftside, Bool_t rotonly)
{
   if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
   if (leftside) {
      fRotationMatrix[6]=-fRotationMatrix[6];
      fRotationMatrix[7]=-fRotationMatrix[7];
      fRotationMatrix[8]=-fRotationMatrix[8];
   } else {
      fRotationMatrix[2]=-fRotationMatrix[2];
      fRotationMatrix[5]=-fRotationMatrix[5];
      fRotationMatrix[8]=-fRotationMatrix[8];
   }
   SetBit(kGeoRotation);
   SetBit(kGeoReflection, !IsReflection());
}

////////////////////////////////////////////////////////////////////////////////
/// Save a primitive as a C++ statement(s) on output stream "out".

void TGeoHMatrix::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
{
   if (TestBit(kGeoSavePrimitive)) return;
   const Double_t *tr = fTranslation;
   const Double_t *rot = fRotationMatrix;
   out << "   // HMatrix: " << GetName() << std::endl;
   out << "   tr[0]  = " << tr[0] << ";    " << "tr[1] = " << tr[1] << ";    " << "tr[2] = " << tr[2] << ";" << std::endl;
   out << "   rot[0] =" << rot[0] << ";    " << "rot[1] = " << rot[1] << ";    " << "rot[2] = " << rot[2] << ";" << std::endl;
   out << "   rot[3] =" << rot[3] << ";    " << "rot[4] = " << rot[4] << ";    " << "rot[5] = " << rot[5] << ";" << std::endl;
   out << "   rot[6] =" << rot[6] << ";    " << "rot[7] = " << rot[7] << ";    " << "rot[8] = " << rot[8] << ";" << std::endl;
   char *name = GetPointerName();
   out << "   TGeoHMatrix *" << name << " = new TGeoHMatrix(\"" << GetName() << "\");" << std::endl;
   out << "   " << name << "->SetTranslation(tr);" << std::endl;
   out << "   " << name << "->SetRotation(rot);" << std::endl;
   if (IsTranslation()) out << "   " << name << "->SetBit(TGeoMatrix::kGeoTranslation);" << std::endl;
   if (IsRotation()) out << "   " << name << "->SetBit(TGeoMatrix::kGeoRotation);" << std::endl;
   if (IsReflection()) out << "   " << name << "->SetBit(TGeoMatrix::kGeoReflection);" << std::endl;
   TObject::SetBit(kGeoSavePrimitive);
}