Skip to content

[math] implement genvector custom axis rotation #18462

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 41 additions & 3 deletions math/genvector/inc/Math/GenVector/VectorUtil.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 <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;
Expand All @@ -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 <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;
Expand All @@ -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 <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;
Expand All @@ -332,6 +338,38 @@ namespace ROOT {
vrot.SetXYZ(x2, y2, v.Z());
return vrot;
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not implementing an access operator in the Rotation3D class instead?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see e.g. #18480

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Monica for the swift review.
Wrt using the Rotation3D class: It's maybe a bit less efficient to define an AxisAngle, then a TRotation3D, and finally calling Rotate, than doing the direct multiplication?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You have a point. However, I agree with the author of the forum post you mentioned, it's counter-intuitive that a TRotation3D does not work with Rotate...

/**
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;
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;
}


/**
Expand Down
3 changes: 2 additions & 1 deletion math/physics/src/TRotation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,7 @@ TVector3 class:
#include "TRotation.h"
#include "TMath.h"
#include "TQuaternion.h"
#include <cmath>

ClassImp(TRotation);

Expand Down Expand Up @@ -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");
Expand Down
Loading