[math] implement genvector custom axis rotation (root-project#18462)

* [math] implement genvector custom axis rotation

Fixes https://root-forum.cern.ch/t/tvector3-rotate-to-arbitrary-rotation-using-xyzvector/63244

* [skip-ci] alignment

* [math] wrap angle into 2pi

as suggested by mdessole

* [math] for consistency, do the same check in TRotation

* Apply formatting

* [nfc] format 3d matrix aligned

* [test] add test case for GenVector Rotate around axis by angle

as suggested by hageboeck

* [test] fix dupe variable

* [test] fix variable unused

* [math] Update pi computation for Windows build

Co-authored-by: Bertrand Bellenot <bellenot@users.noreply.github.com>

---------

Co-authored-by: mdessole <monicadessole@gmail.com>
Co-authored-by: Monica Dessole <36501030+mdessole@users.noreply.github.com>
Co-authored-by: Bertrand Bellenot <bellenot@users.noreply.github.com>

Author: 3 people

math/genvector/inc/Math/GenVector/VectorUtil.h

@@ -275,15 +275,16 @@ namespace ROOT {
 
          // rotation and transformations
 
-
          /**
           rotation along X axis for a generic vector by an Angle alpha
           returning a new vector.
-          The only pre requisite on the Vector is that it has to implement the X() , Y() and Z()
+          The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
           and SetXYZ methods.
           */
          template <class Vector>
          Vector RotateX(const Vector & v, double alpha) {
+            if (std::fmod(alpha, 2 * M_PI) == 0.)
+               return v;
             using std::sin;
             double sina = sin(alpha);
             using std::cos;
@@ -298,11 +299,13 @@ namespace ROOT {
          /**
           rotation along Y axis for a generic vector by an Angle alpha
           returning a new vector.
-          The only pre requisite on the Vector is that it has to implement the X() , Y() and Z()
+          The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
           and SetXYZ methods.
           */
          template <class Vector>
          Vector RotateY(const Vector & v, double alpha) {
+            if (std::fmod(alpha, 2 * M_PI) == 0.)
+               return v;
             using std::sin;
             double sina = sin(alpha);
             using std::cos;
@@ -317,11 +320,13 @@ namespace ROOT {
          /**
           rotation along Z axis for a generic vector by an Angle alpha
           returning a new vector.
-          The only pre requisite on the Vector is that it has to implement the X() , Y() and Z()
+          The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
           and SetXYZ methods.
           */
          template <class Vector>
          Vector RotateZ(const Vector & v, double alpha) {
+            if (std::fmod(alpha, 2 * M_PI) == 0.)
+               return v;
             using std::sin;
             double sina = sin(alpha);
             using std::cos;
@@ -333,6 +338,40 @@ namespace ROOT {
             return vrot;
          }
 
+         /**
+          rotation along a custom axis for a generic vector by an Angle alpha (in rad)
+          returning a new vector.
+          The only pre requisite on the Vector is that it has to implement the X(), Y() and Z()
+          and SetXYZ methods.
+          */
+         template <class Vector>
+         Vector Rotate(const Vector &v, double alpha, const Vector &axis)
+         {
+            if (std::fmod(alpha, 2 * M_PI) == 0.)
+               return v;
+            const double ll = std::sqrt(axis.X() * axis.X() + axis.Y() * axis.Y() + axis.Z() * axis.Z());
+            if (ll == 0.)
+               GenVector::Throw("Axis Vector has zero magnitude");
+            const double sa = std::sin(alpha);
+            const double ca = std::cos(alpha);
+            const double dx = axis.X() / ll;
+            const double dy = axis.Y() / ll;
+            const double dz = axis.Z() / ll;
+            // clang-format off
+            const double rot00 = (1 - ca) * dx * dx + ca     , rot01 = (1 - ca) * dx * dy - sa * dz, rot02 = (1 - ca) * dx * dz + sa * dy,
+                         rot10 = (1 - ca) * dy * dx + sa * dz, rot11 = (1 - ca) * dy * dy + ca     , rot12 = (1 - ca) * dy * dz - sa * dx,
+                         rot20 = (1 - ca) * dz * dx - sa * dy, rot21 = (1 - ca) * dz * dy + sa * dx, rot22 = (1 - ca) * dz * dz + ca     ;
+            // clang-format on
+            const double xX = v.X();
+            const double yY = v.Y();
+            const double zZ = v.Z();
+            const double x2 = rot00 * xX + rot01 * yY + rot02 * zZ;
+            const double y2 = rot10 * xX + rot11 * yY + rot12 * zZ;
+            const double z2 = rot20 * xX + rot21 * yY + rot22 * zZ;
+            Vector vrot;
+            vrot.SetXYZ(x2, y2, z2);
+            return vrot;
+         }
 
          /**
           rotation on a generic vector using a generic rotation class.

math/genvector/test/testGenVector.cxx

@@ -466,6 +466,21 @@ int testRotations3D() {
   iret |= compare(qr1.Y(), qr2.Y(),"y diff",10 );
   iret |= compare(qr1.Z(), qr2.Z(),"z diff",10 );
 
+  // test TVector3-like Rotate function around some axis by an angle. Test case taken from cwessel:
+  // https://root-forum.cern.ch/t/tvector3-rotate-to-arbitrary-rotation-using-xyzvector/63244/7
+  XYZVector ag(17, -4, -23);
+  double angle = 100;
+  XYZVector axisg(-23.4, 1.7, -0.3);
+  XYZVector rotated = Rotate(ag, angle, axisg);
+  // should be equivalent to:
+  // TVector3 at(17, -4, -23);
+  // TVector3 axist(-23.4, 1.7, -0.3);
+  // at.Rotate(angle, axist);
+  // at.Print();
+  // (17.856456,8.106555,-21.199782)
+  iret |= compare(rotated.X(), 17.856456, "x diff", 1e10);
+  iret |= compare(rotated.Y(), 8.106555, "y diff", 1e10);
+  iret |= compare(rotated.Z(), -21.199782, "z diff", 1e10);
 
   if (iret == 0) std::cout << "\tOK\n";
   else std::cout << "\t FAILED\n";

math/physics/src/TRotation.cxx

@@ -185,6 +185,7 @@ TVector3 class:
 #include "TRotation.h"
 #include "TMath.h"
 #include "TQuaternion.h"
+#include <cmath>
 
 ClassImp(TRotation);
 
@@ -323,7 +324,7 @@ TRotation::TRotation(const TQuaternion & Q) {
 /// Rotate along an axis.
 
 TRotation & TRotation::Rotate(Double_t a, const TVector3& axis) {
-   if (a != 0.0) {
+   if (std::fmod(a, 2 * TMath::Pi()) != 0.) {
       Double_t ll = axis.Mag();
       if (ll == 0.0) {
          Warning("Rotate(angle,axis)"," zero axis");