TVector3::RotateUz

General Chat
1

I recently stumbled over the method TVector3::RotateUz(), which did not work as expected for us. I understand that the class TVector3 is legacy, but it is still in the ROOT distribution, and we are making use of it.

The documentation says:

v1.RotateUz(direction) transforms from the rotated frame (z’ parallel to direction…) to the (x,y,z,) frame.

Implementation in TVector3.cxx:

void TVector3::RotateUz(const TVector3& NewUzVector) {
   Double_t u1 = NewUzVector.fX;
   Double_t u2 = NewUzVector.fY;
   Double_t u3 = NewUzVector.fZ;
   Double_t up = u1*u1 + u2*u2;
 
   if (up) {
      up = TMath::Sqrt(up);
      Double_t px = fX,  py = fY,  pz = fZ;
      fX = (u1*u3*px - u2*py + u1*up*pz)/up;
      fY = (u2*u3*px + u1*py + u2*up*pz)/up;
      fZ = (u3*u3*px -    px + u3*up*pz)/up;
   } else if (u3 < 0.) { fX = -fX; fZ = -fZ; }      // phi=0  teta=pi
   else {};
}

Now consider the special case when NewUzVector = (u1, 0, u3) is in the x-z plane. Then we should have a rotation around the y = y’ axis, which should give

fX = u3 * px + u1 * pz
fY = py
fZ = -u1 * px * u3 * pz

But from the implementation above, we get with up = sqrt(u1*u1) = abs(u1):

fX = u1 / abs(u1) * u3 * px + u1 * pz
fY = u1 / abs(u1) * py
fZ = -abs(u1) * px + u3 * pz  (exploiting the normalisation, i.e. u1*u1+u3*u3 = 1)

This corresponds to the expectation only if u1 >= 0, i.e., for rotations between 0 and pi. In other cases, fX and fY are mirrored (rotated by pi) in the xy plane.

Can you confirm my findings, or do I misinterpret the method RotateUz?

2

Hi @friese,

thank you for your question. I would strongly recommend that you try to avoid using the TVector3 legacy class and instead try with ROOT::Math::XYZVector.

While I tag @moneta for some further insights, I invite you to have a look at this thread from a few years ago, where a similar issue was discussed: Rotating vectors with ROOT::Math::XYZVector.

Cheers,
Marta

3

Dear Marta,

thanks for the reply. Yes, we will use XYZVector for new implementations, but we also have to deal with legacy code (in this case inherited through FairRoot).

Thanks also for pointing to the other thread, which I saw before but which does not answer to my problem. Surely, there are ways to achieve the functionality of RotateUz without using the method itself, e.g., using TVector3::Rotate:

TVector3 rot_axis(-fBeamDirection.Unit().Y(), fBeamDirection.Unit().X(), 0);
auto rot_angle = TMath::ACos(fBeamDirection.Unit().Z());
mom.Rotate(rot_angle, rot_axis);

or the equivalent with XYZVector.

I posted this here primarily to ascertain the implementation to be faulty, and if yes, give other users the chance to avoid it in case they encounter the same problem.

Best wishes,

Volker

1 Like
4

Wrt meaning of RotateUz: Meaning of RotateUz - #9 by xavierPrudent

If there is a bug in this function, it should be also reported to the Geant4 version:

5

My hypothesis is that RotateUz only guarantees that the new coordinate frame is the ‘local frame’, ie z’ is aligned to the direction, the x’ and y’ relative position in the local frame do not “change the physics”? So y = y’ is not guaranteed even if u2=0.
The only condition fulfilled is that x’ is in the theta-plane and that y’ is in xy plane and perpendicular to theta plane.

see [skip-ci] Improve RotateUz documentation by ferdymercury · Pull Request #18857 · root-project/root · GitHub and rotateUz looks fishy - #4 by ferhue - Physics Processes, Models and Cross Sections - Geant4 Forum

6

Thank for looking into the matter and communicating also in the GEANT4 forum. I also tried to get deeper insight.

Maybe it is worthwile to explain where we are using this function and stumbled over its behaviour. When transporting an input event with GEANT, we have to rotate the event (all particle momenta within the event) from the coordinate system of the event generator (UrQMD in our case), which usually defines the z axis as beam axis, to the experiment coordinate system, accounting for a finite beam divergence. For this, we apply TVector3::RotateUz to each particle momentum, with the beam direction as argument, before we put them onto the GEANT stack.

Now, the implementation of the function uses two rotations. In polar coordinates, the direction vector NewUzVector can be written as u = (sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)), with 0 <= theta <= pi, and -pi < phi < pi. This vector can be rotated into the unit vector in z by first rotating around the z axis by -phi, which brings u into the x-z plane, and then rotating around the (new) y axis by theta.

What we want is the inverse transformation, so first rotating by -theta around y, and then by phi around z. The implementation follows this as @ferhue has explained in the GEANT4 forum. So, it appears mathematically correct.

But is does not behave as we expect, as can be seen in the special case of e.g., u2 = 0, i.e., u is in the x-z plane. Looking into the formulae, one gets for the transformation matrix in the limit of theta->0

A = (
1  0  0
0  1  0
0  0  1
)

i.e. identity for ux > 0 (phi = 0), but

A = (
-1  0  0
 0 -1  0
 0  0  1
)

for x <0 (phi = pi). So here, the vector is rotated by pi in x-y.

In other words, the transformation has a discontinuity at theta=0. But rotations should be continuous transformations, should they not?

Probably, the limes theta-> 0 is mathematically not legal, since for theta = 0, phi is not defined. But the fact remains that one gets very distinct transformations for infinitesimally small ux depending on the sign of ux.

What one surely expects or wants to have in the case of u being in the x-z plane is just a rotation by theta around the y axis. The function RotateUz does not provide this. To me, it seems an artifact caused by restricting theta to positive values; so, if ux < 0, the rotation by pi in x-y is applied to get theta positive.

One remark still: This is a purely mathematical function, so arguments like “rotations in x-y do not change the physics” should not apply.

7

I agree with you, the discontinuity stems from the fact that they programmed the matrix
B as non-canonical rotation matrix with elements 12 and 21 being |sin(angle)| and -|sin(angle)|, respectively.

So in short, it’s not a fully geometrically interpretable rotation matrix…

The question is now how to proceed:

a) Document better the function warning the user that this is not a (fully) geometrically interpretable rotation matrix. It just ensures that the vector is expressed with the z’ aligned to direction, and x and y are not guaranteed to follow from a senseful rotation. Likewise, one cannot define the inverse rotation matrix, the inverse transformation is not gonna be a angle-standard rotation matrix either.
Suggest using a different function called RotateUz2 that does things better if you want to interpret the chain of rotation angle multiplications that led to this result

or …

b) Convincing Geant4 to reimplement this function differently and then backport the changes to ROOT.

@mdessole what’s your opinion?

8

Btw this reminds me a bit about the gimbal lock Xsens Knowledge base

9

The puzzling thing is that the applied transformation matrix

u1*u3/up    -u2/up    u1
u2*u3/up     u1/up    u2
-up          0        u3

is orthonormal and its determinant is unity - so it IS a rotation matrix. Still it has this weird behaviour…

10

Yes, I agree that the matrix is a mathematical rotation matrix, but it violates the interpretation as a geometrical rotation with an angle theta, since we have for matrix B:
R(-theta) = R(theta)
whereas for a standard rotation matrix, we would have:
R(-theta) = R_transposed(theta)

11

If I remember correctly, rotateUz() was added to CLHEP for needs of Geant4. This functions transfer coordinates of a point/vector from its current coordinate system to a new coordinate system where Z-axis is directed along NewUzVector.

This operation is not strictly defined: there is some uncertainty about the direction of the new X and Y axes. The rotation matrix used in the current implementation is a true rotation matrix, so I don’t see any compelling reason to revise the function.

12

I agree, thanks for the feedback. Clarifying the documentation is thus the best way to go.

13
14

Thanks for the feedback and for taking care. I see that the documentation of TVector3::RotateUz is already updated.

For our purposes, we will not use this function, then. The “correct” transformation (shortest rotation path) is provided by rotation around the unit vector which is perpendicular to u and the z axis: (u2/up, -u1/up, 0) by the angle -theta, where theta = arccos(u3). This can be achieved by using TVector3::Rotate(angle, axis). In the case of u2 = 0, the transformation reduces to the desired one:

 u3    0   u1
  0    1    0
-u1    0   u3
15

Maybe one could add the description of this “correct” transformation to the official “documentation of TVector3::RotateUz” (so that other people will be able to find it easily).

16

see [skip-ci][math] mention gimbal lock by ferdymercury · Pull Request #18987 · root-project/root · GitHub