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.

Sign up for GitHub

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

Jump to bottom

refactor: FPE safe eta/theta conversion #3788

Merged
merged 17 commits into from
Nov 7, 2024
Merged

refactor: FPE safe eta/theta conversion #3788

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
f6f0767
refactor: FPE safe eta/theta conversion
andiwand Oct 28, 2024
5e130b7
experiment with FPEs
andiwand Oct 28, 2024
2c22146
fix track selector tests
andiwand Oct 28, 2024
e1c3261
use std::floating_point
andiwand Oct 28, 2024
f645945
pr feedback
andiwand Oct 30, 2024
eaf6094
Merge branch 'main' into fpe-safe-eta-theta-conversion
andiwand Oct 31, 2024
2a411cd
provide clamped impl
andiwand Nov 1, 2024
a1cbdbb
test corners
andiwand Nov 1, 2024
aa33d08
Apply suggestions from code review
andiwand Nov 1, 2024
1e56b8c
pr feedback
andiwand Nov 1, 2024
3fd43bc
Merge branch 'fpe-safe-eta-theta-conversion' of github.com:andiwand/a…
andiwand Nov 1, 2024
91e82cb
nuke clamped conversion
andiwand Nov 1, 2024
e8d9d18
drop corner tests
andiwand Nov 1, 2024
5a220d7
Apply suggestions from code review
andiwand Nov 7, 2024
345e426
Update Core/include/Acts/Utilities/AngleHelpers.hpp
andiwand Nov 7, 2024
e5550eb
Update Core/include/Acts/Utilities/AngleHelpers.hpp
andiwand Nov 7, 2024
f92f144
Merge branch 'main' into fpe-safe-eta-theta-conversion
andiwand Nov 7, 2024
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
View file
Open in desktop
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> - 1e-6f;
andiwand marked this conversation as resolved.
Show resolved Hide resolved

static constexpr float minEta = -80.0f;
static constexpr float maxEta = 80.0f;
andiwand marked this conversation as resolved.
Show resolved Hide resolved
};

template <>
struct EtaThetaConversionTraits<double> {
static constexpr double minTheta = 1e-12;
static constexpr double maxTheta = std::numbers::pi - 1e-12;
andiwand marked this conversation as resolved.
Show resolved Hide resolved

static constexpr double minEta = -700.0;
static constexpr double maxEta = 700.0;
andiwand marked this conversation as resolved.
Show resolved Hide resolved
};

/// 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
View file
Open in desktop
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
View file
Open in desktop
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