feat: Add AngleHelpers.hpp to help with eta/theta conversions (acts…

…-project#3767) Adds central conversion between $\eta$ and $\theta$
AJPfleger · Oct 22, 2024 · be9777b · be9777b
1 parent 2e9eefc
commit be9777b
Show file tree

Hide file tree

Showing 14 changed files with 126 additions and 45 deletions.
diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp
@@ -12,11 +12,11 @@
 #include "Acts/EventData/TrackStateType.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 
 #include <cmath>
 #include <functional>
 #include <limits>
-#include <numeric>
 #include <ostream>
 #include <vector>
 
@@ -419,7 +419,7 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
 
   auto absEta = [&]() {
     if (_absEta == kUnset) {
-      _eta = -std::log(std::tan(theta / 2));
+      _eta = AngleHelpers::etaFromTheta(theta);
       _absEta = std::abs(_eta);
     }
     return _absEta;

diff --git a/Core/include/Acts/Utilities/AngleHelpers.hpp b/Core/include/Acts/Utilities/AngleHelpers.hpp
@@ -0,0 +1,35 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// 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/.
+
+#pragma once
+
+#include <cmath>
+
+namespace Acts::AngleHelpers {
+
+/// Calculate the pseudorapidity from the polar angle theta.
+///
+/// @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));
+}
+
+/// Calculate the polar angle theta from the pseudorapidity.
+///
+/// @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) {
+  return 2 * std::atan(std::exp(-eta));
+}
+
+}  // namespace Acts::AngleHelpers
diff --git a/Core/src/Vertexing/ImpactPointEstimator.cpp b/Core/src/Vertexing/ImpactPointEstimator.cpp
@@ -9,10 +9,10 @@
 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
 
 #include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Propagator/PropagatorOptions.hpp"
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
@@ -526,7 +526,7 @@ Result<std::pair<double, double>> ImpactPointEstimator::getLifetimeSignOfTrack(
   const double theta = params[BoundIndices::eBoundTheta];
 
   double vs = std::sin(std::atan2(direction[1], direction[0]) - phi) * d0;
-  double eta = -std::log(std::tan(theta / 2.));
+  double eta = AngleHelpers::etaFromTheta(theta);
   double dir_eta = VectorHelpers::eta(direction);
 
   double zs = (dir_eta - eta) * z0;

diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp b/Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp
@@ -11,7 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Definitions/ParticleData.hpp"
-#include "ActsExamples/EventData/SimHit.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
 
@@ -35,8 +35,8 @@ ParametricParticleGenerator::ParametricParticleGenerator(const Config& cfg)
       m_cosThetaMax(std::nextafter(std::cos(m_cfg.thetaMax),
                                    std::numeric_limits<double>::max())),
       // in case we force uniform eta generation
-      m_etaMin(-std::log(std::tan(0.5 * m_cfg.thetaMin))),
-      m_etaMax(-std::log(std::tan(0.5 * m_cfg.thetaMax))) {}
+      m_etaMin(Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMin)),
+      m_etaMax(Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMax)) {}
 
 std::pair<SimVertexContainer, SimParticleContainer>
 ParametricParticleGenerator::operator()(RandomEngine& rng) {

diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
@@ -203,7 +203,7 @@ class BranchStopper {
           } else if constexpr (std::is_same_v<
                                    T, Acts::TrackSelector::EtaBinnedConfig>) {
             double theta = trackState.parameters()[Acts::eBoundTheta];
-            double eta = -std::log(std::tan(0.5 * theta));
+            double eta = Acts::AngleHelpers::etaFromTheta(theta);
             return config.hasCuts(eta) ? &config.getCuts(eta) : nullptr;
           }
         },

diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp
@@ -40,14 +40,6 @@ struct LoopHook : public Acts::ExaTrkXHook {
   }
 };
 
-// TODO do we have these function in the repo somewhere?
-float theta(float r, float z) {
-  return std::atan2(r, z);
-}
-float eta(float r, float z) {
-  return -std::log(std::tan(0.5 * theta(r, z)));
-}
-
 }  // namespace
 
 ActsExamples::TrackFindingAlgorithmExaTrkX::TrackFindingAlgorithmExaTrkX(
@@ -187,7 +179,7 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithmExaTrkX::execute(
         break; case NF::eZ:           f[ift] = sp.z();
         break; case NF::eX:           f[ift] = sp.x();
         break; case NF::eY:           f[ift] = sp.y();
-        break; case NF::eEta:         f[ift] = eta(std::hypot(sp.x(), sp.y()), sp.z());
+        break; case NF::eEta:         f[ift] = Acts::VectorHelpers::eta(Acts::Vector3{sp.x(), sp.y(), sp.z()});
         break; case NF::eClusterX:    f[ift] = cl1->sizeLoc0;
         break; case NF::eClusterY:    f[ift] = cl1->sizeLoc1;
         break; case NF::eCellSum:     f[ift] = cl1->sumActivations();
@@ -198,8 +190,8 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithmExaTrkX::execute(
         break; case NF::eCluster2Phi: f[ift] = std::atan2(cl2->globalPosition[Acts::ePos1], cl2->globalPosition[Acts::ePos0]);
         break; case NF::eCluster1Z:   f[ift] = cl1->globalPosition[Acts::ePos2];
         break; case NF::eCluster2Z:   f[ift] = cl2->globalPosition[Acts::ePos2];
-        break; case NF::eCluster1Eta: f[ift] = eta(std::hypot(cl1->globalPosition[Acts::ePos0], cl1->globalPosition[Acts::ePos1]), cl1->globalPosition[Acts::ePos2]);
-        break; case NF::eCluster2Eta: f[ift] = eta(std::hypot(cl2->globalPosition[Acts::ePos0], cl2->globalPosition[Acts::ePos1]), cl2->globalPosition[Acts::ePos2]);
+        break; case NF::eCluster1Eta: f[ift] = Acts::VectorHelpers::eta(Acts::Vector3{cl1->globalPosition[Acts::ePos0], cl1->globalPosition[Acts::ePos1], cl1->globalPosition[Acts::ePos2]});
+        break; case NF::eCluster2Eta: f[ift] = Acts::VectorHelpers::eta(Acts::Vector3{cl2->globalPosition[Acts::ePos0], cl2->globalPosition[Acts::ePos1], cl2->globalPosition[Acts::ePos2]});
       }
       // clang-format on
 

diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSelector.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSelector.cpp
@@ -9,6 +9,7 @@
 #include "ActsExamples/TruthTracking/TrackParameterSelector.hpp"
 
 #include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 #include "ActsExamples/EventData/Track.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 
@@ -41,7 +42,7 @@ ActsExamples::ProcessCode ActsExamples::TrackParameterSelector::execute(
   };
   auto isValidTrack = [&](const auto& trk) {
     const auto theta = trk.template get<Acts::eBoundTheta>();
-    const auto eta = -std::log(std::tan(theta / 2));
+    const auto eta = Acts::AngleHelpers::etaFromTheta(theta);
     // define charge selection
     return within(trk.transverseMomentum(), m_cfg.ptMin, m_cfg.ptMax) &&
            within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) &&

diff --git a/Examples/Python/src/Generators.cpp b/Examples/Python/src/Generators.cpp
@@ -7,25 +7,21 @@
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
 #include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Definitions/PdgParticle.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/EventData/SimParticle.hpp"
-#include "ActsExamples/Framework/RandomNumbers.hpp"
 #include "ActsExamples/Generators/EventGenerator.hpp"
 #include "ActsExamples/Generators/MultiplicityGenerators.hpp"
 #include "ActsExamples/Generators/ParametricParticleGenerator.hpp"
 #include "ActsExamples/Generators/VertexGenerators.hpp"
 
-#include <array>
 #include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <memory>
-#include <optional>
 #include <string>
 #include <utility>
-#include <vector>
 
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
@@ -36,16 +32,6 @@ class IReader;
 
 namespace py = pybind11;
 
-namespace {
-double thetaToEta(double theta) {
-  assert(theta != 0);
-  return -1 * std::log(std::tan(theta / 2.));
-}
-double etaToTheta(double eta) {
-  return 2 * std::atan(std::exp(-eta));
-}
-}  // namespace
-
 namespace Acts::Python {
 
 void addGenerators(Context& ctx) {
@@ -218,12 +204,12 @@ void addGenerators(Context& ctx) {
         .def_property(
             "eta",
             [](Config& cfg) {
-              return std::pair{thetaToEta(cfg.thetaMin),
-                               thetaToEta(cfg.thetaMax)};
+              return std::pair{Acts::AngleHelpers::etaFromTheta(cfg.thetaMin),
+                               Acts::AngleHelpers::etaFromTheta(cfg.thetaMax)};
             },
             [](Config& cfg, std::pair<double, double> value) {
-              cfg.thetaMin = etaToTheta(value.first);
-              cfg.thetaMax = etaToTheta(value.second);
+              cfg.thetaMin = Acts::AngleHelpers::thetaFromEta(value.first);
+              cfg.thetaMax = Acts::AngleHelpers::thetaFromEta(value.second);
             });
   }
 

diff --git a/Examples/Scripts/fullMaterial.C b/Examples/Scripts/fullMaterial.C
@@ -13,6 +13,8 @@
  *      Author: jhrdinka
  */
 
+#include "Acts/Utilities/AngleHelpers.hpp"
+
 #include "TFile.h"
 #include "TH1F.h"
 #include "TH2F.h"
@@ -111,7 +113,7 @@ fullMaterial(std::string inFile,
   for (auto& mrecord : mrecords) {
     std::vector<Acts::MaterialStep> steps = mrecord.materialSteps();
     float                           theta = mrecord.theta();
-    float                           eta   = -log(tan(theta * 0.5));
+    float                           eta   = Acts::AngleHelpers::etaFromTheta(theta);
     float                           phi   = VectorHelpers::phi(mrecord);
 
     float thickness = 0.;

diff --git a/Examples/Scripts/momentumDistributions.C b/Examples/Scripts/momentumDistributions.C
@@ -6,6 +6,8 @@
 // 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 "Acts/Utilities/AngleHelpers.hpp"
+
 #include "TFile.h"
 #include "TH1F.h"
 #include "TH2F.h"
@@ -106,7 +108,7 @@ void momentumDistributions(std::string inFile, std::string treeName,
 
     for (int j = 0; j < x->size(); j++) {
       float hitTheta = std::atan2(std::hypot(x->at(j), y->at(j)), z->at(j));
-      hitsEta->Fill(-log(tan(hitTheta * 0.5)));
+      hitsEta->Fill(Acts::AngleHelpers::etaFromTheta(hitTheta));
       hitsTheta->Fill(hitTheta);
       hitsZ->Fill(z->at(j));
     }

diff --git a/Plugins/Hashing/include/Acts/Plugins/Hashing/HashingTraining.ipp b/Plugins/Hashing/include/Acts/Plugins/Hashing/HashingTraining.ipp
@@ -8,7 +8,9 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 
+#include <cmath>
 #include <cstdint>
 
 #include <annoy/annoylib.h>
@@ -49,7 +51,7 @@ AnnoyModel HashingTrainingAlgorithm<SpacePointContainer>::execute(
     Scalar y = spacePoint->y() / Acts::UnitConstants::mm;
 
     // Helix transform
-    Scalar phi = atan2(y, x);
+    Scalar phi = std::atan2(y, x);
 
     std::vector<double> vec(f);
     // Avoid potential null pointer dereference
@@ -59,9 +61,9 @@ AnnoyModel HashingTrainingAlgorithm<SpacePointContainer>::execute(
     if (f >= 2) {
       Scalar z = spacePoint->z() / Acts::UnitConstants::mm;
       Scalar r2 = x * x + y * y;
-      Scalar rho = sqrt(r2 + z * z);
-      Scalar theta = acos(z / rho);
-      Scalar eta = -log(tan(0.5 * theta));
+      Scalar rho = std::sqrt(r2 + z * z);
+      Scalar theta = std::acos(z / rho);
+      Scalar eta = Acts::AngleHelpers::etaFromTheta(theta);
       vec[1] = eta;
     }
 

diff --git a/Tests/UnitTests/Core/Utilities/AngleHelpersTests.cpp b/Tests/UnitTests/Core/Utilities/AngleHelpersTests.cpp
@@ -0,0 +1,29 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// 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/unit_test.hpp>
+
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
+
+#include <numbers>
+
+using Acts::AngleHelpers::etaFromTheta;
+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_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Utilities/CMakeLists.txt b/Tests/UnitTests/Core/Utilities/CMakeLists.txt
@@ -56,6 +56,8 @@ add_unittest(ParticleData ParticleDataTests.cpp)
 add_unittest(Zip ZipTests.cpp)
 add_unittest(TransformRange TransformRangeTests.cpp)
 add_unittest(MathHelpers MathHelpersTests.cpp)
+add_unittest(AngleHelpers AngleHelpersTests.cpp)
+add_unittest(VectorHelpers VectorHelpersTests.cpp)
 
 add_unittest(TrackHelpers TrackHelpersTests.cpp)
 add_unittest(GraphViz GraphVizTests.cpp)
diff --git a/Tests/UnitTests/Core/Utilities/VectorHelpersTests.cpp b/Tests/UnitTests/Core/Utilities/VectorHelpersTests.cpp
@@ -0,0 +1,30 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// 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/unit_test.hpp>
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Utilities/VectorHelpers.hpp"
+
+#include <numbers>
+
+using Acts::VectorHelpers::eta;
+using Acts::VectorHelpers::theta;
+
+BOOST_AUTO_TEST_SUITE(AngleHelpers)
+
+BOOST_AUTO_TEST_CASE(EtaFromVector) {
+  CHECK_CLOSE_ABS(0.0, eta(Acts::Vector3{1, 0, 0}), 1e-6);
+}
+
+BOOST_AUTO_TEST_CASE(ThetaFromVector) {
+  CHECK_CLOSE_ABS(std::numbers::pi / 2, theta(Acts::Vector3{1, 0, 0}), 1e-6);
+}
+
+BOOST_AUTO_TEST_SUITE_END()