Skip to content

Commit

Permalink
refactor: FPE safe eta/theta conversion (acts-project#3788)
Browse files Browse the repository at this point in the history 
<!-- This is an auto-generated comment: release notes by coderabbit.ai -->
## Summary by CodeRabbit

- **New Features**
	- Introduced a new structure for angle conversion traits, enhancing precision and safety for angle calculations.
	- Added safety checks in angle conversion functions to handle edge cases more gracefully.

- **Bug Fixes**
	- Updated angle conversion functions to prevent floating-point exceptions by implementing bounds checking.

- **Tests**
	- Added new data-driven test cases to validate the robustness of angle conversion functions, ensuring consistent and valid outputs across a range of inputs.
<!-- end of auto-generated comment: release notes by coderabbit.ai -->
  • Loading branch information
@andiwand @Rosie-Hasan
andiwand authored and Rosie-Hasan committed Nov 13, 2024
1 parent 39815c6 commit 37db0d4
Show file tree
Hide file tree
Showing 3 changed files with 148 additions and 41 deletions.
44 changes: 43 additions & 1 deletion Core/include/Acts/Utilities/AngleHelpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,26 +9,68 @@
#pragma once

#include <cmath>
#include <numbers>

namespace Acts::AngleHelpers {

template <typename Scalar>
struct EtaThetaConversionTraits {};

template <>
struct EtaThetaConversionTraits<float> {
static constexpr float minTheta = 1e-6f;
static constexpr float maxTheta = std::numbers::pi_v<float> - minTheta;

static constexpr float maxEta = 80.0f;
static constexpr float minEta = -maxEta;
};

template <>
struct EtaThetaConversionTraits<double> {
static constexpr double minTheta = 1e-12;
static constexpr double maxTheta = std::numbers::pi - minTheta;

static constexpr double maxEta = 700.0;
static constexpr double minEta = -maxEta;
};

/// Calculate the pseudorapidity from the polar angle theta.
///
/// This function aims to be FPE safe and returns infinity for theta values
/// outside the floating point precision range.
///
/// @param theta is the polar angle in radian towards the z-axis.
///
/// @return the pseudorapidity towards the z-axis.
template <typename Scalar>
Scalar etaFromTheta(Scalar theta) {
return -std::log(std::tan(0.5 * theta));
if (theta <= EtaThetaConversionTraits<Scalar>::minTheta) {
return std::numeric_limits<Scalar>::infinity();
}
if (theta >= EtaThetaConversionTraits<Scalar>::maxTheta) {
return -std::numeric_limits<Scalar>::infinity();
}

return -std::log(std::tan(theta / 2));
}

/// Calculate the polar angle theta from the pseudorapidity.
///
/// This function aims to be FPE safe and returns 0/pi for eta values outside
/// the floating point precision range.
///
/// @param eta is the pseudorapidity towards the z-axis.
///
/// @return the polar angle in radian towards the z-axis.
template <typename Scalar>
Scalar thetaFromEta(Scalar eta) {
if (eta <= EtaThetaConversionTraits<Scalar>::minEta) {
return std::numbers::pi_v<Scalar>;
}
if (eta >= EtaThetaConversionTraits<Scalar>::maxEta) {
return 0;
}

return 2 * std::atan(std::exp(-eta));
}

Expand Down
73 changes: 35 additions & 38 deletions Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/TrackFinding/TrackSelector.hpp"
#include "Acts/Utilities/AngleHelpers.hpp"

#include <limits>

Expand Down Expand Up @@ -83,10 +84,6 @@ struct MockTrack {
TrackStateRange trackStatesReversed() const { return {}; }
};

double thetaFromEta(double eta) {
return 2 * std::atan(std::exp(-eta));
}

BOOST_AUTO_TEST_SUITE(TrackSelectorTests)

std::vector<double> etaValues{-5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5,
Expand All @@ -97,7 +94,7 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
TrackSelector::EtaBinnedConfig cfgBase;

MockTrack baseTrack{};
baseTrack.m_theta = thetaFromEta(eta);
baseTrack.m_theta = AngleHelpers::thetaFromEta(eta);
baseTrack.m_phi = 0.5;
baseTrack.m_pt = 0.5;
baseTrack.m_loc0 = 0.5;
Expand Down Expand Up @@ -206,9 +203,9 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).etaMin = {-1.0};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(-1.1);
track.m_theta = AngleHelpers::thetaFromEta(-1.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand All @@ -218,9 +215,9 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).etaMax = {1.0};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(1.1);
track.m_theta = AngleHelpers::thetaFromEta(1.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand All @@ -231,11 +228,11 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).etaMax = {1.0};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(-1.1);
track.m_theta = AngleHelpers::thetaFromEta(-1.1);
BOOST_CHECK(!selector.isValidTrack(track));
track.m_theta = thetaFromEta(1.1);
track.m_theta = AngleHelpers::thetaFromEta(1.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand All @@ -245,14 +242,14 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).absEtaMin = {0.2};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(-0.5);
track.m_theta = AngleHelpers::thetaFromEta(-0.5);
BOOST_CHECK(selector.isValidTrack(track));

track.m_theta = thetaFromEta(0.1);
track.m_theta = AngleHelpers::thetaFromEta(0.1);
BOOST_CHECK(!selector.isValidTrack(track));
track.m_theta = thetaFromEta(-0.1);
track.m_theta = AngleHelpers::thetaFromEta(-0.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand All @@ -262,14 +259,14 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).absEtaMax = {1.0};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(-0.5);
track.m_theta = AngleHelpers::thetaFromEta(-0.5);
BOOST_CHECK(selector.isValidTrack(track));

track.m_theta = thetaFromEta(1.1);
track.m_theta = AngleHelpers::thetaFromEta(1.1);
BOOST_CHECK(!selector.isValidTrack(track));
track.m_theta = thetaFromEta(-1.1);
track.m_theta = AngleHelpers::thetaFromEta(-1.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand All @@ -280,19 +277,19 @@ BOOST_DATA_TEST_CASE(TestSingleBinCase, bdata::make(etaValues), eta) {
cfg.cutSets.at(0).absEtaMax = {1.0};
TrackSelector selector{cfg};
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(selector.isValidTrack(track));
track.m_theta = thetaFromEta(-0.5);
track.m_theta = AngleHelpers::thetaFromEta(-0.5);
BOOST_CHECK(selector.isValidTrack(track));

track.m_theta = thetaFromEta(0.1);
track.m_theta = AngleHelpers::thetaFromEta(0.1);
BOOST_CHECK(!selector.isValidTrack(track));
track.m_theta = thetaFromEta(-0.1);
track.m_theta = AngleHelpers::thetaFromEta(-0.1);
BOOST_CHECK(!selector.isValidTrack(track));

track.m_theta = thetaFromEta(1.1);
track.m_theta = AngleHelpers::thetaFromEta(1.1);
BOOST_CHECK(!selector.isValidTrack(track));
track.m_theta = thetaFromEta(-1.1);
track.m_theta = AngleHelpers::thetaFromEta(-1.1);
BOOST_CHECK(!selector.isValidTrack(track));
}

Expand Down Expand Up @@ -373,27 +370,27 @@ BOOST_AUTO_TEST_CASE(TestSingleBinEtaCutByBinEdge) {
BOOST_TEST_INFO_SCOPE(selector.config());

MockTrack track{};
track.m_theta = thetaFromEta(0.0);
track.m_theta = AngleHelpers::thetaFromEta(0.0);
BOOST_CHECK(!selector.isValidTrack(track));

track.m_theta = thetaFromEta(0.5);
track.m_theta = AngleHelpers::thetaFromEta(0.5);
BOOST_CHECK(!selector.isValidTrack(track));

// cannot easily check on-edge behavior because of floating point arithmetic
// (it won't be exactly 1.0 in selector)
track.m_theta = thetaFromEta(1.01);
track.m_theta = AngleHelpers::thetaFromEta(1.01);
BOOST_CHECK(selector.isValidTrack(track));

track.m_theta = thetaFromEta(1.5);
track.m_theta = AngleHelpers::thetaFromEta(1.5);
BOOST_CHECK(selector.isValidTrack(track));

track.m_theta = thetaFromEta(2.0);
track.m_theta = AngleHelpers::thetaFromEta(2.0);
BOOST_CHECK(!selector.isValidTrack(track));
}

BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
MockTrack baseTrack{};
baseTrack.m_theta = thetaFromEta(1.0);
baseTrack.m_theta = AngleHelpers::thetaFromEta(1.0);
baseTrack.m_phi = 0.5;
baseTrack.m_pt = 0.5;
baseTrack.m_loc0 = 0.5;
Expand Down Expand Up @@ -425,7 +422,7 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
{
// exactly at zero
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(0.0);
track.m_theta = AngleHelpers::thetaFromEta(0.0);

BOOST_CHECK(selector.isValidTrack(track));

Expand All @@ -439,7 +436,7 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
{
// first bin
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(1.0);
track.m_theta = AngleHelpers::thetaFromEta(1.0);

BOOST_CHECK(selector.isValidTrack(track));

Expand All @@ -453,8 +450,8 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
{
// first bin edge
MockTrack track = baseTrack;
track.m_theta =
thetaFromEta(2.0 - std::numeric_limits<double>::epsilon());
track.m_theta = AngleHelpers::thetaFromEta(
2.0 - std::numeric_limits<double>::epsilon());

BOOST_CHECK(selector.isValidTrack(track));

Expand All @@ -468,7 +465,7 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
{
// second bin lower edge
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(2.0);
track.m_theta = AngleHelpers::thetaFromEta(2.0);

BOOST_CHECK(selector.isValidTrack(track));

Expand All @@ -488,7 +485,7 @@ BOOST_AUTO_TEST_CASE(TestMultiBinCuts) {
{
// second bin
MockTrack track = baseTrack;
track.m_theta = thetaFromEta(666.0);
track.m_theta = AngleHelpers::thetaFromEta(10.0);

track.*prop = -1.1;
BOOST_CHECK(selector.isValidTrack(track));
Expand Down
72 changes: 70 additions & 2 deletions Tests/UnitTests/Core/Utilities/AngleHelpersTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,24 +6,92 @@
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.

#include <boost/test/data/test_case.hpp>
#include <boost/test/unit_test.hpp>

#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
#include "Acts/Utilities/AngleHelpers.hpp"

#include <cmath>
#include <numbers>

namespace bd = boost::unit_test::data;

using Acts::AngleHelpers::etaFromTheta;
using Acts::AngleHelpers::EtaThetaConversionTraits;
using Acts::AngleHelpers::thetaFromEta;

BOOST_AUTO_TEST_SUITE(AngleHelpers)

BOOST_AUTO_TEST_CASE(EtaThetaConversion) {
CHECK_CLOSE_ABS(0.0, etaFromTheta(std::numbers::pi / 2), 1e-6);

CHECK_CLOSE_ABS(1.0, etaFromTheta(thetaFromEta(1.0)), 1e-6);

CHECK_CLOSE_ABS(1.0, thetaFromEta(etaFromTheta(1.0)), 1e-6);
}

BOOST_DATA_TEST_CASE(EtaFromThetaRobustness, bd::xrange(0, 1000, 1), exponent) {
{
// check right

float thetaRight = exponent < 30 ? std::pow(10.0f, -1.0f * exponent) : 0.0f;

float etaRight = etaFromTheta<float>(thetaRight);
BOOST_CHECK(!std::isnan(etaRight));

// check left

float thetaLeft = std::numbers::pi_v<float> - thetaRight;

float etaLeft = etaFromTheta<float>(thetaLeft);
BOOST_CHECK(!std::isnan(etaLeft));
}

{
// check right

double thetaRight = exponent < 300 ? std::pow(10.0, -1.0 * exponent) : 0.0;

double etaRight = etaFromTheta<double>(thetaRight);
BOOST_CHECK(!std::isnan(etaRight));

// check left

double thetaLeft = std::numbers::pi - thetaRight;

double etaLeft = etaFromTheta<double>(thetaLeft);
BOOST_CHECK(!std::isnan(etaLeft));
}
}

BOOST_DATA_TEST_CASE(ThetaFromEtaRobustness, bd::xrange(1.0, 1000.0, 1.0),
etaRight) {
{
// check right

float thetaRight = thetaFromEta<float>(etaRight);
BOOST_CHECK(!std::isnan(thetaRight));

// check left

float etaLeft = -etaRight;

float thetaLeft = thetaFromEta<float>(etaLeft);
BOOST_CHECK(!std::isnan(thetaLeft));
}

{
// check right

double thetaRight = thetaFromEta<double>(etaRight);
BOOST_CHECK(!std::isnan(thetaRight));

// check left

double etaLeft = -etaRight;

double thetaLeft = thetaFromEta<double>(etaLeft);
BOOST_CHECK(!std::isnan(thetaLeft));
}
}

BOOST_AUTO_TEST_SUITE_END()

0 comments on commit 37db0d4

Please sign in to comment.