From 4eea530e03fd8b3a951b2ed2a61ba2be19cb2f52 Mon Sep 17 00:00:00 2001 From: ferdymercury Date: Tue, 22 Apr 2025 15:19:29 +0200 Subject: [PATCH 1/5] [math] implement genvector custom axis rotation Fixes https://root-forum.cern.ch/t/tvector3-rotate-to-arbitrary-rotation-using-xyzvector/63244 --- .../genvector/inc/Math/GenVector/VectorUtil.h | 44 +++++++++++++++++-- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/math/genvector/inc/Math/GenVector/VectorUtil.h b/math/genvector/inc/Math/GenVector/VectorUtil.h index dbd7661c009a4..d3b341a9092f9 100644 --- a/math/genvector/inc/Math/GenVector/VectorUtil.h +++ b/math/genvector/inc/Math/GenVector/VectorUtil.h @@ -279,11 +279,13 @@ namespace ROOT { /** 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 Vector RotateX(const Vector & v, double alpha) { + if (alpha == 0.) + return v; using std::sin; double sina = sin(alpha); using std::cos; @@ -298,11 +300,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 Vector RotateY(const Vector & v, double alpha) { + if (alpha == 0.) + return v; using std::sin; double sina = sin(alpha); using std::cos; @@ -317,11 +321,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 Vector RotateZ(const Vector & v, double alpha) { + if (alpha == 0.) + return v; using std::sin; double sina = sin(alpha); using std::cos; @@ -332,6 +338,38 @@ namespace ROOT { vrot.SetXYZ(x2, y2, v.Z()); 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 + Vector Rotate(const Vector & v, double alpha, const Vector & axis) { + if (alpha == 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; + 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; + 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; + } /** From 3b7bbffd9a646b966ac401eeb0599fbb9999155f Mon Sep 17 00:00:00 2001 From: ferdymercury Date: Tue, 22 Apr 2025 15:21:38 +0200 Subject: [PATCH 2/5] [skip-ci] alignment --- math/genvector/inc/Math/GenVector/VectorUtil.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/math/genvector/inc/Math/GenVector/VectorUtil.h b/math/genvector/inc/Math/GenVector/VectorUtil.h index d3b341a9092f9..f39bdc095b23a 100644 --- a/math/genvector/inc/Math/GenVector/VectorUtil.h +++ b/math/genvector/inc/Math/GenVector/VectorUtil.h @@ -358,7 +358,7 @@ namespace ROOT { const double dy = axis.Y()/ll; const double dz = axis.Z()/ll; 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, + 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; const double xX = v.X(); const double yY = v.Y(); From df7b30a1d4b34bf7242d9fd76b289f0a8273c53b Mon Sep 17 00:00:00 2001 From: ferdymercury Date: Wed, 23 Apr 2025 15:26:33 +0200 Subject: [PATCH 3/5] [math] wrap angle into 2pi as suggested by mdessole --- math/genvector/inc/Math/GenVector/VectorUtil.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/math/genvector/inc/Math/GenVector/VectorUtil.h b/math/genvector/inc/Math/GenVector/VectorUtil.h index f39bdc095b23a..018844adb1a13 100644 --- a/math/genvector/inc/Math/GenVector/VectorUtil.h +++ b/math/genvector/inc/Math/GenVector/VectorUtil.h @@ -284,7 +284,7 @@ namespace ROOT { */ template Vector RotateX(const Vector & v, double alpha) { - if (alpha == 0.) + if (std::fmod(alpha, 2*M_PI ) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -305,7 +305,7 @@ namespace ROOT { */ template Vector RotateY(const Vector & v, double alpha) { - if (alpha == 0.) + if (std::fmod(alpha, 2*M_PI ) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -326,7 +326,7 @@ namespace ROOT { */ template Vector RotateZ(const Vector & v, double alpha) { - if (alpha == 0.) + if (std::fmod(alpha, 2*M_PI ) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -347,7 +347,7 @@ namespace ROOT { */ template Vector Rotate(const Vector & v, double alpha, const Vector & axis) { - if (alpha == 0.) + 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.) From ce3059f9fa324bf089bbb55ee5fdb35f8cca3dad Mon Sep 17 00:00:00 2001 From: ferdymercury Date: Wed, 23 Apr 2025 15:31:43 +0200 Subject: [PATCH 4/5] [math] for consistency, do the same check in TRotation --- math/physics/src/TRotation.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/math/physics/src/TRotation.cxx b/math/physics/src/TRotation.cxx index cbe0a49e26bf4..61764b2587817 100644 --- a/math/physics/src/TRotation.cxx +++ b/math/physics/src/TRotation.cxx @@ -185,6 +185,7 @@ TVector3 class: #include "TRotation.h" #include "TMath.h" #include "TQuaternion.h" +#include 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*M_PI) != 0.) { Double_t ll = axis.Mag(); if (ll == 0.0) { Warning("Rotate(angle,axis)"," zero axis"); From debbea58f97f472d54102e409ce254dc7739ad5a Mon Sep 17 00:00:00 2001 From: mdessole Date: Wed, 30 Apr 2025 15:51:32 +0200 Subject: [PATCH 5/5] Apply formatting --- .../genvector/inc/Math/GenVector/VectorUtil.h | 39 ++++++++++--------- math/physics/src/TRotation.cxx | 2 +- 2 files changed, 21 insertions(+), 20 deletions(-) diff --git a/math/genvector/inc/Math/GenVector/VectorUtil.h b/math/genvector/inc/Math/GenVector/VectorUtil.h index 018844adb1a13..27c92898a2c1f 100644 --- a/math/genvector/inc/Math/GenVector/VectorUtil.h +++ b/math/genvector/inc/Math/GenVector/VectorUtil.h @@ -275,7 +275,6 @@ namespace ROOT { // rotation and transformations - /** rotation along X axis for a generic vector by an Angle alpha returning a new vector. @@ -284,7 +283,7 @@ namespace ROOT { */ template Vector RotateX(const Vector & v, double alpha) { - if (std::fmod(alpha, 2*M_PI ) == 0.) + if (std::fmod(alpha, 2 * M_PI) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -305,7 +304,7 @@ namespace ROOT { */ template Vector RotateY(const Vector & v, double alpha) { - if (std::fmod(alpha, 2*M_PI ) == 0.) + if (std::fmod(alpha, 2 * M_PI) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -326,7 +325,7 @@ namespace ROOT { */ template Vector RotateZ(const Vector & v, double alpha) { - if (std::fmod(alpha, 2*M_PI ) == 0.) + if (std::fmod(alpha, 2 * M_PI) == 0.) return v; using std::sin; double sina = sin(alpha); @@ -338,7 +337,7 @@ namespace ROOT { vrot.SetXYZ(x2, y2, v.Z()); return vrot; } - + /** rotation along a custom axis for a generic vector by an Angle alpha (in rad) returning a new vector. @@ -346,32 +345,34 @@ namespace ROOT { and SetXYZ methods. */ template - Vector Rotate(const Vector & v, double alpha, const Vector & axis) { - if (std::fmod(alpha, 2*M_PI ) == 0.) + 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()); + 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; - 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; + const double dx = axis.X() / ll; + const double dy = axis.Y() / ll; + const double dz = axis.Z() / ll; + 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; 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; + 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); + vrot.SetXYZ(x2, y2, z2); return vrot; } - /** rotation on a generic vector using a generic rotation class. The only requirement on the vector is that implements the diff --git a/math/physics/src/TRotation.cxx b/math/physics/src/TRotation.cxx index 61764b2587817..a65f6fd909ffe 100644 --- a/math/physics/src/TRotation.cxx +++ b/math/physics/src/TRotation.cxx @@ -324,7 +324,7 @@ TRotation::TRotation(const TQuaternion & Q) { /// Rotate along an axis. TRotation & TRotation::Rotate(Double_t a, const TVector3& axis) { - if (std::fmod(a, 2*M_PI) != 0.) { + if (std::fmod(a, 2 * M_PI) != 0.) { Double_t ll = axis.Mag(); if (ll == 0.0) { Warning("Rotate(angle,axis)"," zero axis");