diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml index 5686d842b66..edd31d051a4 100644 --- a/.github/workflows/checks.yml +++ b/.github/workflows/checks.yml @@ -116,6 +116,16 @@ jobs: - name: Check run: > CI/check_spelling + math_macros: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - name: Check + run: > + CI/check_math_macros.py . --exclude "thirdparty/*" missing_includes: runs-on: ubuntu-latest steps: diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4c6a7b166fc..b0a213e607d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -37,8 +37,6 @@ clang_tidy: -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_C_COMPILER=clang -DPython_EXECUTABLE=$(which python3) - -DACTS_RUN_CLANG_TIDY=ON - -DACTS_BUILD_ODD=OFF # Main clang-tidy run during cmake compilation - CI/clang_tidy/run_clang_tidy.sh clang-tidy build @@ -168,7 +166,7 @@ build_exatrkx: # - git clone $CLONE_URL src # - cd src # - git checkout $HEAD_SHA -# - pip3 install -r Examples/Python/tests/requirements_ubuntu2004.txt +# - pip3 install -r Examples/Python/tests/requirements.txt # - nvidia-smi # - pytest -rFsv -k test_exatrkx diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 76af4847b7f..68c78341c5d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -93,3 +93,11 @@ repos: name: Leftover conflict markers language: system entry: git diff --staged --check + + - repo: local + hooks: + - id: math_macros + name: math_macros + language: system + entry: CI/check_math_macros.py + files: \.(cpp|hpp|ipp|cu|cuh)$ diff --git a/CI/check_math_macros.py b/CI/check_math_macros.py new file mode 100755 index 00000000000..29da233372e --- /dev/null +++ b/CI/check_math_macros.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 + +from pathlib import Path +import os +import argparse +from fnmatch import fnmatch +import re +import sys + + +math_constants = [ + ("M_PI", "std::numbers::pi"), + ("M_PI_2", "std::numbers::pi / 2."), + ("M_PI_4", "std::numbers::pi / 4."), + ("M_1_PI", "std::numbers::inv_pi"), + ("M_2_PI", "2. * std::numbers::inv_pi"), + ("M_2_SQRTPI", "2. * std::numbers::inv_sqrtpi"), + ("M_E", "std::numbers::e"), + ("M_LOG2E", "std::numbers::log2e"), + ("M_LOG10E", "std::numbers::log10e"), + ("M_LN2", "std::numbers::ln2"), + ("M_LN10", "std::numbers::ln10"), + ("M_SQRT2", "std::numbers::sqrt2"), + ("M_SQRT1_2", "1. / std::numbers::sqrt2"), + ("M_SQRT3", "std::numbers::sqrt3"), + ("M_INV_SQRT3", "std::numbers::inv_sqrt3"), + ("M_EGAMMA", "std::numbers::egamma"), + ("M_PHI", "std::numbers::phi"), +] + + +github = "GITHUB_ACTIONS" in os.environ + + +def handle_file( + file: Path, fix: bool, math_const: tuple[str, str] +) -> list[tuple[int, str]]: + ex = re.compile(rf"(? 0: + changed_lines.append((i, oline)) + + if fix and len(changed_lines) > 0: + file.write_text("\n".join(lines) + "\n") + + return changed_lines + + +def main(): + p = argparse.ArgumentParser() + p.add_argument("input", nargs="+") + p.add_argument("--fix", action="store_true", help="Attempt to fix M_* macros.") + p.add_argument("--exclude", "-e", action="append", default=[]) + + args = p.parse_args() + + exit_code = 0 + + inputs = [] + + if len(args.input) == 1 and os.path.isdir(args.input[0]): + # walk over all files + for root, _, files in os.walk(args.input[0]): + root = Path(root) + for filename in files: + # get the full path of the file + filepath = root / filename + if filepath.suffix not in ( + ".hpp", + ".cpp", + ".ipp", + ".h", + ".C", + ".c", + ".cu", + ".cuh", + ): + continue + + if any([fnmatch(str(filepath), e) for e in args.exclude]): + continue + + inputs.append(filepath) + else: + for file in args.input: + inputs.append(Path(file)) + + for filepath in inputs: + for math_const in math_constants: + changed_lines = handle_file( + file=filepath, fix=args.fix, math_const=math_const + ) + if len(changed_lines) > 0: + exit_code = 1 + print() + print(filepath) + for i, oline in changed_lines: + print(f"{i}: {oline}") + + if github: + print( + f"::error file={filepath},line={i+1},title=Do not use macro {math_const[0]}::Replace {math_const[0]} with std::{math_const[1]}" + ) + + if exit_code == 1 and github: + print(f"::info You will need in each flagged file #include ") + + return exit_code + + +if "__main__" == __name__: + sys.exit(main()) diff --git a/CI/physmon/config/pythia8_ttbar.yml b/CI/physmon/config/pythia8_ttbar.yml index 301aba1c7e5..6a71fe03988 100644 --- a/CI/physmon/config/pythia8_ttbar.yml +++ b/CI/physmon/config/pythia8_ttbar.yml @@ -65,3 +65,8 @@ exclude: - particle - generation - sub_particle + - e_loss + - total_x0 + - total_l0 + - number_of_hits + - outcome diff --git a/CI/physmon/workflows/physmon_simulation.py b/CI/physmon/workflows/physmon_simulation.py index 37874663187..206890540b8 100755 --- a/CI/physmon/workflows/physmon_simulation.py +++ b/CI/physmon/workflows/physmon_simulation.py @@ -81,8 +81,7 @@ preSelectParticles=None, postSelectParticles=ParticleSelectorConfig(removeSecondaries=True), inputParticles="particles_input", - outputParticlesInitial="particles_initial_fatras", - outputParticlesFinal="particles_final_fatras", + outputParticles="particles_fatras", outputSimHits="simhits_fatras", outputDirRoot=tp / "fatras", ) @@ -99,8 +98,7 @@ killAfterTime=25 * u.ns, killSecondaries=True, inputParticles="particles_input", - outputParticlesInitial="particles_initial_geant4", - outputParticlesFinal="particles_final_geant4", + outputParticles="particles_geant4", outputSimHits="simhits_geant4", outputDirRoot=tp / "geant4", ) diff --git a/CMakeLists.txt b/CMakeLists.txt index f5209e0d210..601264badc0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -290,8 +290,12 @@ if(ACTS_SETUP_BOOST) endif() if(Boost_VERSION VERSION_EQUAL "1.85.0") + set(_boost_version_severity WARNING) + if(ACTS_BUILD_EXAMPLES) + set(_boost_version_severity FATAL_ERROR) + endif() message( - WARNING + ${_boost_version_severity} "Boost 1.85.0 is known to be broken (https://github.com/boostorg/container/issues/273). Please use a different version." ) endif() diff --git a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp index c5c0293b8ad..fbccb611d1a 100644 --- a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp +++ b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -103,8 +104,8 @@ class ScoreBasedAmbiguityResolution { double pTMin = 0 * UnitConstants::GeV; double pTMax = 1e5 * UnitConstants::GeV; - double phiMin = -M_PI * UnitConstants::rad; - double phiMax = M_PI * UnitConstants::rad; + double phiMin = -std::numbers::pi * UnitConstants::rad; + double phiMax = std::numbers::pi * UnitConstants::rad; double etaMin = -5; double etaMax = 5; diff --git a/Core/include/Acts/Definitions/Units.hpp b/Core/include/Acts/Definitions/Units.hpp index c41aa7dd830..65d7cee775d 100644 --- a/Core/include/Acts/Definitions/Units.hpp +++ b/Core/include/Acts/Definitions/Units.hpp @@ -8,6 +8,8 @@ #pragma once +#include + namespace Acts { /// @verbatim embed:rst:leading-slashes @@ -170,7 +172,7 @@ constexpr double h = 3600.0 * s; // Angles, native unit radian constexpr double mrad = 1e-3; constexpr double rad = 1.0; -constexpr double degree = 0.017453292519943295; // = M_PI / 180.0 * rad; +constexpr double degree = std::numbers::pi / 180. / rad; // Energy/mass/momentum, native unit GeV constexpr double GeV = 1.0; constexpr double eV = 1e-9 * GeV; diff --git a/Core/include/Acts/EventData/TransformationHelpers.hpp b/Core/include/Acts/EventData/TransformationHelpers.hpp index 39240bf469e..e31e29fd3c0 100644 --- a/Core/include/Acts/EventData/TransformationHelpers.hpp +++ b/Core/include/Acts/EventData/TransformationHelpers.hpp @@ -15,6 +15,8 @@ #include "Acts/Utilities/Result.hpp" #include "Acts/Utilities/detail/periodic.hpp" +#include + namespace Acts { class Surface; @@ -25,8 +27,9 @@ class Surface; /// @return Reflected bound track parameters vector inline BoundVector reflectBoundParameters(const BoundVector& boundParams) { BoundVector reflected = boundParams; - auto [phi, theta] = detail::normalizePhiTheta( - boundParams[eBoundPhi] - M_PI, M_PI - boundParams[eBoundTheta]); + auto [phi, theta] = + detail::normalizePhiTheta(boundParams[eBoundPhi] - std::numbers::pi, + std::numbers::pi - boundParams[eBoundTheta]); reflected[eBoundPhi] = phi; reflected[eBoundTheta] = theta; reflected[eBoundQOverP] = -boundParams[eBoundQOverP]; diff --git a/Core/include/Acts/EventData/detail/CalculateResiduals.hpp b/Core/include/Acts/EventData/detail/CalculateResiduals.hpp deleted file mode 100644 index 33d800d16e6..00000000000 --- a/Core/include/Acts/EventData/detail/CalculateResiduals.hpp +++ /dev/null @@ -1,96 +0,0 @@ -// 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 "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Utilities/detail/periodic.hpp" - -#include - -#include - -namespace Acts::detail { - -/// Residuals between bound reference parameters and a measured subspace. -/// -/// @tparam index_container_t SequenceContainer for measured indices -/// @tparam measured_t Measured parameters vector type -/// @tparam residuals_t Residuals vector type -/// @param[in] size Size of the measured parameters subspace -/// @param[in] indices Indices of measured subspace, must have `size` entries -/// @param[in] reference Reference bound parameters -/// @param[in] measured Measured parameters subspace -/// @param[out] residuals Resulting residuals in the measured subspace -/// -/// @note The separate `size` parameter is also used to allow the selection of -/// the correct residuals methods depending on the parameters type. -template -inline void calculateResiduals(BoundIndices size, - const index_container_t& indices, - const BoundVector& reference, - const Eigen::MatrixBase& measured, - Eigen::MatrixBase& residuals) { - using OutputScalar = typename residuals_t::Scalar; - - EIGEN_STATIC_ASSERT_VECTOR_ONLY(measured_t); - EIGEN_STATIC_ASSERT_VECTOR_ONLY(residuals_t); - assert((size <= eBoundSize) && "Measured subspace is too large"); - assert((size <= measured.size()) && "Inconsistent measured size"); - assert((size <= residuals.size()) && "Inconsistent residuals size"); - - for (std::size_t i = 0; i < size; ++i) { - std::size_t fullIndex = indices[i]; - // this is neither elegant nor smart but it is the simplest solution. - // - // only phi must be handled specially here. the theta limits can only be - // correctly handled if phi is updated, too. since we can not ensure that - // both are available, it is probably less error-prone to treat theta as a - // regular, unrestricted parameter. - if (fullIndex == eBoundPhi) { - residuals[i] = difference_periodic( - measured[i], reference[fullIndex], - static_cast(2 * M_PI)); - } else { - residuals[i] = measured[i] - reference[fullIndex]; - } - } -} - -/// Residuals between free reference parameters and a measured subspace. -/// -/// @tparam index_container_t SequenceContainer for measured inidices -/// @tparam measured_t Measured parameters vector type -/// @tparam residuals_t Residuals vector type -/// @param[in] size Size of the measured parameters subspace -/// @param[in] indices Indices of measured subspace, must have `size` entries -/// @param[in] reference Reference free parameters -/// @param[in] measured Measured parameters subspace -/// @param[out] residuals Resulting residuals in the measured subspace -/// -/// @note The separate `size` parameter is also used to allow the selection of -/// the correct residuals methods depending on the parameters type. -template -inline void calculateResiduals(FreeIndices size, - const index_container_t& indices, - const FreeVector& reference, - const Eigen::MatrixBase& measured, - Eigen::MatrixBase& residuals) { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(measured_t); - EIGEN_STATIC_ASSERT_VECTOR_ONLY(residuals_t); - assert((size <= eFreeSize) && "Measured subspace is too large"); - assert((size <= measured.size()) && "Inconsistent measured size"); - assert((size <= residuals.size()) && "Inconsistent residuals size"); - - for (std::size_t i = 0; i < size; ++i) { - // all free parameters are unrestricted. no need to call parameter traits - residuals[i] = measured[i] - reference[indices[i]]; - } -} - -} // namespace Acts::detail diff --git a/Core/include/Acts/EventData/detail/ParameterTraits.hpp b/Core/include/Acts/EventData/detail/ParameterTraits.hpp index ae5eb8eeb30..ff3cdd81835 100644 --- a/Core/include/Acts/EventData/detail/ParameterTraits.hpp +++ b/Core/include/Acts/EventData/detail/ParameterTraits.hpp @@ -13,6 +13,7 @@ #include #include #include +#include namespace Acts::detail { @@ -106,12 +107,12 @@ struct CyclicParameterTraits { // // The functions names are chosen to be consistent w/ std::numeric_limits struct PhiBoundParameterLimits { - static constexpr double lowest() { return -M_PI; } - static constexpr double max() { return M_PI; } + static constexpr double lowest() { return -std::numbers::pi; } + static constexpr double max() { return std::numbers::pi; } }; struct ThetaBoundParameterLimits { static constexpr double lowest() { return 0; } - static constexpr double max() { return M_PI; } + static constexpr double max() { return std::numbers::pi; } }; // Traits implementation structs for single parameters. diff --git a/Core/include/Acts/EventData/detail/TestSourceLink.hpp b/Core/include/Acts/EventData/detail/TestSourceLink.hpp index 34a234bbf4d..affec84aafa 100644 --- a/Core/include/Acts/EventData/detail/TestSourceLink.hpp +++ b/Core/include/Acts/EventData/detail/TestSourceLink.hpp @@ -139,14 +139,13 @@ void testSourceLinkCalibratorReturn( trackState.allocateCalibrated(2); trackState.template calibrated<2>() = sl.parameters; trackState.template calibratedCovariance<2>() = sl.covariance; - trackState.template setSubspaceIndices( - std::array{sl.indices[0], sl.indices[1]}); + trackState.setSubspaceIndices(std::array{sl.indices[0], sl.indices[1]}); } else if (sl.indices[0] != Acts::eBoundSize) { trackState.allocateCalibrated(1); trackState.template calibrated<1>() = sl.parameters.head<1>(); trackState.template calibratedCovariance<1>() = sl.covariance.topLeftCorner<1, 1>(); - trackState.template setSubspaceIndices(std::array{sl.indices[0]}); + trackState.setSubspaceIndices(std::array{sl.indices[0]}); } else { throw std::runtime_error( "Tried to extract measurement from invalid TestSourceLink"); diff --git a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp index 002af04cd59..7f368442f40 100644 --- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp +++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -105,8 +106,9 @@ class CylinderVolumeBounds : public VolumeBounds { /// @param bevelMinZ The bevel angle, in radians, for the negative side /// @param bevelMaxZ The bevel angle, in radians, for the positive side CylinderVolumeBounds(ActsScalar rmin, ActsScalar rmax, ActsScalar halfz, - ActsScalar halfphi = M_PI, ActsScalar avgphi = 0., - ActsScalar bevelMinZ = 0., ActsScalar bevelMaxZ = 0.); + ActsScalar halfphi = std::numbers::pi_v, + ActsScalar avgphi = 0., ActsScalar bevelMinZ = 0., + ActsScalar bevelMaxZ = 0.); /// Constructor - from a fixed size array /// diff --git a/Core/include/Acts/Geometry/PortalShell.hpp b/Core/include/Acts/Geometry/PortalShell.hpp index 7405a64ad0a..fe9f606c576 100644 --- a/Core/include/Acts/Geometry/PortalShell.hpp +++ b/Core/include/Acts/Geometry/PortalShell.hpp @@ -34,12 +34,10 @@ class PortalShellBase { /// Virtusl destructor virtual ~PortalShellBase() = default; - /// Connect a volume to the outer side of all portal shells. Which "side" is - /// "outer" depends on the volume type. - /// This method essentially creates a @c TrivialPortalLink on the unconnected - /// side of each portal that is part of the chell + /// Fill the open slots of the shell with a @c TrivialPortalLink + /// to the given @p volume. /// @param volume The volume to connect - virtual void connectOuter(TrackingVolume& volume) = 0; + virtual void fill(TrackingVolume& volume) = 0; /// Get the number of portals in the shell. This number depends on the volume /// type @@ -88,8 +86,8 @@ class CylinderPortalShell : public PortalShellBase { /// @param face The face to set the portal virtual void setPortal(std::shared_ptr portal, Face face) = 0; - /// @copydoc PortalShellBase::connectOuter - void connectOuter(TrackingVolume& volume) override; + /// @copydoc PortalShellBase::fill + void fill(TrackingVolume& volume) override; }; /// Output stream operator for the CylinderPortalShell::Face enum diff --git a/Core/include/Acts/Geometry/SurfaceArrayCreator.hpp b/Core/include/Acts/Geometry/SurfaceArrayCreator.hpp index bc52d700c72..584e8ee6992 100644 --- a/Core/include/Acts/Geometry/SurfaceArrayCreator.hpp +++ b/Core/include/Acts/Geometry/SurfaceArrayCreator.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -257,7 +258,7 @@ class SurfaceArrayCreator { // ...so by injecting them into atan2, we get the angle between them auto dPhi = std::atan2(sin_dPhi_n2, cos_dPhi_n2); - return std::abs(dPhi) < M_PI / 180.; + return std::abs(dPhi) < std::numbers::pi / 180.; } if (bValue == Acts::BinningValue::binZ) { diff --git a/Core/include/Acts/Geometry/VolumeBounds.hpp b/Core/include/Acts/Geometry/VolumeBounds.hpp index 3b5810ff4c8..28bcc029166 100644 --- a/Core/include/Acts/Geometry/VolumeBounds.hpp +++ b/Core/include/Acts/Geometry/VolumeBounds.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -33,12 +34,14 @@ struct OrientedSurface { // Planar definitions to help construct the boundary surfaces static const Transform3 s_planeXY = Transform3::Identity(); -static const Transform3 s_planeYZ = AngleAxis3(0.5 * M_PI, Vector3::UnitY()) * - AngleAxis3(0.5 * M_PI, Vector3::UnitZ()) * - Transform3::Identity(); -static const Transform3 s_planeZX = AngleAxis3(-0.5 * M_PI, Vector3::UnitX()) * - AngleAxis3(-0.5 * M_PI, Vector3::UnitZ()) * - Transform3::Identity(); +static const Transform3 s_planeYZ = + AngleAxis3(std::numbers::pi / 2., Vector3::UnitY()) * + AngleAxis3(std::numbers::pi / 2., Vector3::UnitZ()) * + Transform3::Identity(); +static const Transform3 s_planeZX = + AngleAxis3(-std::numbers::pi / 2., Vector3::UnitX()) * + AngleAxis3(-std::numbers::pi / 2., Vector3::UnitZ()) * + Transform3::Identity(); /// @class VolumeBounds /// diff --git a/Core/include/Acts/Propagator/RiddersPropagator.ipp b/Core/include/Acts/Propagator/RiddersPropagator.ipp index cef7c79bd4f..1b1764a13d7 100644 --- a/Core/include/Acts/Propagator/RiddersPropagator.ipp +++ b/Core/include/Acts/Propagator/RiddersPropagator.ipp @@ -8,6 +8,8 @@ #include "Acts/Definitions/TrackParametrization.hpp" +#include + template template auto Acts::RiddersPropagator::propagate( @@ -150,8 +152,8 @@ bool Acts::RiddersPropagator::inconsistentDerivativesOnDisc( for (unsigned int j = 0; j < derivatives.size(); j++) { // If there is at least one with a similar angle then it seems to work // properly - if (i != j && - std::abs(derivatives[i](1) - derivatives[j](1)) < 0.5 * M_PI) { + if (i != j && std::abs(derivatives[i](1) - derivatives[j](1)) < + std::numbers::pi / 2.) { jumpedAngle = false; break; } @@ -179,8 +181,8 @@ Acts::RiddersPropagator::wiggleParameter( // Treatment for theta if (param == eBoundTheta) { const double current_theta = start.template get(); - if (current_theta + h > M_PI) { - h = M_PI - current_theta; + if (current_theta + h > std::numbers::pi) { + h = std::numbers::pi - current_theta; } if (current_theta + h < 0) { h = -current_theta; @@ -202,10 +204,14 @@ Acts::RiddersPropagator::wiggleParameter( if (param == eBoundPhi) { double phi0 = nominal(Acts::eBoundPhi); double phi1 = r.endParameters->parameters()(Acts::eBoundPhi); - if (std::abs(phi1 + 2. * M_PI - phi0) < std::abs(phi1 - phi0)) { - derivatives.back()[Acts::eBoundPhi] = (phi1 + 2. * M_PI - phi0) / h; - } else if (std::abs(phi1 - 2. * M_PI - phi0) < std::abs(phi1 - phi0)) { - derivatives.back()[Acts::eBoundPhi] = (phi1 - 2. * M_PI - phi0) / h; + if (std::abs(phi1 + 2. * std::numbers::pi - phi0) < + std::abs(phi1 - phi0)) { + derivatives.back()[Acts::eBoundPhi] = + (phi1 + 2. * std::numbers::pi - phi0) / h; + } else if (std::abs(phi1 - 2. * std::numbers::pi - phi0) < + std::abs(phi1 - phi0)) { + derivatives.back()[Acts::eBoundPhi] = + (phi1 - 2. * std::numbers::pi - phi0) / h; } } } diff --git a/Core/include/Acts/Propagator/detail/LoopProtection.hpp b/Core/include/Acts/Propagator/detail/LoopProtection.hpp index a5efe15306e..62c7e363423 100644 --- a/Core/include/Acts/Propagator/detail/LoopProtection.hpp +++ b/Core/include/Acts/Propagator/detail/LoopProtection.hpp @@ -11,6 +11,8 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Utilities/Logger.hpp" +#include + namespace Acts::detail { /// Estimate the loop protection limit @@ -46,7 +48,8 @@ void setupLoopProtection(propagator_state_t& state, const stepper_t& stepper, // Transverse component at start is taken for the loop protection const double p = stepper.absoluteMomentum(state.stepping); // Calculate the full helix path - const double helixPath = state.options.direction * 2 * M_PI * p / B; + const double helixPath = + state.options.direction * 2 * std::numbers::pi * p / B; // And set it as the loop limit if it overwrites the internal limit const double loopLimit = state.options.loopFraction * helixPath; const double previousLimit = pathAborter.internalLimit; diff --git a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp index f3850e102ce..f4f56946c84 100644 --- a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp +++ b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp @@ -22,105 +22,6 @@ #include namespace Acts { -/// @todo: -/// 1) Implement the simple Line and Circle fit based on Taubin Circle fit -/// 2) Implement the simple Line and Parabola fit (from HPS reconstruction by -/// Robert Johnson) - -/// Estimate the track parameters on the xy plane from at least three space -/// points. It assumes the trajectory projection on the xy plane is a circle, -/// i.e. the magnetic field is along global z-axis. -/// -/// The method is based on V. Karimaki NIM A305 (1991) 187-191: -/// https://doi.org/10.1016/0168-9002(91)90533-V -/// - no weights are used in Karimaki's fit; d0 is the distance of the point of -/// closest approach to the origin, 1/R is the curvature, phi is the angle of -/// the direction propagation (counter clockwise as positive) at the point of -/// cloest approach. -/// -/// @tparam spacepoint_iterator_t The type of space point iterator -/// -/// @param spBegin is the begin iterator for the space points -/// @param spEnd is the end iterator for the space points -/// @param logger A logger instance -/// -/// @return optional bound track parameters with the estimated d0, phi and 1/R -/// stored with the indices, eBoundLoc0, eBoundPhi and eBoundQOverP, -/// respectively. The bound parameters with other indices are set to zero. -template -std::optional estimateTrackParamsFromSeed( - spacepoint_iterator_t spBegin, spacepoint_iterator_t spEnd, - const Logger& logger = getDummyLogger()) { - // Check the number of provided space points - std::size_t numSP = std::distance(spBegin, spEnd); - if (numSP < 3) { - ACTS_ERROR("At least three space points are required."); - return std::nullopt; - } - - ActsScalar x2m = 0., xm = 0.; - ActsScalar xym = 0.; - ActsScalar y2m = 0., ym = 0.; - ActsScalar r2m = 0., r4m = 0.; - ActsScalar xr2m = 0., yr2m = 0.; - - for (spacepoint_iterator_t it = spBegin; it != spEnd; it++) { - if (*it == nullptr) { - ACTS_ERROR("Empty space point found. This should not happen."); - return std::nullopt; - } - - const auto& sp = *it; - - ActsScalar x = sp->x(); - ActsScalar y = sp->y(); - ActsScalar r2 = x * x + y * y; - x2m += x * x; - xm += x; - xym += x * y; - y2m += y * y; - ym += y; - r2m += r2; - r4m += r2 * r2; - xr2m += x * r2; - yr2m += y * r2; - numSP++; - } - x2m = x2m / numSP; - xm = xm / numSP; - xym = xym / numSP; - y2m = y2m / numSP; - ym = ym / numSP; - r2m = r2m / numSP; - r4m = r4m / numSP; - xr2m = xr2m / numSP; - yr2m = yr2m / numSP; - - ActsScalar Cxx = x2m - xm * xm; - ActsScalar Cxy = xym - xm * ym; - ActsScalar Cyy = y2m - ym * ym; - ActsScalar Cxr2 = xr2m - xm * r2m; - ActsScalar Cyr2 = yr2m - ym * r2m; - ActsScalar Cr2r2 = r4m - r2m * r2m; - - ActsScalar q1 = Cr2r2 * Cxy - Cxr2 * Cyr2; - ActsScalar q2 = Cr2r2 * (Cxx - Cyy) - Cxr2 * Cxr2 + Cyr2 * Cyr2; - - ActsScalar phi = 0.5 * std::atan(2 * q1 / q2); - ActsScalar k = (std::sin(phi) * Cxr2 - std::cos(phi) * Cyr2) * (1. / Cr2r2); - ActsScalar delta = -k * r2m + std::sin(phi) * xm - std::cos(phi) * ym; - - ActsScalar rho = (2 * k) / (std::sqrt(1 - 4 * delta * k)); - ActsScalar d0 = (2 * delta) / (1 + std::sqrt(1 - 4 * delta * k)); - - // Initialize the bound parameters vector - BoundVector params = BoundVector::Zero(); - params[eBoundLoc0] = d0; - params[eBoundPhi] = phi; - params[eBoundQOverP] = rho; - - return params; -} /// Estimate the full track parameters from three space points /// diff --git a/Core/include/Acts/Seeding/GbtsDataStorage.hpp b/Core/include/Acts/Seeding/GbtsDataStorage.hpp index 64ba0bff1b9..c28b0cc2413 100644 --- a/Core/include/Acts/Seeding/GbtsDataStorage.hpp +++ b/Core/include/Acts/Seeding/GbtsDataStorage.hpp @@ -14,6 +14,7 @@ #include #include +#include #include namespace Acts { @@ -117,12 +118,12 @@ class GbtsEtaBin { // float phi = pN->m_sp.phi(); // float phi = (std::atan(pN->m_sp.x() / pN->m_sp.y())); float phi = pN->m_spGbts.phi(); - if (phi <= M_PI - dphi) { + if (phi <= std::numbers::pi_v - dphi) { continue; } - m_vPhiNodes.push_back( - std::pair(phi - 2 * M_PI, nIdx)); + m_vPhiNodes.push_back(std::pair( + phi - static_cast(2. * std::numbers::pi), nIdx)); } for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) { @@ -134,11 +135,11 @@ class GbtsEtaBin { for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) { GbtsNode *pN = m_vn.at(nIdx); float phi = pN->m_spGbts.phi(); - if (phi >= -M_PI + dphi) { + if (phi >= -std::numbers::pi_v + dphi) { break; } - m_vPhiNodes.push_back( - std::pair(phi + 2 * M_PI, nIdx)); + m_vPhiNodes.push_back(std::pair( + phi + static_cast(2. * std::numbers::pi), nIdx)); } } diff --git a/Core/include/Acts/Seeding/SeedFinderConfig.hpp b/Core/include/Acts/Seeding/SeedFinderConfig.hpp index bb20f9fead0..93dacaeafd8 100644 --- a/Core/include/Acts/Seeding/SeedFinderConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFinderConfig.hpp @@ -15,6 +15,7 @@ #include #include +#include #include namespace Acts { @@ -33,8 +34,8 @@ struct SeedFinderConfig { /// Geometry Settings + Detector ROI /// (r, z, phi) range for limiting location of all measurements and grid /// creation - float phiMin = -M_PI; - float phiMax = M_PI; + float phiMin = -std::numbers::pi_v; + float phiMax = std::numbers::pi_v; float zMin = -2800 * Acts::UnitConstants::mm; float zMax = 2800 * Acts::UnitConstants::mm; float rMax = 600 * Acts::UnitConstants::mm; diff --git a/Core/include/Acts/Seeding/SeedFinderGbts.ipp b/Core/include/Acts/Seeding/SeedFinderGbts.ipp index ab557d86977..ed2fe157acf 100644 --- a/Core/include/Acts/Seeding/SeedFinderGbts.ipp +++ b/Core/include/Acts/Seeding/SeedFinderGbts.ipp @@ -9,7 +9,6 @@ // SeedFinderGbts.ipp // TODO: update to C++17 style -#include "Acts/Definitions/Algebra.hpp" //for M_PI #include "Acts/Geometry/Extent.hpp" #include "Acts/Seeding/SeedFilter.hpp" #include "Acts/Seeding/SeedFinder.hpp" @@ -23,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -54,7 +54,7 @@ void SeedFinderGbts::loadSpacePoints( m_storage->addSpacePoint(gbtssp, (m_config.m_useClusterWidth > 0)); } - m_config.m_phiSliceWidth = 2 * M_PI / m_config.m_nMaxPhiSlice; + m_config.m_phiSliceWidth = 2 * std::numbers::pi / m_config.m_nMaxPhiSlice; m_storage->sortByPhi(); @@ -390,10 +390,10 @@ void SeedFinderGbts::runGbts_TrackFinder( float dPhi = pNS->m_p[3] - Phi1; - if (dPhi < -M_PI) { - dPhi += 2 * M_PI; - } else if (dPhi > M_PI) { - dPhi -= 2 * M_PI; + if (dPhi < -std::numbers::pi_v) { + dPhi += static_cast(2 * std::numbers::pi); + } else if (dPhi > std::numbers::pi_v) { + dPhi -= static_cast(2 * std::numbers::pi); } if (dPhi < -m_config.cut_dphi_max || dPhi > m_config.cut_dphi_max) { diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp index 594e35b44cc..a4129fad556 100644 --- a/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp @@ -65,7 +65,7 @@ class SeedFinderOrthogonal { */ ~SeedFinderOrthogonal() = default; - SeedFinderOrthogonal() = delete; + SeedFinderOrthogonal() = default; SeedFinderOrthogonal(const SeedFinderOrthogonal &) = delete; SeedFinderOrthogonal &operator=( @@ -254,7 +254,7 @@ class SeedFinderOrthogonal { /** * @brief The logger */ - std::unique_ptr m_logger{nullptr}; + std::unique_ptr m_logger{getDummyLogger().clone()}; }; } // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp index 4bb710abc78..ac0b51f5d34 100644 --- a/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp @@ -14,6 +14,7 @@ #include "Acts/Utilities/Delegate.hpp" #include +#include namespace Acts { // forward declaration to avoid cyclic dependence @@ -28,8 +29,8 @@ struct SeedFinderOrthogonalConfig { /// Seeding parameters for geometry settings and detector ROI // Limiting location of all measurements - float phiMin = -M_PI; - float phiMax = M_PI; + float phiMin = -std::numbers::pi_v; + float phiMax = std::numbers::pi_v; /// limiting location of measurements float zMin = -2800 * Acts::UnitConstants::mm; float zMax = 2800 * Acts::UnitConstants::mm; diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp index a2a214fa39b..dcbe2289ffb 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp @@ -71,7 +71,7 @@ struct CylindricalSpacePointGridConfig { // maximum phi value for phiAxis construction float phiMax = std::numbers::pi_v; // Multiplicator for the number of phi-bins. The minimum number of phi-bins - // depends on min_pt, magnetic field: 2*M_PI/(minPT particle phi-deflection). + // depends on min_pt, magnetic field: 2*pi/(minPT particle phi-deflection). // phiBinDeflectionCoverage is a multiplier for this number. If // numPhiNeighbors (in the configuration of the BinFinders) is configured to // return 1 neighbor on either side of the current phi-bin (and you want to diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp index b99cdbdfaf9..2a8f3dad41d 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp @@ -7,6 +7,7 @@ // file, You can obtain one at https://mozilla.org/MPL/2.0/. #include +#include template Acts::CylindricalSpacePointGrid @@ -91,16 +92,14 @@ Acts::CylindricalSpacePointGridCreator::createGrid( // divide 2pi by angle delta to get number of phi-bins // size is always 2pi even for regions of interest - phiBins = static_cast(std::ceil(2 * M_PI / deltaPhi)); + phiBins = static_cast(std::ceil(2 * std::numbers::pi / deltaPhi)); // need to scale the number of phi bins accordingly to the number of // consecutive phi bins in the seed making step. // Each individual bin should be approximately a fraction (depending on this // number) of the maximum expected azimutal deflection. // set protection for large number of bins, by default it is large - if (phiBins > config.maxPhiBins) { - phiBins = config.maxPhiBins; - } + phiBins = std::min(phiBins, config.maxPhiBins); } Acts::Axis phiAxis( diff --git a/Core/include/Acts/Surfaces/AnnulusBounds.hpp b/Core/include/Acts/Surfaces/AnnulusBounds.hpp index f1283ec63b8..942c20dc764 100644 --- a/Core/include/Acts/Surfaces/AnnulusBounds.hpp +++ b/Core/include/Acts/Surfaces/AnnulusBounds.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -231,7 +232,7 @@ inline double AnnulusBounds::phiMax() const { } inline bool AnnulusBounds::coversFullAzimuth() const { - return (std::abs((get(eMinPhiRel) - get(eMaxPhiRel)) - M_PI) < + return (std::abs((get(eMinPhiRel) - get(eMaxPhiRel)) - std::numbers::pi) < s_onSurfaceTolerance); } diff --git a/Core/include/Acts/Surfaces/ConeBounds.hpp b/Core/include/Acts/Surfaces/ConeBounds.hpp index 55b2203686c..d0d684016a8 100644 --- a/Core/include/Acts/Surfaces/ConeBounds.hpp +++ b/Core/include/Acts/Surfaces/ConeBounds.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -55,7 +56,7 @@ class ConeBounds : public SurfaceBounds { /// @param halfphi is the half opening angle (default is pi) /// @param avphi is the phi value around which the bounds are opened /// (default=0) - ConeBounds(double alpha, bool symm, double halfphi = M_PI, + ConeBounds(double alpha, bool symm, double halfphi = std::numbers::pi, double avphi = 0.) noexcept(false); /// Constructor - open cone with alpha, minz and maxz, by @@ -67,7 +68,8 @@ class ConeBounds : public SurfaceBounds { /// @param halfphi is the half opening angle (default is pi) /// @param avphi is the phi value around which the bounds are opened /// (default=0) - ConeBounds(double alpha, double minz, double maxz, double halfphi = M_PI, + ConeBounds(double alpha, double minz, double maxz, + double halfphi = std::numbers::pi, double avphi = 0.) noexcept(false); /// Constructor - from parameters array @@ -141,14 +143,14 @@ inline std::vector ConeBounds::values() const { } inline void ConeBounds::checkConsistency() noexcept(false) { - if (get(eAlpha) < 0. || get(eAlpha) >= M_PI) { + if (get(eAlpha) < 0. || get(eAlpha) >= std::numbers::pi) { throw std::invalid_argument("ConeBounds: invalid open angle."); } if (get(eMinZ) > get(eMaxZ) || std::abs(get(eMinZ) - get(eMaxZ)) < s_epsilon) { throw std::invalid_argument("ConeBounds: invalid z range setup."); } - if (get(eHalfPhiSector) < 0. || abs(eHalfPhiSector) > M_PI) { + if (get(eHalfPhiSector) < 0. || abs(eHalfPhiSector) > std::numbers::pi) { throw std::invalid_argument("ConeBounds: invalid phi sector setup."); } if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) { diff --git a/Core/include/Acts/Surfaces/ConeSurface.hpp b/Core/include/Acts/Surfaces/ConeSurface.hpp index a2b2572bba8..2f88bfe9532 100644 --- a/Core/include/Acts/Surfaces/ConeSurface.hpp +++ b/Core/include/Acts/Surfaces/ConeSurface.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include namespace Acts { @@ -61,7 +62,7 @@ class ConeSurface : public RegularSurface { /// @param zmax is the z range over which the cone spans /// @param halfPhi is the opening angle for cone ssectors ConeSurface(const Transform3& transform, double alpha, double zmin, - double zmax, double halfPhi = M_PI); + double zmax, double halfPhi = std::numbers::pi); /// Constructor from HepTransform and ConeBounds /// diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp index 51508055e06..d144c8b995b 100644 --- a/Core/include/Acts/Surfaces/CylinderBounds.hpp +++ b/Core/include/Acts/Surfaces/CylinderBounds.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -69,11 +70,11 @@ class CylinderBounds : public SurfaceBounds { /// @param avgPhi (optional) The phi value from which the opening angle spans /// @param bevelMinZ (optional) The bevel on the negative z side /// @param bevelMaxZ (optional) The bevel on the positive z sid The bevel on the positive z side - CylinderBounds(double r, double halfZ, double halfPhi = M_PI, + CylinderBounds(double r, double halfZ, double halfPhi = std::numbers::pi, double avgPhi = 0., double bevelMinZ = 0., double bevelMaxZ = 0.) noexcept(false) : m_values({r, halfZ, halfPhi, avgPhi, bevelMinZ, bevelMaxZ}), - m_closed(std::abs(halfPhi - M_PI) < s_epsilon) { + m_closed(std::abs(halfPhi - std::numbers::pi) < s_epsilon) { checkConsistency(); } @@ -82,7 +83,8 @@ class CylinderBounds : public SurfaceBounds { /// @param values The parameter values CylinderBounds(const std::array& values) noexcept(false) : m_values(values), - m_closed(std::abs(values[eHalfPhiSector] - M_PI) < s_epsilon) { + m_closed(std::abs(values[eHalfPhiSector] - std::numbers::pi) < + s_epsilon) { checkConsistency(); } diff --git a/Core/include/Acts/Surfaces/CylinderSurface.hpp b/Core/include/Acts/Surfaces/CylinderSurface.hpp index 09b953afce1..53c52f13f50 100644 --- a/Core/include/Acts/Surfaces/CylinderSurface.hpp +++ b/Core/include/Acts/Surfaces/CylinderSurface.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include namespace Acts { @@ -57,7 +58,7 @@ class CylinderSurface : public RegularSurface { /// @param bevelMinZ (optional) The bevel on the negative z side /// @param bevelMaxZ (optional) The bevel on the positive z sid The bevel on the positive z side CylinderSurface(const Transform3& transform, double radius, double halfz, - double halfphi = M_PI, double avphi = 0., + double halfphi = std::numbers::pi, double avphi = 0., double bevelMinZ = 0., double bevelMaxZ = 0.); /// Constructor from Transform3 and CylinderBounds arguments diff --git a/Core/include/Acts/Surfaces/DiscSurface.hpp b/Core/include/Acts/Surfaces/DiscSurface.hpp index e2c10df4023..d2db14fc566 100644 --- a/Core/include/Acts/Surfaces/DiscSurface.hpp +++ b/Core/include/Acts/Surfaces/DiscSurface.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include namespace Acts { @@ -64,7 +65,7 @@ class DiscSurface : public RegularSurface { /// @param hphisec The opening angle of the disc surface and is optional /// the default is a full disc DiscSurface(const Transform3& transform, double rmin, double rmax, - double hphisec = M_PI); + double hphisec = std::numbers::pi); /// Constructor for Discs from Transform3, \f$ r_{min}, r_{max}, hx_{min}, /// hx_{max} \f$ diff --git a/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp index a3e4e5ea28f..32c63c2aef6 100644 --- a/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp +++ b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -52,7 +53,7 @@ class DiscTrapezoidBounds : public DiscBounds { /// @param avgPhi average phi value /// @param stereo optional stero angle applied DiscTrapezoidBounds(double halfXminR, double halfXmaxR, double minR, - double maxR, double avgPhi = M_PI_2, + double maxR, double avgPhi = std::numbers::pi / 2., double stereo = 0.) noexcept(false); /// Constructor - from fixed size array diff --git a/Core/include/Acts/Surfaces/EllipseBounds.hpp b/Core/include/Acts/Surfaces/EllipseBounds.hpp index aec08d10381..ea90e3d2a88 100644 --- a/Core/include/Acts/Surfaces/EllipseBounds.hpp +++ b/Core/include/Acts/Surfaces/EllipseBounds.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -57,7 +58,8 @@ class EllipseBounds : public PlanarBounds { /// @param halfPhi spanning phi sector (is set to pi as default) /// @param averagePhi average phi (is set to 0. as default) EllipseBounds(double innerRx, double innerRy, double outerRx, double outerRy, - double halfPhi = M_PI, double averagePhi = 0.) noexcept(false) + double halfPhi = std::numbers::pi, + double averagePhi = 0.) noexcept(false) : m_values({innerRx, innerRy, outerRx, outerRy, halfPhi, averagePhi}), m_boundingBox(m_values[eInnerRy], m_values[eOuterRy]) { checkConsistency(); @@ -134,7 +136,7 @@ inline void EllipseBounds::checkConsistency() noexcept(false) { get(eOuterRy) <= 0.) { throw std::invalid_argument("EllipseBounds: invalid along y axis."); } - if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > M_PI) { + if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > std::numbers::pi) { throw std::invalid_argument("EllipseBounds: invalid phi sector setup."); } if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) { diff --git a/Core/include/Acts/Surfaces/RadialBounds.hpp b/Core/include/Acts/Surfaces/RadialBounds.hpp index 0c7aea81400..ef3e8996572 100644 --- a/Core/include/Acts/Surfaces/RadialBounds.hpp +++ b/Core/include/Acts/Surfaces/RadialBounds.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -47,7 +48,7 @@ class RadialBounds : public DiscBounds { /// @param maxR The outer radius /// @param halfPhi The half opening angle (Pi for full angular coverage) /// @param avgPhi The average phi for the disc/ring sector - RadialBounds(double minR, double maxR, double halfPhi = M_PI, + RadialBounds(double minR, double maxR, double halfPhi = std::numbers::pi, double avgPhi = 0.) noexcept(false) : m_values({minR, maxR, halfPhi, avgPhi}) { checkConsistency(); @@ -142,7 +143,7 @@ inline double RadialBounds::rMax() const { } inline bool RadialBounds::coversFullAzimuth() const { - return (get(eHalfPhiSector) == M_PI); + return (get(eHalfPhiSector) == std::numbers::pi); } inline bool RadialBounds::insideRadialBounds(double R, double tolerance) const { @@ -167,7 +168,7 @@ inline void RadialBounds::checkConsistency() noexcept(false) { if (get(eMinR) < 0. || get(eMaxR) <= 0. || get(eMinR) > get(eMaxR)) { throw std::invalid_argument("RadialBounds: invalid radial setup"); } - if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > M_PI) { + if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > std::numbers::pi) { throw std::invalid_argument("RadialBounds: invalid phi sector setup."); } if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) { diff --git a/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp b/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp index dd18f3e800e..4384c17139c 100644 --- a/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp +++ b/Core/include/Acts/Surfaces/detail/VerticesHelper.hpp @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -29,10 +30,11 @@ namespace Acts::detail::VerticesHelper { /// @param quarterSegments number of segments used to approximate a segment quarter /// /// @return a vector of generated phi values -std::vector phiSegments(ActsScalar phiMin = -M_PI, - ActsScalar phiMax = M_PI, - const std::vector& phiRefs = {}, - unsigned int quarterSegments = 2u); +std::vector phiSegments( + ActsScalar phiMin = -std::numbers::pi_v, + ActsScalar phiMax = std::numbers::pi_v, + const std::vector& phiRefs = {}, + unsigned int quarterSegments = 2u); /// Helper method to create a regular 2 or 3 D segment /// between two phi values with a given number of segments @@ -83,11 +85,11 @@ std::vector segmentVertices( /// @param quarterSegments number of segments used to approximate a segment quarter /// /// @return a vector of 2d-vectors -std::vector ellipsoidVertices(ActsScalar innerRx, ActsScalar innerRy, - ActsScalar outerRx, ActsScalar outerRy, - ActsScalar avgPhi = 0., - ActsScalar halfPhi = M_PI, - unsigned int quarterSegments = 2u); +std::vector ellipsoidVertices( + ActsScalar innerRx, ActsScalar innerRy, ActsScalar outerRx, + ActsScalar outerRy, ActsScalar avgPhi = 0., + ActsScalar halfPhi = std::numbers::pi_v, + unsigned int quarterSegments = 2u); /// Construct vertices on an disc/wheel-like bound object. /// @@ -98,10 +100,10 @@ std::vector ellipsoidVertices(ActsScalar innerRx, ActsScalar innerRy, /// @param quarterSegments number of segments used to approximate a segment quarter /// /// @return a vector of 2d-vectors -std::vector circularVertices(ActsScalar innerR, ActsScalar outerR, - ActsScalar avgPhi = 0., - ActsScalar halfPhi = M_PI, - unsigned int quarterSegments = 2u); +std::vector circularVertices( + ActsScalar innerR, ActsScalar outerR, ActsScalar avgPhi = 0., + ActsScalar halfPhi = std::numbers::pi_v, + unsigned int quarterSegments = 2u); /// Check if the point is inside the polygon w/o any tolerances. /// diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp index 459ab996857..5a2cf53df27 100644 --- a/Core/include/Acts/TrackFinding/TrackSelector.hpp +++ b/Core/include/Acts/TrackFinding/TrackSelector.hpp @@ -304,24 +304,32 @@ inline TrackSelector::Config& TrackSelector::Config::pt(double min, inline std::ostream& operator<<(std::ostream& os, const TrackSelector::Config& cuts) { - auto print = [&](const char* name, const auto& min, const auto& max) { + // for printing cuts set up with `within` + auto printMinMax = [&](const char* name, const auto& min, const auto& max) { os << " - " << min << " <= " << name << " < " << max << "\n"; }; + // for printing cuts set up with `checkMin` + auto printMin = [&](const char* name, const auto& min) { + os << " - " << min << " <= " << name << "\n"; + }; + // for printing cuts set up with `checkMax` + auto printMax = [&](const char* name, const auto& max) { + os << " - " << name << " <= " << max << "\n"; + }; - print("loc0", cuts.loc0Min, cuts.loc0Max); - print("loc1", cuts.loc1Min, cuts.loc1Max); - print("time", cuts.timeMin, cuts.timeMax); - print("phi", cuts.phiMin, cuts.phiMax); - print("eta", cuts.etaMin, cuts.etaMax); - print("absEta", cuts.absEtaMin, cuts.absEtaMax); - print("pt", cuts.ptMin, cuts.ptMax); - print("nHoles", 0, cuts.maxHoles); - print("nOutliers", 0, cuts.maxOutliers); - print("nHoles + nOutliers", 0, cuts.maxHolesAndOutliers); - print("nSharedHits", 0, cuts.maxSharedHits); - print("chi2", 0.0, cuts.maxChi2); - os << " - " << cuts.minMeasurements << " <= nMeasurements\n"; - + printMinMax("loc0", cuts.loc0Min, cuts.loc0Max); + printMinMax("loc1", cuts.loc1Min, cuts.loc1Max); + printMinMax("time", cuts.timeMin, cuts.timeMax); + printMinMax("phi", cuts.phiMin, cuts.phiMax); + printMinMax("eta", cuts.etaMin, cuts.etaMax); + printMinMax("absEta", cuts.absEtaMin, cuts.absEtaMax); + printMinMax("pt", cuts.ptMin, cuts.ptMax); + printMax("nHoles", cuts.maxHoles); + printMax("nOutliers", cuts.maxOutliers); + printMax("nHoles + nOutliers", cuts.maxHolesAndOutliers); + printMax("nSharedHits", cuts.maxSharedHits); + printMax("chi2", cuts.maxChi2); + printMin("nMeasurements", cuts.minMeasurements); return os; } diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 666ebf6aa23..5d0c5601bea 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -86,7 +87,7 @@ struct BetheHeitlerApproxSingleCmp { ret[0].weight = 1.0; - const double c = x / std::log(2); + const double c = x / std::numbers::ln2; ret[0].mean = std::pow(2, -c); ret[0].var = std::pow(3, -c) - std::pow(4, -c); diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index 968843fa479..8cd67c38c2b 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -283,6 +283,84 @@ struct ScatteringProperties { bool m_materialIsValid; }; +/// @brief A container to manage all properties of a gx2f system +/// +/// This struct manages the mathematical infrastructure for the gx2f. It +/// initializes and maintains the extended aMatrix and extended bVector. +struct Gx2fSystem { + public: + /// @brief Constructor to initialize matrices and vectors to zero based on specified dimensions. + /// + /// @param nDims Number of dimensions for the extended matrix and vector. + explicit Gx2fSystem(std::size_t nDims) + : m_nDims{nDims}, + m_aMatrix{Eigen::MatrixXd::Zero(nDims, nDims)}, + m_bVector{Eigen::VectorXd::Zero(nDims)} {} + + // Accessor for nDims (const reference). + std::size_t nDims() const { return m_nDims; } + + // Accessor for chi2 + double chi2() const { return m_chi2; } + + // Modifier for chi2 + double& chi2() { return m_chi2; } + + // Accessor for the matrix. + const Eigen::MatrixXd& aMatrix() const { return m_aMatrix; } + + // Accessor for a modifiable reference to the matrix. + Eigen::MatrixXd& aMatrix() { return m_aMatrix; } + + // Accessor for the vector. + const Eigen::VectorXd& bVector() const { return m_bVector; } + + // Accessor for a modifiable reference to the vector. + Eigen::VectorXd& bVector() { return m_bVector; } + + // Accessor for NDF + std::size_t ndf() const { return m_ndf; } + + // Modifier for NDF + std::size_t& ndf() { return m_ndf; } + + // It automatically deduces if we want to fit e.g. q/p and adjusts itself + // later. We have only 3 cases, because we always have l0, l1, phi, theta: + // - 4: no magnetic field -> q/p is empty + // - 5: no time measurement -> time is not fittable + // - 6: full fit + std::size_t findRequiredNdf() { + std::size_t ndfSystem = 0; + if (m_aMatrix(4, 4) == 0) { + ndfSystem = 4; + } else if (m_aMatrix(5, 5) == 0) { + ndfSystem = 5; + } else { + ndfSystem = 6; + } + + return ndfSystem; + } + + bool isWellDefined() { return m_ndf > findRequiredNdf(); } + + private: + /// Number of dimensions of the (extended) system + std::size_t m_nDims; + + /// Sum of chi-squared values. + double m_chi2 = 0.; + + /// Extended matrix for accumulation. + Eigen::MatrixXd m_aMatrix; + + /// Extended vector for accumulation. + Eigen::VectorXd m_bVector; + + /// Number of degrees of freedom of the system + std::size_t m_ndf = 0u; +}; + /// @brief Process measurements and fill the aMatrix and bVector /// /// The function processes each measurement for the GX2F Actor fitting process. @@ -292,15 +370,12 @@ struct ScatteringProperties { /// @tparam kMeasDim Number of dimensions of the measurement /// @tparam track_state_t The type of the track state /// -/// @param aMatrixExtended The aMatrix sums over the second derivatives -/// @param bVectorExtended The bVector sums over the first derivatives -/// @param chi2sum The total chi2 of the system +/// @param extendedSystem All parameters of the current equation system /// @param jacobianFromStart The Jacobian matrix from start to the current state /// @param trackState The track state to analyse /// @param logger A logger instance template -void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, - Eigen::VectorXd& bVectorExtended, double& chi2sum, +void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const std::vector& jacobianFromStart, const track_state_t& trackState, const Logger& logger) { @@ -317,13 +392,10 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, } // Create an extended Jacobian. This one contains only eBoundSize rows, - // because the rest is irrelevant + // because the rest is irrelevant. We fill it in the next steps. // TODO make dimsExtendedParams template with unrolling - const std::size_t dimsExtendedParams = aMatrixExtended.rows(); - - // We create an empty Jacobian and fill it in the next steps Eigen::MatrixXd extendedJacobian = - Eigen::MatrixXd::Zero(eBoundSize, dimsExtendedParams); + Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); // This part of the Jacobian comes from the material-less propagation extendedJacobian.topLeftCorner() = @@ -360,13 +432,14 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, const ActsVector residual = measurement - projPredicted; // Finally contribute to chi2sum, aMatrix, and bVector - chi2sum += (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); + extendedSystem.chi2() += + (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); - aMatrixExtended += + extendedSystem.aMatrix() += (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) .eval(); - bVectorExtended += + extendedSystem.bVector() += (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) .eval() .transpose(); @@ -410,17 +483,14 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, /// /// @tparam track_state_t The type of the track state /// -/// @param aMatrixExtended The aMatrix sums over the second derivatives -/// @param bVectorExtended The bVector sums over the first derivatives -/// @param chi2sum The total chi2 of the system +/// @param extendedSystem All parameters of the current equation system /// @param nMaterialsHandled How many materials we already handled. Used for the offset. /// @param scatteringMap The scattering map, containing all scattering angles and covariances /// @param trackState The track state to analyse /// @param logger A logger instance template void addMaterialToGx2fSums( - Eigen::MatrixXd& aMatrixExtended, Eigen::VectorXd& bVectorExtended, - double& chi2sum, const std::size_t nMaterialsHandled, + Gx2fSystem& extendedSystem, const std::size_t nMaterialsHandled, const std::unordered_map& scatteringMap, const track_state_t& trackState, const Logger& logger) { @@ -444,18 +514,18 @@ void addMaterialToGx2fSums( const ActsScalar invCov = scatteringMapId->second.invCovarianceMaterial(); // Phi contribution - aMatrixExtended(deltaPosition, deltaPosition) += + extendedSystem.aMatrix()(deltaPosition, deltaPosition) += invCov * sinThetaLoc * sinThetaLoc; - bVectorExtended(deltaPosition, 0) -= + extendedSystem.bVector()(deltaPosition, 0) -= invCov * scatteringAngles[eBoundPhi] * sinThetaLoc; - chi2sum += invCov * scatteringAngles[eBoundPhi] * sinThetaLoc * - scatteringAngles[eBoundPhi] * sinThetaLoc; + extendedSystem.chi2() += invCov * scatteringAngles[eBoundPhi] * sinThetaLoc * + scatteringAngles[eBoundPhi] * sinThetaLoc; // Theta Contribution - aMatrixExtended(deltaPosition + 1, deltaPosition + 1) += invCov; - bVectorExtended(deltaPosition + 1, 0) -= + extendedSystem.aMatrix()(deltaPosition + 1, deltaPosition + 1) += invCov; + extendedSystem.bVector()(deltaPosition + 1, 0) -= invCov * scatteringAngles[eBoundTheta]; - chi2sum += + extendedSystem.chi2() += invCov * scatteringAngles[eBoundTheta] * scatteringAngles[eBoundTheta]; ACTS_VERBOSE( @@ -494,19 +564,16 @@ void addMaterialToGx2fSums( /// /// @tparam track_proxy_t The type of the track proxy /// -/// @param track A mutable track proxy to operate on -/// @param aMatrixExtended The aMatrix, summing over second derivatives for the fitting system -/// @param bVectorExtended The bVector, summing over first derivatives for the fitting system -/// @param chi2sum Accumulated chi2 value of the system -/// @param countNdf The number of degrees of freedom counted so far +/// @param track A constant track proxy to inspect +/// @param extendedSystem All parameters of the current equation system /// @param multipleScattering Flag to consider multiple scattering in the calculation -/// @param scatteringMap Map of geometry identifiers to scattering properties, containing all scattering angles and covariances +/// @param scatteringMap Map of geometry identifiers to scattering properties, +/// containing scattering angles and validation status /// @param geoIdVector A vector to store geometry identifiers for tracking processed elements /// @param logger A logger instance template void fillGx2fSystem( - const track_proxy_t track, Eigen::MatrixXd& aMatrixExtended, - Eigen::VectorXd& bVectorExtended, double& chi2sum, std::size_t& countNdf, + const track_proxy_t track, Gx2fSystem& extendedSystem, const bool multipleScattering, const std::unordered_map& scatteringMap, @@ -559,11 +626,11 @@ void fillGx2fSystem( "Found measurement with less than 1 or more than 6 dimension(s)."); } - countNdf += measDim; + extendedSystem.ndf() += measDim; visit_measurement(measDim, [&](auto N) { - addMeasurementToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, - jacobianFromStart, trackState, logger); + addMeasurementToGx2fSums(extendedSystem, jacobianFromStart, + trackState, logger); }); } @@ -574,15 +641,59 @@ void fillGx2fSystem( jacobianFromStart.emplace_back(BoundMatrix::Identity()); // Add the material contribution to the system - addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, - geoIdVector.size(), scatteringMap, trackState, - logger); + addMaterialToGx2fSums(extendedSystem, geoIdVector.size(), scatteringMap, + trackState, logger); geoIdVector.emplace_back(geoId); } } } +/// @brief Count the valid material states in a track for scattering calculations. +/// +/// This function counts the valid material surfaces encountered in a track +/// by examining each track state. The count is based on the presence of +/// material flags and the availability of scattering information for each +/// surface. +/// +/// @tparam track_proxy_t The type of the track proxy +/// +/// @param track A constant track proxy to inspect +/// @param scatteringMap Map of geometry identifiers to scattering properties, +/// containing scattering angles and validation status +/// @param logger A logger instance +template +std::size_t countMaterialStates( + const track_proxy_t track, + const std::unordered_map& + scatteringMap, + const Logger& logger) { + std::size_t nMaterialSurfaces = 0; + ACTS_DEBUG("Count the valid material surfaces."); + for (const auto& trackState : track.trackStates()) { + const auto typeFlags = trackState.typeFlags(); + const bool stateHasMaterial = typeFlags.test(TrackStateFlag::MaterialFlag); + + if (!stateHasMaterial) { + continue; + } + + // Get and store geoId for the current material surface + const GeometryIdentifier geoId = trackState.referenceSurface().geometryId(); + + const auto scatteringMapId = scatteringMap.find(geoId); + assert(scatteringMapId != scatteringMap.end() && + "No scattering angles found for material surface."); + if (!scatteringMapId->second.materialIsValid()) { + continue; + } + + nMaterialSurfaces++; + } + + return nMaterialSurfaces; +} + /// @brief Calculate and update the covariance of the fitted parameters /// /// This function calculates the covariance of the fitted parameters using @@ -592,13 +703,11 @@ void fillGx2fSystem( /// no qop/time fit) /// /// @param fullCovariancePredicted The covariance matrix to update -/// @param aMatrixExtended The matrix containing the coefficients of the linear system. -/// @param ndfSystem The number of degrees of freedom, determining the size of meaning full block +/// @param extendedSystem All parameters of the current equation system /// /// @return deltaParams The calculated delta parameters. void updateGx2fCovarianceParams(BoundMatrix& fullCovariancePredicted, - Eigen::MatrixXd& aMatrixExtended, - const std::size_t ndfSystem); + Gx2fSystem& extendedSystem); /// Global Chi Square fitter (GX2F) implementation. /// @@ -1128,8 +1237,6 @@ class Gx2Fitter { BoundVector deltaParams = BoundVector::Zero(); double chi2sum = 0; double oldChi2sum = std::numeric_limits::max(); - BoundMatrix aMatrix = BoundMatrix::Zero(); - BoundVector bVector = BoundVector::Zero(); // We need to create a temporary track container. We create several times a // new track and delete it after updating the parameters. However, if we @@ -1145,10 +1252,6 @@ class Gx2Fitter { // and used for the final track std::size_t tipIndex = Acts::MultiTrajectoryTraits::kInvalid; - // Here we will store, the ndf of the system. It automatically deduces if we - // want to fit e.g. q/p and adjusts itself later. - std::size_t ndfSystem = std::numeric_limits::max(); - // The scatteringMap stores for each visited surface their scattering // properties std::unordered_map scatteringMap; @@ -1203,12 +1306,12 @@ class Gx2Fitter { // existing states, but this needs some more thinking. trackContainerTemp.clear(); - auto propagationResult = m_propagator.template propagate(propagatorState); + auto propagationResult = m_propagator.propagate(propagatorState); // Run the fitter - auto result = m_propagator.template makeResult(std::move(propagatorState), - propagationResult, - propagatorOptions, false); + auto result = + m_propagator.makeResult(std::move(propagatorState), propagationResult, + propagatorOptions, false); if (!result.ok()) { ACTS_ERROR("Propagation failed: " << result.error()); @@ -1242,48 +1345,20 @@ class Gx2Fitter { track.tipIndex() = tipIndex; track.linkForward(); - // This goes up for each measurement (for each dimension) - std::size_t countNdf = 0; - // Count the material surfaces, to set up the system. In the multiple // scattering case, we need to extend our system. - std::size_t nMaterialSurfaces = 0; - if (multipleScattering) { - ACTS_DEBUG("Count the valid material surfaces."); - for (const auto& trackState : track.trackStates()) { - const auto typeFlags = trackState.typeFlags(); - const bool stateHasMaterial = - typeFlags.test(TrackStateFlag::MaterialFlag); - - if (!stateHasMaterial) { - continue; - } - - // Get and store geoId for the current material surface - const GeometryIdentifier geoId = - trackState.referenceSurface().geometryId(); - - const auto scatteringMapId = scatteringMap.find(geoId); - assert(scatteringMapId != scatteringMap.end() && - "No scattering angles found for material surface."); - if (!scatteringMapId->second.materialIsValid()) { - continue; - } - - nMaterialSurfaces++; - } - } + const std::size_t nMaterialSurfaces = + multipleScattering + ? countMaterialStates(track, scatteringMap, *m_addToSumLogger) + : 0u; // We need 6 dimensions for the bound parameters and 2 * nMaterialSurfaces // dimensions for the scattering angles. const std::size_t dimsExtendedParams = eBoundSize + 2 * nMaterialSurfaces; - // Set to zero before filling - chi2sum = 0; - Eigen::MatrixXd aMatrixExtended = - Eigen::MatrixXd::Zero(dimsExtendedParams, dimsExtendedParams); - Eigen::VectorXd bVectorExtended = - Eigen::VectorXd::Zero(dimsExtendedParams); + // System that we fill with the information gathered by the actor and + // evaluate later + Gx2fSystem extendedSystem{dimsExtendedParams}; // This vector stores the IDs for each visited material in order. We use // it later for updating the scattering angles. We cannot use @@ -1291,22 +1366,10 @@ class Gx2Fitter { // all stored material in each propagation. std::vector geoIdVector; - fillGx2fSystem(track, aMatrixExtended, bVectorExtended, chi2sum, countNdf, - multipleScattering, scatteringMap, geoIdVector, - *m_addToSumLogger); - - // Get required number of degrees of freedom ndfSystem. - // We have only 3 cases, because we always have l0, l1, phi, theta - // 4: no magnetic field -> q/p is empty - // 5: no time measurement -> time not fittable - // 6: full fit - if (aMatrixExtended(4, 4) == 0) { - ndfSystem = 4; - } else if (aMatrixExtended(5, 5) == 0) { - ndfSystem = 5; - } else { - ndfSystem = 6; - } + fillGx2fSystem(track, extendedSystem, multipleScattering, scatteringMap, + geoIdVector, *m_addToSumLogger); + + chi2sum = extendedSystem.chi2(); // This check takes into account the evaluated dimensions of the // measurements. To fit, we need at least NDF+1 measurements. However, we @@ -1317,50 +1380,45 @@ class Gx2Fitter { // We skip the check during the first iteration, since we cannot guarantee // to hit all/enough measurement surfaces with the initial parameter // guess. - if ((nUpdate > 0) && (ndfSystem + 1 > countNdf)) { + if ((nUpdate > 0) && !extendedSystem.isWellDefined()) { ACTS_INFO("Not enough measurements. Require " - << ndfSystem + 1 << ", but only " << countNdf - << " could be used."); + << extendedSystem.findRequiredNdf() + 1 << ", but only " + << extendedSystem.ndf() << " could be used."); return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // get back the Bound vector components - aMatrix = aMatrixExtended.topLeftCorner().eval(); - bVector = bVectorExtended.topLeftCorner().eval(); - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - aMatrixExtended.colPivHouseholderQr().solve(bVectorExtended); + extendedSystem.aMatrix().colPivHouseholderQr().solve( + extendedSystem.bVector()); deltaParams = deltaParamsExtended.topLeftCorner().eval(); ACTS_VERBOSE("aMatrix:\n" - << aMatrix << "\n" + << extendedSystem.aMatrix() << "\n" << "bVector:\n" - << bVector << "\n" + << extendedSystem.bVector() << "\n" << "deltaParams:\n" << deltaParams << "\n" << "deltaParamsExtended:\n" << deltaParamsExtended << "\n" << "oldChi2sum = " << oldChi2sum << "\n" - << "chi2sum = " << chi2sum); + << "chi2sum = " << extendedSystem.chi2()); if ((gx2fOptions.relChi2changeCutOff != 0) && (nUpdate > 0) && - (std::abs(chi2sum / oldChi2sum - 1) < + (std::abs(extendedSystem.chi2() / oldChi2sum - 1) < gx2fOptions.relChi2changeCutOff)) { ACTS_INFO("Abort with relChi2changeCutOff after " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax << " iterations."); - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } - if (chi2sum > oldChi2sum + 1e-5) { + if (extendedSystem.chi2() > oldChi2sum + 1e-5) { ACTS_DEBUG("chi2 not converging monotonically"); - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } @@ -1378,8 +1436,7 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::DidNotConverge; } - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } @@ -1397,7 +1454,7 @@ class Gx2Fitter { } } - oldChi2sum = chi2sum; + oldChi2sum = extendedSystem.chi2(); } ACTS_DEBUG("Finished to iterate"); ACTS_VERBOSE("Final parameters: " << params.parameters().transpose()); @@ -1445,12 +1502,12 @@ class Gx2Fitter { auto& r = propagatorState.template get>(); r.fittedStates = &trackContainer.trackStateContainer(); - auto propagationResult = m_propagator.template propagate(propagatorState); + auto propagationResult = m_propagator.propagate(propagatorState); // Run the fitter - auto result = m_propagator.template makeResult(std::move(propagatorState), - propagationResult, - propagatorOptions, false); + auto result = + m_propagator.makeResult(std::move(propagatorState), propagationResult, + propagatorOptions, false); if (!result.ok()) { ACTS_ERROR("Propagation failed: " << result.error()); diff --git a/Core/include/Acts/TrackFitting/KalmanFitter.hpp b/Core/include/Acts/TrackFitting/KalmanFitter.hpp index b0fe848c17b..12b70c8a562 100644 --- a/Core/include/Acts/TrackFitting/KalmanFitter.hpp +++ b/Core/include/Acts/TrackFitting/KalmanFitter.hpp @@ -1233,7 +1233,7 @@ class KalmanFitter { track_container_t& trackContainer) const -> Result { auto propagatorState = - m_propagator.template makeState(sParameters, propagatorOptions); + m_propagator.makeState(sParameters, propagatorOptions); auto& kalmanResult = propagatorState.template get>(); diff --git a/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp index 019865f49c1..b7ea835e460 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp @@ -15,6 +15,7 @@ #include "Acts/Utilities/detail/periodic.hpp" #include +#include #include #include @@ -90,10 +91,10 @@ auto gaussianMixtureCov(const components_t components, // Apply corrections for cyclic coordinates auto handleCyclicCov = [&l = pars_l, &m = mean, &diff = diff](auto desc) { - diff[desc.idx] = - difference_periodic(l[desc.idx] / desc.constant, - m[desc.idx] / desc.constant, 2 * M_PI) * - desc.constant; + diff[desc.idx] = difference_periodic(l[desc.idx] / desc.constant, + m[desc.idx] / desc.constant, + 2 * std::numbers::pi) * + desc.constant; }; std::apply([&](auto... dsc) { (handleCyclicCov(dsc), ...); }, angleDesc); @@ -174,10 +175,10 @@ auto gaussianMixtureMeanCov(const components_t components, &weight = weight_l, &mean = mean](auto desc) { const auto delta = (ref[desc.idx] - pars[desc.idx]) / desc.constant; - if (delta > M_PI) { - mean[desc.idx] += (2 * M_PI) * weight * desc.constant; - } else if (delta < -M_PI) { - mean[desc.idx] -= (2 * M_PI) * weight * desc.constant; + if (delta > std::numbers::pi) { + mean[desc.idx] += 2. * std::numbers::pi * weight * desc.constant; + } else if (delta < -std::numbers::pi) { + mean[desc.idx] -= 2. * std::numbers::pi * weight * desc.constant; } }; @@ -187,9 +188,9 @@ auto gaussianMixtureMeanCov(const components_t components, mean /= sumOfWeights; auto wrap = [&](auto desc) { - mean[desc.idx] = - wrap_periodic(mean[desc.idx] / desc.constant, -M_PI, 2 * M_PI) * - desc.constant; + mean[desc.idx] = wrap_periodic(mean[desc.idx] / desc.constant, + -std::numbers::pi, 2 * std::numbers::pi) * + desc.constant; }; std::apply([&](auto... dsc) { (wrap(dsc), ...); }, angleDesc); diff --git a/Core/include/Acts/Utilities/AngleHelpers.hpp b/Core/include/Acts/Utilities/AngleHelpers.hpp index 02feec2323e..8a38aeffc94 100644 --- a/Core/include/Acts/Utilities/AngleHelpers.hpp +++ b/Core/include/Acts/Utilities/AngleHelpers.hpp @@ -9,26 +9,68 @@ #pragma once #include +#include namespace Acts::AngleHelpers { +template +struct EtaThetaConversionTraits {}; + +template <> +struct EtaThetaConversionTraits { + static constexpr float minTheta = 1e-6f; + static constexpr float maxTheta = std::numbers::pi_v - minTheta; + + static constexpr float maxEta = 80.0f; + static constexpr float minEta = -maxEta; +}; + +template <> +struct EtaThetaConversionTraits { + 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 Scalar etaFromTheta(Scalar theta) { - return -std::log(std::tan(0.5 * theta)); + if (theta <= EtaThetaConversionTraits::minTheta) { + return std::numeric_limits::infinity(); + } + if (theta >= EtaThetaConversionTraits::maxTheta) { + return -std::numeric_limits::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 Scalar thetaFromEta(Scalar eta) { + if (eta <= EtaThetaConversionTraits::minEta) { + return std::numbers::pi_v; + } + if (eta >= EtaThetaConversionTraits::maxEta) { + return 0; + } + return 2 * std::atan(std::exp(-eta)); } diff --git a/Core/include/Acts/Utilities/BinAdjustmentVolume.hpp b/Core/include/Acts/Utilities/BinAdjustmentVolume.hpp index 7e115e7b607..bef7bff3e83 100644 --- a/Core/include/Acts/Utilities/BinAdjustmentVolume.hpp +++ b/Core/include/Acts/Utilities/BinAdjustmentVolume.hpp @@ -19,6 +19,7 @@ #include "Acts/Geometry/Volume.hpp" #include "Acts/Utilities/BinUtility.hpp" +#include #include namespace Acts { @@ -94,8 +95,8 @@ BinUtility adjustBinUtility(const BinUtility& bu, // The parameters from the cutout cylinder bounds double minR = cBounds.get(CutoutCylinderVolumeBounds::eMinR); double maxR = cBounds.get(CutoutCylinderVolumeBounds::eMaxR); - double minPhi = -M_PI; - double maxPhi = M_PI; + double minPhi = -std::numbers::pi; + double maxPhi = std::numbers::pi; double minZ = -cBounds.get(CutoutCylinderVolumeBounds::eHalfLengthZ); double maxZ = cBounds.get(CutoutCylinderVolumeBounds::eHalfLengthZ); // Retrieve the binning data diff --git a/Core/include/Acts/Utilities/Frustum.ipp b/Core/include/Acts/Utilities/Frustum.ipp index 5854ef2d94a..bca40c35464 100644 --- a/Core/include/Acts/Utilities/Frustum.ipp +++ b/Core/include/Acts/Utilities/Frustum.ipp @@ -8,6 +8,8 @@ #include "Acts/Utilities/VectorHelpers.hpp" +#include + template Acts::Frustum::Frustum(const VertexType& origin, const VertexType& dir, @@ -17,13 +19,13 @@ Acts::Frustum::Frustum(const VertexType& origin, using rotation_t = Eigen::Rotation2D; static_assert(SIDES == 2, "2D frustum can only have 2 sides"); - assert(opening_angle < M_PI); + assert(opening_angle < std::numbers::pi_v); translation_t translation(origin); value_type angle = VectorHelpers::phi(dir); Eigen::Rotation2D rot(angle); - value_type normal_angle = 0.5 * M_PI - 0.5 * opening_angle; + value_type normal_angle = std::numbers::pi / 2. - opening_angle / 2.; VertexType normal1 = rotation_t(normal_angle) * VertexType::UnitX(); VertexType normal2 = rotation_t(-normal_angle) * VertexType::UnitX(); @@ -37,7 +39,7 @@ Acts::Frustum::Frustum(const VertexType& origin, requires(DIM == 3) : m_origin(origin) { static_assert(SIDES > 2, "3D frustum must have 3 or more sides"); - assert(opening_angle < M_PI); + assert(opening_angle < std::numbers::pi_v); using angle_axis_t = Eigen::AngleAxis; const VertexType ldir = VertexType::UnitZ(); @@ -48,7 +50,7 @@ Acts::Frustum::Frustum(const VertexType& origin, m_normals[0] = ldir; - const value_type phi_sep = 2 * M_PI / sides; + const value_type phi_sep = 2. * std::numbers::pi_v / sides; transform_type rot; rot = angle_axis_t(phi_sep, ldir); diff --git a/Core/include/Acts/Utilities/GraphViz.hpp b/Core/include/Acts/Utilities/GraphViz.hpp index 02cf9487d8d..8716fca7d00 100644 --- a/Core/include/Acts/Utilities/GraphViz.hpp +++ b/Core/include/Acts/Utilities/GraphViz.hpp @@ -92,7 +92,7 @@ std::ostream& operator<<(std::ostream& os, const Style& style); struct Node { std::string id; - std::string label; + std::string label = ""; Shape shape = Shape::Ellipse; std::vector