rotateUz looks fishy

[ferhue]

Geant4 implements a rotateUz function, which has been copy-pasted into ROOT (some reorganization of factors apart).

A ROOT user found out that the implementation of this function does not take into account the ‘quadrant’, ie the result might be mirrored if the x direction is negative and y direction is 0.
see: TVector3::RotateUz - General Chat - ROOT Forum

Could a developer look at it and confirm this extent?
Here is the relevant snippet:

Blaming geant4/source/externals/clhep/src/ThreeVector.cc at master ·...

Alternatively: could someone described what the equation was to reach this formula? What was the product of two individual rotation matrices assumed to reach that code implementation? Thanks in advance.

[maire1]

I think it is a theta rotation around y axis, followed by a phi rotation around z axis.
Theta is assumed to be in [0, pi] , phi in [0, twopi]

I used examples/extended/eventgenerator/particleGun, case 3, to check that rotateUz() works properly

[ferhue]

Thanks a lot for checking. Let’s assume phi=0, ie uy=0..

The question: why is then theta restricted to 0 to pi? And not -pi to pi.
It’s not the same to describe a vector in polar coordinates than to perform a rotation of one full reference frame in 3D to another I guess?

In other words, i guess because in a local frame it does not matter how X and Y are oriented, I think that rotateUz does a two rotations plus one XY reflection (sometimes).

[ferhue]

I found out the equation behind the scenes:

C = A x B

A = 
u1/up  -u2/up 0
u2/up   u1/up 0
0       0     1

B =
u3  0 up
0   1 0
-up 0 u3

But up is always a positive number (theta is in 0, pi, as you said). So the sign of u1 and u2 is totally ignored for creating the B rotation (not for A).

In the case of u1=0, I do not see clearly if this rotation makes sense. (Why should it be rotated around y-axis?)

In the case of u2=0, and u1<0, A becomes a rotation by pi but B stays the same independently of u1<0 or u1>0, ie B is rotating in the ‘wrong sense’ of the y-axis right-hand rule.

[ferhue]

Related: [math] Improve TVector3::RotateUz documentation · root-project/root@0ae85ac · GitHub

[ferhue]

Related post about the discontinuity bug of this pseudo-rotation matrix:

TVector3::RotateUz

@gcosmo @bmorgan @atolosad @mkelsey what’s your opinion?

Could this behavior be better documented?

[mkelsey]

This feels to me like it’s nothing more than gimbal lock, which is inescapable with polar coodinates. Phi is intrinsically undefined at theta=0, etc. etc. You’re right that the behaviour could be better documented, but it’s not obvious to me at least where that documention would/should live.

[gcosmo]

rotateUz is a method of Hep3Vector in CLHEP. Its documentation can be found in the CLHEP User Guide, so any update should eventually be applied there.

[evc]

The rotation matrix used in this function is:

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

Let me explain how it was calculated. The columns of a rotation matrix represent the coordinates of the new axes: the first column is the new X-axis, the second is the new Y-axis , the third is the new Z-axis.

In the case of rotateUz(), the new Z was specified, so the third column is (u1,u2,u3). To fill the elements of the second column we can take any unit vector orthogonal to the new Z-axis. After that, the first column can be calculated as a cross product of the second and the third columns.

New Y-axis is calculated as a cross product of the old and new Z axes:

(0,0,1)x(u1,u2,u3) = (-u2,u1,0)

after normalisation it becomes (-u2/up,u1/up,0)

Respectively the new X-axis:

(-u2/up,u1/up,0)x(u1,u2,u3) = (u1*u3/up,u2*u3/up,-up)

There is a special case when the new Z-axis is collinear to the old Z-axis. This case is considered separately.

[ferhue]

Thanks a lot! That makes it much clearer than speaking about the composition of two rotations.