GeoDefinitions - Add function to extract the euler angles. Fix bugs in the intrinsic calculation of GeoRotation. Add corresponding unit test

Hi everybody,

first of all happy new year. My first project of the year is to take the AGDD-based Toroid plugin, load it and then dump the transient GeoModel tree in memory into a GeoModelXML like format. If that's done then we can say, wuuuhu cyaaa unmaintained AGDD code and in fact, any other piece of Geometry can be translated into the same format and we can do the same gestures to the depreciated code as to AGDD. But that's only neccessary, if it's as unmaintained. Anyhow, as a first step, I'd like to extract the six defining paramaters of a GeoTransform. The translation is a no-brainer. For the three Euler angles {\theta}, {\phi} and {\psi} with their associated rotation matrices

R_{x}(\theta) = \begin{pmatrix}
    1 &      0      & 0 \\
    0 & \cos (\theta) & -\sin(\theta) \\
    0 & \sin(\theta)  & \cos(\theta)
\end{pmatrix},
R_{y}(\psi) = \begin{pmatrix}
    \cos (\psi) & 0 & -\sin(\psi) \\
     0            & 1 & 0             \\
    \sin(\psi)  & 0 & \cos(\psi)
\end{pmatrix},
R_{z}(\phi) = \begin{pmatrix}     
     \cos (\phi) & -\sin(\phi)   & 0  \\    
     \sin(\phi)  & \cos(\phi)    & 0  \\
        0        &    0          & 1
\end{pmatrix},

a bit more work is needed. The GeoRotation implements kind of a concatenation of the rotation matrices, where I'm currently trying to figure out the respective order. If I interpret the constrcutor

  GeoRotation::GeoRotation(double phi, double theta, double psi)
  {
    double sinPhi   = std::sin( phi   ), cosPhi   = std::cos( phi   );
    double sinTheta = std::sin( theta ), cosTheta = std::cos( theta );
    double sinPsi   = std::sin( psi   ), cosPsi   = std::cos( psi   );
    
    this->operator()(0,0) =   cosPsi * cosPhi - cosTheta * sinPhi * sinPsi;
    this->operator()(0,1) =   cosPsi * sinPhi + cosTheta * cosPhi * sinPsi;
    this->operator()(0,2) =   sinPsi * sinTheta;
    
    this->operator()(1,0) = - sinPsi * cosPhi - cosTheta * sinPhi * cosPsi;
    this->operator()(1,1) = - sinPsi * sinPhi + cosTheta * cosPhi * cosPsi;
    this->operator()(1,2) =   cosPsi * sinTheta;
    
    this->operator()(2,0) =   sinTheta * sinPhi;
    this->operator()(2,1) = - sinTheta * cosPhi;
    this->operator()(2,2) =   cosTheta;
  }

correctly, then

 \text{GeoRot}(\theta,\psi,\phi) = \begin{pmatrix} 
   \cos(\psi)\cos(\phi) - \cos(\theta)\sin(\phi)\sin(\psi)  & -\sin(\psi) \cos(\phi) -\cos(\theta)\sin(\phi)\cos(\psi) & \sin(\theta)\sin(\phi) \\
    \cos(\psi) \sin(\phi) + \cos(\theta)\cos(\phi)\sin(\psi) & -\sin(\psi)\sin(\phi) + \cos(\theta)\cos(\phi) \cos(\psi) & -\sin(\theta)\cos(\phi) \\
   \sin(\psi)\sin(\theta) & \cos(\psi)\sin(\theta) & \cos(\theta)
\end{pmatrix}

is the corresponding matrix equation. Next, let's see which of the three {R_{i}} is obtained by setting two of the parameters to zero.

  \text{GeoRot}(0, 0,\phi) = \begin{pmatrix} 
   \cos(0)\cos(\phi) - \cos(0)\sin(\phi)\sin(0)  & -\sin(0) \cos(\phi) -\cos(0)\sin(\phi)\cos(0) & \sin(0)\sin(\phi) \\
    \cos(0) \sin(\phi) + \cos(0)\cos(\phi)\sin(0) & -\sin(0)\sin(\phi) + \cos(0)\cos(\phi) \cos(0) & -\sin(0)\cos(\phi) \\
   \sin(0)\sin(0) & \cos(0)\sin(0) & \cos(0)
\end{pmatrix}  = \begin{pmatrix}     
     \cos (\phi) & -\sin(\phi)   & 0  \\    
     \sin(\phi)  & \cos(\phi)    & 0  \\
        0        &    0          & 1
\end{pmatrix} = R_{z}(\phi)

okay, {\phi} is a rotation around the z-axis, Let's take a look at {\theta}:

 \text{GeoRot}(\theta,0, 0) = \begin{pmatrix} 
   \cos(0)\cos(0) - \cos(\theta)\sin(0)\sin(0)  & -\sin(0) \cos(0) -\cos(\theta)\sin(0)\cos(0) & \sin(\theta)\sin(0) \\
    \cos(0) \sin(0) + \cos(\theta)\cos(0)\sin(0) & -\sin(0)\sin(0) + \cos(\theta)\cos(0) \cos(0) & -\sin(\theta)\cos(0) \\
   \sin(0)\sin(\theta) & \cos(0)\sin(\theta) & \cos(\theta)
\end{pmatrix} = \begin{pmatrix}
    1 &      0      & 0 \\
    0 & \cos (\theta) & -\sin(\theta) \\
    0 & \sin(\theta)  & \cos(\theta)
\end{pmatrix} = R_{x}(\theta)

Indeed, that's the rotation around the x-axis. Last but not least, the angle {\psi} will be then the y-axis rotation.

\text{GeoRot}(0,\psi,0) = \begin{pmatrix} 
   \cos(\psi)\cos(0) - \cos(0)\sin(0)\sin(\psi)  & -\sin(\psi) \cos(0) -\cos(0)\sin(0)\cos(\psi) & \sin(0)\sin(0) \\
    \cos(\psi) \sin(0) + \cos(0)\cos(0)\sin(\psi) & -\sin(\psi)\sin(0) + \cos(0)\cos(0) \cos(\psi) & -\sin(0)\cos(0) \\
   \sin(\psi)\sin(0) & \cos(\psi)\sin(0) & \cos(0)
\end{pmatrix}  = 
 \begin{pmatrix} 
   \cos(\psi)  & -\sin(\psi)  & 0 \\
    \sin(\psi) &  \cos(\psi)  & 0\\
   0           & 0            & 1
\end{pmatrix} = R_{x}(\psi)

That's another x-axis rotation. I was quite confused about this result, so I've written a quick unit test to repeat the exercise above

1: testEulerAngles() 19 The GeoTrf::GeoRotation around the x-axis with angle 315 
1:  0.707107 -0.707107         0
1:  0.707107  0.707107         0
1:        -0        -0         1
1: testEulerAngles() 27 The GeoTrf::GeoRotation around the y-axis with angle 315 
1:         1         0        -0
1:        -0  0.707107 -0.707107
1:        -0  0.707107  0.707107
1: testEulerAngles() 35 The GeoTrf::GeoRotation around the z-axis with angle 315 
1:  0.707107 -0.707107        -0
1:  0.707107  0.707107         0
1:         0        -0         1

Indeed, in 2 out of the three cases, the rotation matrix is actually the same. In athena, GeoRotations around the three 'axes' are used in the inner detector GeoModel. There're also three usages of the GeoRotation in GeoModelATLAS.

There's also another interesting bug that setting 2 angles to non-zero does not neccessarily construct a valid rotation matrix

Tagging: @boudreau, @tsulaia, @nstyles, @jchapman, @jcatmore, @mbandier, @sroe

Wolframalpha

Wolframalpha

Wolframalpha

Wolframalpha

Wolframalpha

Wolframalpha