Skip to content
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

refactor: FPE safe eta/theta conversion #3788

Merged
merged 17 commits into from
Nov 7, 2024
Merged
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: 43 additions & 1 deletion Core/include/Acts/Utilities/AngleHelpers.hpp
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 {};

andiwand marked this conversation as resolved.
Show resolved Hide resolved
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));
}
andiwand marked this conversation as resolved.
Show resolved Hide resolved

/// 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));
andiwand marked this conversation as resolved.
Show resolved Hide resolved
}

Expand Down
73 changes: 35 additions & 38 deletions Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp
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());
andiwand marked this conversation as resolved.
Show resolved Hide resolved

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
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));
}
}
andiwand marked this conversation as resolved.
Show resolved Hide resolved
andiwand marked this conversation as resolved.
Show resolved Hide resolved

BOOST_AUTO_TEST_SUITE_END()
Loading