Wrong result when ROOT::Math::AxisAngle applied to ROOT::Math::Polar3D<double>

[karuboniru]

Check duplicate issues.

• Checked for duplicates

Description

I identified this issue when trying to re-write some of my code using GenVector, and the result seems strange that I identified that it seems that the rotation result don't make sense.

In the following code, I am creating different but equal vectors in different coordinate system and apply a same rotation to both. The results after rotation are not the same, which I think should be caused a bug (or I am having a very wrong assumption of how the GenVector interfaces works? It would be also helpful if you can point out the issue even it is on my side.)

Reproducer

#include <Math/AxisAngle.h>
#include <Math/Polar3D.h>
#include <Math/Vector3D.h>
#include <iostream>
#include <sstream>

std::string format_vec(const auto &vec) {
  std::stringstream ss;
  ss << "(" << vec.x() << ", " << vec.y() << ", " << vec.z() << ")"
     << " R = " << vec.R();
  return ss.str();
}

int main() {
  // Create a rotation transformation using axis-angle representation
  const ROOT::Math::AxisAngle the_rec_rot(ROOT::Math::XYZVector{1, 1, 0}, 0.5);

  // Define a vector in polar coordinates and use it to initialize
  // an equivalent vector in Cartesian coordinates
  const ROOT::Math::Polar3D<double> to_rot_polar{1, 1, 1};
  const ROOT::Math::XYZVector to_rot_cartesian{to_rot_polar};

  std::cout << "Polar Vector: " << format_vec(to_rot_polar) << "\n"
            << "Cartesian Vector: " << format_vec(to_rot_cartesian) << '\n';

  // Apply the rotation transformation to both vector representations
  // The results should be identical
  auto rot_polar = the_rec_rot(to_rot_polar);
  auto rot_cartesian = the_rec_rot(to_rot_cartesian);

  // But the results are not identical...
  std::cout << "Rotated Polar: " << format_vec(rot_polar) << '\n'
            << "Rotated Cartesian: " << format_vec(rot_cartesian) << "\n";
}

c++ $(root-config --cflags --glibs) test.cxx -o test -O2 -lGenVector && ./test

Get output

Polar Vector: (0.454649, 0.708073, 0.540302) R = 1
Cartesian Vector: (0.454649, 0.708073, 0.540302) R = 1
Rotated Polar: (0.269919, 0.169252, 0.570379) R = 0.653326
Rotated Cartesian: (0.653326, 0.509396, 0.560072) R = 1

ROOT version

root --version                                                                                                           
ROOT Version: 6.32.02
Built for linuxx8664gcc on Jun 19 2024, 05:05:49
From heads/master@tags/v6-32-02

Installation method

/cvmfs/sft.cern.ch/lcg/views/LCG_106/x86_64-el9-gcc13-opt/setup.sh

Operating system

Red Hat Enterprise Linux release 9.6 (Plow)

Additional context

I did not test more rotation classes.

No response

[ferdymercury]

The bug in ROOT is that the rotation is done correctly, but then the return value is a:
const ROOT::Math::Polar3D<double>{0.653326, 0.509396, 0.560072};
which then results in:
(0.269919, 0.169252, 0.570379) R = 0.653326
whereas the function should return rather:
const ROOT::Math::Polar3D<double> ReexpressInPolarR_Theta_Phi({0.653326, 0.509396, 0.560072});

See AxisAngle.h line 202

Pinging @mdessole and @hageboeck for ideas.

[karuboniru]

See AxisAngle.h line 202

Seems more than this, at least rotation with ROOT::Math::Rotation3D is also affected (or they have some code in common?)

#include <Math/AxisAngle.h>
#include <Math/Polar3D.h>
#include <Math/Vector3D.h>
#include <iostream>
#include <sstream>

std::string format_vec(const auto &vec) {
  std::stringstream ss;
  ss << "(" << vec.x() << ", " << vec.y() << ", " << vec.z() << ")"
     << " R = " << vec.R();
  return ss.str();
}

int main() {
  // Create a rotation transformation using axis-angle representation
  const ROOT::Math::AxisAngle the_rec_rot(ROOT::Math::XYZVector{1, 1, 0}, 0.5);
  const ROOT::Math::Rotation3D the_rec_rot_general{the_rec_rot};

  // Define a vector in polar coordinates and use it to initialize
  // an equivalent vector in Cartesian coordinates
  const ROOT::Math::Polar3D<double> to_rot_polar{1, 1, 1};
  const ROOT::Math::XYZVector to_rot_cartesian{to_rot_polar};

  std::cout << "Polar Vector:      " << format_vec(to_rot_polar) << "\n"
            << "Cartesian Vector:  " << format_vec(to_rot_cartesian) << '\n';

  // Apply the rotation transformation to both vector representations
  // The results should be identical
  const auto rot_polar = the_rec_rot_general(to_rot_polar);
  const auto rot_cartesian = the_rec_rot_general(to_rot_cartesian);

  // But the results are not identical...
  std::cout << "Rotated Polar:     " << format_vec(rot_polar) << '\n'
            << "Rotated Cartesian: " << format_vec(rot_cartesian) << "\n";
}

[karuboniru]

Thing like this seems everywhere (multiple times)

root/math/genvector/inc/Math/GenVector/AxisAngle.h

Lines 239 to 240 in a21a966

	DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);
	return ForeignVector ( rxyz.X(), rxyz.Y(), rxyz.Z() );

https://github.com/search?q=repo%3Aroot-project%2Froot%20ForeignVector%20(%20rxyz.X()%2C%20rxyz.Y()%2C%20rxyz.Z()%20)%3B&type=code

   template <class ForeignVector>
   ForeignVector
   operator() (const  ForeignVector & v) const {
      DisplacementVector3D< Cartesian3D<double> > xyz(v);
      DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);
      return ForeignVector ( rxyz.X(), rxyz.Y(), rxyz.Z() );
   }

I think it should be something like

   template <class ForeignVector>
   ForeignVector
   operator() (const  ForeignVector & v) const {
      DisplacementVector3D< Cartesian3D<double> > xyz(v);
      DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);
      return ForeignVector ( rxyz );
   }

As long as ForeignVector can always be initialized from DisplacementVector3D

[ferdymercury]

The issue is that ForeignVector is defined as:

      Preconditions:  v must implement methods x(), y(), and z()
      and the arbitrary vector type must have a constructor taking (x,y,z)

If we now change the precondition to also implement a constructor taking (v) then that might break some user scripts I guess?

[karuboniru]

At least for me I'd rather face a compilation failure than random wrong result 🤣

[ferdymercury]

Actually, the execution is following this code path:

so it directly thinks it's an XYZVector rather than a ForeignVector

[karuboniru]

Maybe its just the template gets inlined and escaped the debugger's eyes? And the code path you are looking at is during DisplacementVector3D< Cartesian3D<double> > rxyz = operator()(xyz);?

[ferdymercury]

Yes, I think you're right