From bf36f464f98d32fd97dc9833013ed14f1c64166e Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 17 Oct 2024 19:29:54 +0200 Subject: [PATCH 1/8] feat: Allow reflection of track parameters (#3682) Add helper functions to reflect track parameters. This is useful when seed parameters are estimated for the spacepoints in reverse direction which is necessary for strip seeds. --- This pull request introduces functionality to reflect track parameters for both bound and free track parameters in the Acts framework. The key changes include adding new methods to reflect parameters and the corresponding implementation to handle the reflection logic. ### New Methods for Reflecting Track Parameters: * [`Core/include/Acts/EventData/GenericBoundTrackParameters.hpp`](diffhunk://#diff-b3a2aaa0a0f138a1b44e7e6494aae4f8b905ad83716a4253bbea6446623f0eb1R253-R263): Added `reflectInplace` and `reflect` methods to `GenericBoundTrackParameters` to enable in-place reflection and returning reflected parameters. * [`Core/include/Acts/EventData/GenericFreeTrackParameters.hpp`](diffhunk://#diff-35a3d07a1eb05110c620e3cdbc38968456e66270c9f9460472f49e65d8b43acdR177-R187): Added `reflectInplace` and `reflect` methods to `GenericFreeTrackParameters` to enable in-place reflection and returning reflected parameters. ### Helper Functions for Reflection: * [`Core/include/Acts/EventData/TransformationHelpers.hpp`](diffhunk://#diff-73d04b62535d54171256def2e580c39602c9a2ab36db1689ebd74d84b79c127dR21-R32): Introduced `reflectBoundParameters` and `reflectFreeParameters` functions to handle the reflection logic for bound and free track parameters. * [`Core/src/EventData/TransformationHelpers.cpp`](diffhunk://#diff-59ee8d8625a58700aca8eae78a6655d0162483f2644cdd2c77c668493d7d8b87R16-R40): Implemented the `reflectBoundParameters` and `reflectFreeParameters` functions to perform the actual reflection calculations. --- .../EventData/GenericBoundTrackParameters.hpp | 11 ++++++++ .../GenericCurvilinearTrackParameters.hpp | 8 ++++++ .../EventData/GenericFreeTrackParameters.hpp | 12 ++++++++ .../Acts/EventData/TransformationHelpers.hpp | 28 +++++++++++++++++++ .../EventData/BoundTrackParametersTests.cpp | 8 ++++++ .../CurvilinearTrackParametersTests.cpp | 9 ++++++ .../EventData/FreeTrackParametersTests.cpp | 8 ++++++ 7 files changed, 84 insertions(+) diff --git a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp index 27b411a99cc..f624cd88aee 100644 --- a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp @@ -250,6 +250,17 @@ class GenericBoundTrackParameters { return m_surface->referenceFrame(geoCtx, position(geoCtx), momentum()); } + /// Reflect the parameters in place. + void reflectInPlace() { m_params = reflectBoundParameters(m_params); } + + /// Reflect the parameters. + /// @return Reflected parameters. + GenericBoundTrackParameters reflect() const { + GenericBoundTrackParameters reflected = *this; + reflected.reflectInPlace(); + return reflected; + } + private: BoundVector m_params; std::optional m_cov; diff --git a/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp b/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp index a42ac2f116f..6e4ed9bae7a 100644 --- a/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp @@ -111,6 +111,14 @@ class GenericCurvilinearTrackParameters Vector3 position() const { return GenericBoundTrackParameters::position({}); } + + /// Reflect the parameters. + /// @return Reflected parameters. + GenericCurvilinearTrackParameters reflect() const { + GenericCurvilinearTrackParameters reflected = *this; + reflected.reflectInPlace(); + return reflected; + } }; } // namespace Acts diff --git a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp index 1e7846318ef..857d6c328b9 100644 --- a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp @@ -12,6 +12,7 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/TrackParametersConcept.hpp" +#include "Acts/EventData/TransformationHelpers.hpp" #include "Acts/EventData/detail/PrintParameters.hpp" #include "Acts/Utilities/MathHelpers.hpp" #include "Acts/Utilities/UnitVectors.hpp" @@ -175,6 +176,17 @@ class GenericFreeTrackParameters { return m_particleHypothesis; } + /// Reflect the parameters in place. + void reflectInPlace() { m_params = reflectFreeParameters(m_params); } + + /// Reflect the parameters. + /// @return Reflected parameters. + GenericFreeTrackParameters reflect() const { + GenericFreeTrackParameters reflected = *this; + reflected.reflectInPlace(); + return reflected; + } + private: FreeVector m_params; std::optional m_cov; diff --git a/Core/include/Acts/EventData/TransformationHelpers.hpp b/Core/include/Acts/EventData/TransformationHelpers.hpp index 7fa297f4abc..39240bf469e 100644 --- a/Core/include/Acts/EventData/TransformationHelpers.hpp +++ b/Core/include/Acts/EventData/TransformationHelpers.hpp @@ -13,11 +13,39 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Utilities/Result.hpp" +#include "Acts/Utilities/detail/periodic.hpp" namespace Acts { class Surface; +/// Reflect bound track parameters. +/// +/// @param boundParams Bound track parameters vector +/// @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]); + reflected[eBoundPhi] = phi; + reflected[eBoundTheta] = theta; + reflected[eBoundQOverP] = -boundParams[eBoundQOverP]; + return reflected; +} + +/// Reflect free track parameters. +/// +/// @param freeParams Free track parameters vector +/// @return Reflected free track parameters vector +inline FreeVector reflectFreeParameters(const FreeVector& freeParams) { + FreeVector reflected = freeParams; + reflected[eFreeDir0] = -freeParams[eFreeDir0]; + reflected[eFreeDir1] = -freeParams[eFreeDir1]; + reflected[eFreeDir2] = -freeParams[eFreeDir2]; + reflected[eFreeQOverP] = -freeParams[eFreeQOverP]; + return reflected; +} + /// Transform bound track parameters into equivalent free track parameters. /// /// @param surface Surface onto which the input parameters are bound diff --git a/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp index a05d273f7a1..9daf24c8a79 100644 --- a/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp +++ b/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp @@ -76,6 +76,14 @@ void checkParameters(const BoundTrackParameters& params, double l0, double l1, eps); CHECK_CLOSE_OR_SMALL(params.momentum(), p * unitDir, eps, eps); BOOST_CHECK_EQUAL(params.charge(), q); + + // reflection + BoundTrackParameters reflectedParams = params; + reflectedParams.reflectInPlace(); + CHECK_CLOSE_OR_SMALL(params.reflect().parameters(), + reflectedParams.parameters(), eps, eps); + CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(), + params.parameters(), eps, eps); } void runTest(const std::shared_ptr& surface, double l0, diff --git a/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp index 122d4b075d0..143515148d4 100644 --- a/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp +++ b/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp @@ -72,6 +72,15 @@ void checkParameters(const CurvilinearTrackParameters& params, double phi, // curvilinear reference surface CHECK_CLOSE_OR_SMALL(referenceSurface->center(geoCtx), pos, eps, eps); CHECK_CLOSE_OR_SMALL(referenceSurface->normal(geoCtx), unitDir, eps, eps); + + // reflection + CurvilinearTrackParameters reflectedParams = params; + reflectedParams.reflectInPlace(); + CHECK_CLOSE_OR_SMALL(params.reflect().parameters(), + reflectedParams.parameters(), eps, eps); + CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(), + params.parameters(), eps, eps); + // TODO verify reference frame } diff --git a/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp index 71b0dc0a20e..fb7a11f3f68 100644 --- a/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp +++ b/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp @@ -67,6 +67,14 @@ void checkParameters(const FreeTrackParameters& params, const Vector4& pos4, eps); CHECK_CLOSE_OR_SMALL(params.time(), params.template get(), eps, eps); + + // reflection + FreeTrackParameters reflectedParams = params; + reflectedParams.reflectInPlace(); + CHECK_CLOSE_OR_SMALL(params.reflect().parameters(), + reflectedParams.parameters(), eps, eps); + CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(), + params.parameters(), eps, eps); } } // namespace From 057d33d60918ef7e3c37b291827d15dfc6fb04a3 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 17 Oct 2024 21:48:39 +0200 Subject: [PATCH 2/8] feat: Add `estimateTrackParamCovariance` to Core (#3683) Add a `estimateTrackParamCovariance` function to Core to estimate the initial covariance of track parameters after seeding. blocked by - https://github.com/acts-project/acts/pull/3665 --- This pull request introduces a new function for estimating track parameter covariance and refactors existing code to utilize this new function. The changes improve code modularity and maintainability by centralizing the covariance estimation logic. ### New Functionality: * [`Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp`](diffhunk://#diff-1d63df6364bcd57d6eae717e4c8dea355f07a36f5fca41b0e9fefdbf3224c0a0R292-R335): Added `EstimateTrackParamCovarianceConfig` struct and `estimateTrackParamCovariance` function to estimate track parameter covariance. ### Refactoring: * [`Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp`](diffhunk://#diff-3fa61f411627effaf29755fa1d539c90bbd2c5155f04c4fa7a0d7f36d99340f3L38-L79): Removed the `makeInitialCovariance` function and replaced its usage with the new `estimateTrackParamCovariance` function. This change simplifies the code by removing redundant logic. [[1]](diffhunk://#diff-3fa61f411627effaf29755fa1d539c90bbd2c5155f04c4fa7a0d7f36d99340f3L38-L79) [[2]](diffhunk://#diff-3fa61f411627effaf29755fa1d539c90bbd2c5155f04c4fa7a0d7f36d99340f3L174-R144) * [`Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp`](diffhunk://#diff-d4047681edb055a3b4ca13bb837f6fef194fd0c8b3205f337db73e97f4b6de82L131-R141): Updated the covariance calculation to use the new `estimateTrackParamCovariance` function, ensuring consistency across the codebase. ### Miscellaneous: * [`Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp`](diffhunk://#diff-d4047681edb055a3b4ca13bb837f6fef194fd0c8b3205f337db73e97f4b6de82R14): Included the `EstimateTrackParamsFromSeed.hpp` header to access the new covariance estimation function. --- Core/CMakeLists.txt | 1 + .../Seeding/EstimateTrackParamsFromSeed.hpp | 37 +++++++++++++ Core/src/Seeding/CMakeLists.txt | 1 + .../Seeding/EstimateTrackParamsFromSeed.cpp | 50 +++++++++++++++++ .../src/TrackParamsEstimationAlgorithm.cpp | 53 ++++--------------- .../TruthTracking/ParticleSmearing.cpp | 32 ++++------- Examples/Python/tests/root_file_hashes.txt | 8 +-- 7 files changed, 111 insertions(+), 71 deletions(-) create mode 100644 Core/src/Seeding/CMakeLists.txt create mode 100644 Core/src/Seeding/EstimateTrackParamsFromSeed.cpp diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt index eeb6bc8b92a..554656f514c 100644 --- a/Core/CMakeLists.txt +++ b/Core/CMakeLists.txt @@ -110,6 +110,7 @@ add_subdirectory(src/MagneticField) add_subdirectory(src/Material) add_subdirectory(src/Navigation) add_subdirectory(src/Propagator) +add_subdirectory(src/Seeding) add_subdirectory(src/Surfaces) add_subdirectory(src/TrackFinding) add_subdirectory(src/TrackFitting) diff --git a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp index 7a5c597f626..f3850e102ce 100644 --- a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp +++ b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp @@ -287,4 +287,41 @@ std::optional estimateTrackParamsFromSeed( return params; } +/// Configuration for the estimation of the covariance matrix of the track +/// parameters with `estimateTrackParamCovariance`. +struct EstimateTrackParamCovarianceConfig { + /// The initial sigmas for the track parameters + BoundVector initialSigmas = {1. * UnitConstants::mm, + 1. * UnitConstants::mm, + 1. * UnitConstants::degree, + 1. * UnitConstants::degree, + 1. * UnitConstants::e / UnitConstants::GeV, + 1. * UnitConstants::ns}; + + /// The initial relative uncertainty of the q/pt + double initialSigmaPtRel = 0.1; + + /// The inflation factors for the variances of the track parameters + BoundVector initialVarInflation = {1., 1., 1., 1., 1., 1.}; + /// The inflation factor for time uncertainty if the time parameter was not + /// estimated + double noTimeVarInflation = 100.; +}; + +/// Estimate the covariance matrix of the given track parameters based on the +/// provided configuration. The assumption is that we can model the uncertainty +/// of the track parameters as a diagonal matrix with the provided initial +/// sigmas. The inflation factors are used to inflate the initial variances +/// based on the provided configuration. The uncertainty of q/p is estimated +/// based on the relative uncertainty of the q/pt and the theta uncertainty. +/// +/// @param config is the configuration for the estimation +/// @param params is the track parameters +/// @param hasTime is true if the track parameters have time +/// +/// @return the covariance matrix of the track parameters +BoundMatrix estimateTrackParamCovariance( + const EstimateTrackParamCovarianceConfig& config, const BoundVector& params, + bool hasTime); + } // namespace Acts diff --git a/Core/src/Seeding/CMakeLists.txt b/Core/src/Seeding/CMakeLists.txt new file mode 100644 index 00000000000..770037b1dfd --- /dev/null +++ b/Core/src/Seeding/CMakeLists.txt @@ -0,0 +1 @@ +target_sources(ActsCore PRIVATE EstimateTrackParamsFromSeed.cpp) diff --git a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp new file mode 100644 index 00000000000..83d950fc3a1 --- /dev/null +++ b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp @@ -0,0 +1,50 @@ +// 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 "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" + +#include "Acts/Definitions/TrackParametrization.hpp" + +Acts::BoundMatrix Acts::estimateTrackParamCovariance( + const EstimateTrackParamCovarianceConfig& config, const BoundVector& params, + bool hasTime) { + assert((params[eBoundTheta] > 0 && params[eBoundTheta] < M_PI) && + "Theta must be in the range (0, pi)"); + + BoundSquareMatrix result = BoundSquareMatrix::Zero(); + + for (std::size_t i = eBoundLoc0; i < eBoundSize; ++i) { + double sigma = config.initialSigmas[i]; + double variance = sigma * sigma; + + if (i == eBoundQOverP) { + // note that we rely on the fact that sigma theta is already computed + double varianceTheta = result(eBoundTheta, eBoundTheta); + + // transverse momentum contribution + variance += std::pow(config.initialSigmaPtRel * params[eBoundQOverP], 2); + + // theta contribution + variance += + varianceTheta * + std::pow(params[eBoundQOverP] / std::tan(params[eBoundTheta]), 2); + } + + if (i == eBoundTime && !hasTime) { + // Inflate the time uncertainty if no time measurement is available + variance *= config.noTimeVarInflation; + } + + // Inflate the initial covariance + variance *= config.initialVarInflation[i]; + + result(i, i) = variance; + } + + return result; +} diff --git a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp index a54b10d2e7d..c8a30c7147d 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp @@ -32,48 +32,6 @@ namespace ActsExamples { -namespace { - -Acts::BoundSquareMatrix makeInitialCovariance( - const TrackParamsEstimationAlgorithm::Config& config, - const Acts::BoundVector& params, const SimSpacePoint& sp) { - Acts::BoundSquareMatrix result = Acts::BoundSquareMatrix::Zero(); - - for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) { - double sigma = config.initialSigmas[i]; - double variance = sigma * sigma; - - if (i == Acts::eBoundQOverP) { - // note that we rely on the fact that sigma theta is already computed - double varianceTheta = result(Acts::eBoundTheta, Acts::eBoundTheta); - - // transverse momentum contribution - variance += - std::pow(config.initialSigmaPtRel * params[Acts::eBoundQOverP], 2); - - // theta contribution - variance += - varianceTheta * std::pow(params[Acts::eBoundQOverP] / - std::tan(params[Acts::eBoundTheta]), - 2); - } - - // Inflate the time uncertainty if no time measurement is available - if (i == Acts::eBoundTime && !sp.t().has_value()) { - variance *= config.noTimeVarInflation; - } - - // Inflate the initial covariance - variance *= config.initialVarInflation[i]; - - result(i, i) = variance; - } - - return result; -} - -} // namespace - TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm( TrackParamsEstimationAlgorithm::Config cfg, Acts::Logging::Level lvl) : IAlgorithm("TrackParamsEstimationAlgorithm", lvl), m_cfg(std::move(cfg)) { @@ -166,8 +124,15 @@ ProcessCode TrackParamsEstimationAlgorithm::execute( const auto& params = optParams.value(); - Acts::BoundSquareMatrix cov = - makeInitialCovariance(m_cfg, params, *bottomSP); + Acts::EstimateTrackParamCovarianceConfig config{ + .initialSigmas = + Eigen::Map{m_cfg.initialSigmas.data()}, + .initialSigmaPtRel = m_cfg.initialSigmaPtRel, + .initialVarInflation = Eigen::Map{ + m_cfg.initialVarInflation.data()}}; + + Acts::BoundSquareMatrix cov = Acts::estimateTrackParamCovariance( + config, params, bottomSP->t().has_value()); trackParameters.emplace_back(surface->getSharedPtr(), params, cov, m_cfg.particleHypothesis); diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp index 937a5bc1863..cfc62ab9da9 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" +#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Helpers.hpp" @@ -128,31 +129,16 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute( Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero(); if (m_cfg.initialSigmas) { // use the initial sigmas if set - for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) { - double sigma = (*m_cfg.initialSigmas)[i]; - double variance = sigma * sigma; - - if (i == Acts::eBoundQOverP) { - // note that we rely on the fact that sigma theta is already - // computed - double varianceTheta = cov(Acts::eBoundTheta, Acts::eBoundTheta); - // transverse momentum contribution - variance += std::pow( - m_cfg.initialSigmaPtRel * params[Acts::eBoundQOverP], 2); + Acts::EstimateTrackParamCovarianceConfig config{ + .initialSigmas = + Eigen::Map{ + m_cfg.initialSigmas->data()}, + .initialSigmaPtRel = m_cfg.initialSigmaPtRel, + .initialVarInflation = Eigen::Map{ + m_cfg.initialVarInflation.data()}}; - // theta contribution - variance += varianceTheta * - std::pow(params[Acts::eBoundQOverP] / - std::tan(params[Acts::eBoundTheta]), - 2); - } - - // Inflate the initial covariance - variance *= m_cfg.initialVarInflation[i]; - - cov(i, i) = variance; - } + cov = Acts::estimateTrackParamCovariance(config, params, false); } else { // otherwise use the smearing sigmas diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 6ddb1a405cd..39df454b3e8 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -39,16 +39,16 @@ test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0 test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: a8c5c6f6c1e6303b887d47b509b7f71a2ffa5f38638fe46ce5bce76fd20d64ca test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 417f7326e1e1bb4519f1378145ac733bdda6653eb9871fd69e455e0269d996a6 test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e -test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: de11c0868a70ade0dcc80465d4e6dcf1dd7fcf8149603b47ee7d87d862a6534a -test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: f18e9ecce6d9585fd150c5aafc9ac225a5bab342aaab50a28283ba879691af1f +test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: edf0b06ce9ee0e4fcb153e41859af7b5153271de18f49a6842a23ad2d66b7e09 +test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 06d6ae1d05cb611b19df3c59531997c9b0108f5ef6027d76c4827bd2d9edb921 test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 463d6aaed4d869652b5b184940e789cde0fb441bdd135813b85462a515e6480a test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: a8ad83a07b48d4cfcf70d0e6fdc3c8997eb03c1f8c2a7be27ea888b099000d79 test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 247dd581cc177625c0286718261c004e2149536d70c8281dfaf697879a84d76d test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 1b08a80e73aedf5cf38a3a407794b82297bec37f556ad4efcda3489a1b17d4d2 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575 -test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 7adfc2bf5ee35a126b713187dd8b11f4497cf864a4a83e57a40885688974413e -test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 7a9de8a8bd1c09f7b4d1c547f824af6c8123afb044dd429180b0d13e47d6f975 +test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: a9621b535ea2912d172142394f51f68e4e7dc255b32d479d6305fa599152b420 +test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: af1a6bb16a070db7ed8043e2188d56f0034843099fc3c332731c4cf86ba39c57 test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3 From c862d2da097451a197b931c810c3ad667f1586c1 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 18 Oct 2024 00:27:00 +0200 Subject: [PATCH 3/8] fix: Handle Fatras immediate abort edge case (#3744) Currently Fatras might abort before reaching the pointer where we initialize the particle. This leads to desynchronization between initial and final particles which will then lead to warning/error logs downstream ``` ParticleSele WARNING No final particle found for 4|11|169|0|0 ``` The initial call of the actor is now caught be checking the propagation stage which then initializes the result object gracefully. Discovered in https://github.com/acts-project/acts/pull/3742 --- .../include/ActsFatras/Kernel/Simulation.hpp | 10 +++++- .../Kernel/detail/SimulationActor.hpp | 31 +++++++--------- .../Fatras/Kernel/SimulationActorTests.cpp | 36 +++++++++++++++++++ 3 files changed, 57 insertions(+), 20 deletions(-) diff --git a/Fatras/include/ActsFatras/Kernel/Simulation.hpp b/Fatras/include/ActsFatras/Kernel/Simulation.hpp index cbc1115225a..2d8280737f0 100644 --- a/Fatras/include/ActsFatras/Kernel/Simulation.hpp +++ b/Fatras/include/ActsFatras/Kernel/Simulation.hpp @@ -245,6 +245,9 @@ struct Simulation { continue; } + assert(result->particle.particleId() == initialParticle.particleId() && + "Particle id must not change during simulation"); + copyOutputs(result.value(), simulatedParticlesInitial, simulatedParticlesFinal, hits); // since physics processes are independent, there can be particle id @@ -256,6 +259,10 @@ struct Simulation { } } + assert( + (simulatedParticlesInitial.size() == simulatedParticlesFinal.size()) && + "Inconsistent final sizes of the simulated particle containers"); + // the overall function call succeeded, i.e. no fatal errors occurred. // yet, there might have been some particle for which the propagation // failed. thus, the successful result contains a list of failed particles. @@ -284,12 +291,13 @@ struct Simulation { // initial particle state was already pushed to the container before // store final particle state at the end of the simulation particlesFinal.push_back(result.particle); + std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits)); + // move generated secondaries that should be simulated to the output std::copy_if( result.generatedParticles.begin(), result.generatedParticles.end(), std::back_inserter(particlesInitial), [this](const Particle &particle) { return selectParticle(particle); }); - std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits)); } /// Renumber particle ids in the tail of the container. diff --git a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp index cc6d7200868..acbaaaa713d 100644 --- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp +++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp @@ -69,7 +69,17 @@ struct SimulationActor { void act(propagator_state_t &state, stepper_t &stepper, navigator_t &navigator, result_type &result, const Acts::Logger &logger) const { - assert(generator && "The generator pointer must be valid"); + assert(generator != nullptr && "The generator pointer must be valid"); + + if (state.stage == Acts::PropagatorStage::prePropagation) { + // first step is special: there is no previous state and we need to arm + // the decay simulation for all future steps. + result.particle = + makeParticle(initialParticle, state, stepper, navigator); + result.properTimeLimit = + decay.generateProperTimeLimit(*generator, initialParticle); + return; + } // actors are called once more after the propagation terminated if (!result.isAlive) { @@ -82,28 +92,11 @@ struct SimulationActor { return; } - // check if we are still on the start surface and skip if so - if ((navigator.startSurface(state.navigation) != nullptr) && - (navigator.startSurface(state.navigation) == - navigator.currentSurface(state.navigation))) { - return; - } - // update the particle state first. this also computes the proper time which // needs the particle state from the previous step for reference. that means // this must happen for every step (not just on surface) and before // everything, e.g. any interactions that could modify the state. - if (std::isnan(result.properTimeLimit)) { - // first step is special: there is no previous state and we need to arm - // the decay simulation for all future steps. - result.particle = - makeParticle(initialParticle, state, stepper, navigator); - result.properTimeLimit = - decay.generateProperTimeLimit(*generator, initialParticle); - } else { - result.particle = - makeParticle(result.particle, state, stepper, navigator); - } + result.particle = makeParticle(result.particle, state, stepper, navigator); // decay check. needs to happen at every step, not just on surfaces. if (std::isfinite(result.properTimeLimit) && diff --git a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp index e52d586c30b..3306c5702d6 100644 --- a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp +++ b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp @@ -244,7 +244,13 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: surface selection -> one hit, no material -> no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -271,6 +277,7 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) { BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum()); // call.actor again: one more hit, still no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -317,7 +324,13 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: surface selection -> one hit, material -> one secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -346,6 +359,7 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) { tol); // call.actor again: one more hit, one more secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -392,7 +406,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: no surface sel. -> no hit, no material -> no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -419,6 +439,7 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum()); // call.actor again: no hit, still no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -455,7 +476,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) { Fixture f(125_MeV, makeMaterialSurface()); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: no surface sel. -> no hit, material -> one secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -483,6 +510,7 @@ BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) { tol); // call.actor again: still no hit, one more secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -523,7 +551,13 @@ BOOST_AUTO_TEST_CASE(Decay) { // inverse Lorentz factor for proper time dilation: 1/gamma = m/E const auto gammaInv = f.m / f.e; + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // first step w/ defaults leaves particle alive + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -536,6 +570,7 @@ BOOST_AUTO_TEST_CASE(Decay) { BOOST_CHECK_EQUAL(f.result.particle.properTime(), 0_ns); // second step w/ defaults increases proper time + f.state.stage = Acts::PropagatorStage::postStep; f.state.stepping.time += 1_ns; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); @@ -549,6 +584,7 @@ BOOST_AUTO_TEST_CASE(Decay) { CHECK_CLOSE_REL(f.result.particle.properTime(), gammaInv * 1_ns, tol); // third step w/ proper time limit decays the particle + f.state.stage = Acts::PropagatorStage::postStep; f.state.stepping.time += 1_ns; f.result.properTimeLimit = f.result.particle.properTime() + gammaInv * 0.5_ns; f.actor.act(f.state, f.stepper, f.navigator, f.result, From 4643170b29d7724339b9c0cb8d8c3e6adebef1b0 Mon Sep 17 00:00:00 2001 From: Tim Adye Date: Fri, 18 Oct 2024 01:47:42 +0100 Subject: [PATCH 4/8] fix: fix nMeasurements() comment (#3740) comments for mutable and const versions got mixed up. --- Core/include/Acts/EventData/TrackProxy.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Core/include/Acts/EventData/TrackProxy.hpp b/Core/include/Acts/EventData/TrackProxy.hpp index ce674e9cfa4..5939022e2e8 100644 --- a/Core/include/Acts/EventData/TrackProxy.hpp +++ b/Core/include/Acts/EventData/TrackProxy.hpp @@ -324,7 +324,8 @@ class TrackProxy { return std::distance(tsRange.begin(), tsRange.end()); } - /// Return the number of measurements for the track. Const version + /// Return a mutable reference to the number of measurements for the track. + /// Mutable version /// @note Only available if the track proxy is not read-only /// @return The number of measurements unsigned int& nMeasurements() @@ -333,8 +334,7 @@ class TrackProxy { return component(); } - /// Return a mutable reference to the number of measurements for the track. - /// Mutable version + /// Return the number of measurements for the track. Const version /// @return The number of measurements unsigned int nMeasurements() const { return component(); From 66d18541523bb00854ed754d190202b553909f50 Mon Sep 17 00:00:00 2001 From: Tim Adye Date: Fri, 18 Oct 2024 03:23:58 +0100 Subject: [PATCH 5/8] fix: TrackSelector remove broken m_noEtaCuts (#3640) `m_noEtaCuts` wasn't set correctly, and can be removed. Also check that eta-binned cuts don't themselves cut on eta. This was only tested previously in the case of a single bin. Also added `binIndexNoCheck(eta)` which is generally useful and saves having to check `eta` instead of the index. I think it would be better to have `getCuts(eta)` not `throw` but just return the first/last bin, but that would change the API - unless you added a no-`throw` variant. `binIndexNoCheck` at least helps in the meantime. --- .../Acts/TrackFinding/TrackSelector.hpp | 63 ++++++++++--------- 1 file changed, 34 insertions(+), 29 deletions(-) diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp index d0e8526dc35..e782b7feb9d 100644 --- a/Core/include/Acts/TrackFinding/TrackSelector.hpp +++ b/Core/include/Acts/TrackFinding/TrackSelector.hpp @@ -142,7 +142,7 @@ class TrackSelector { std::vector cutSets = {}; /// Eta bin edges for varying cuts by eta - std::vector absEtaEdges = {}; + std::vector absEtaEdges = {0, inf}; /// Get the number of eta bins /// @return Number of eta bins @@ -150,7 +150,7 @@ class TrackSelector { /// Construct an empty (accepts everything) configuration. /// Results in a single cut set and one abs eta bin from 0 to infinity. - EtaBinnedConfig() : cutSets{{}}, absEtaEdges{{0, inf}} {}; + EtaBinnedConfig() : cutSets{{}} {}; /// Constructor to create a config object that is not upper-bounded. /// This is useful to use the "fluent" API to populate the configuration. @@ -163,13 +163,12 @@ class TrackSelector { /// @param absEtaEdgesIn is the vector of eta bin edges EtaBinnedConfig(std::vector absEtaEdgesIn) : absEtaEdges{std::move(absEtaEdgesIn)} { - cutSets.resize(absEtaEdges.size() - 1); + cutSets.resize(nEtaBins()); } /// Auto-converting constructor from a single cut configuration. /// Results in a single absolute eta bin from 0 to infinity. - EtaBinnedConfig(Config cutSet) - : cutSets{std::move(cutSet)}, absEtaEdges{{0, inf}} {} + EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {} /// Add a new eta bin with the given upper bound. /// @param etaMax Upper bound of the new eta bin @@ -195,11 +194,17 @@ class TrackSelector { /// @return True if the configuration has a bin for the given eta bool hasCuts(double eta) const; - /// Get the index of the eta bin for a given eta + /// Get the index of the eta bin for a given eta. + /// throws an exception if Eta is outside the abs eta bin edges. /// @param eta Eta value /// @return Index of the eta bin std::size_t binIndex(double eta) const; + /// Get the index of the eta bin for a given eta + /// @param eta Eta value + /// @return Index of the eta bin, or >= nEtaBins() if Eta is outside the abs eta bin edges. + std::size_t binIndexNoCheck(double eta) const; + /// Get the cuts for a given eta /// @param eta Eta value /// @return Cuts for the given eta @@ -237,8 +242,7 @@ class TrackSelector { private: EtaBinnedConfig m_cfg; - bool m_isUnbinned; - bool m_noEtaCuts; + bool m_isUnbinned = false; }; inline TrackSelector::Config& TrackSelector::Config::loc0(double min, @@ -350,14 +354,22 @@ inline bool TrackSelector::EtaBinnedConfig::hasCuts(double eta) const { } inline std::size_t TrackSelector::EtaBinnedConfig::binIndex(double eta) const { - if (!hasCuts(eta)) { + std::size_t index = binIndexNoCheck(eta); + if (!(index < nEtaBins())) { throw std::invalid_argument{"Eta is outside the abs eta bin edges"}; } + return index; +} +inline std::size_t TrackSelector::EtaBinnedConfig::binIndexNoCheck( + double eta) const { auto binIt = std::upper_bound(absEtaEdges.begin(), absEtaEdges.end(), std::abs(eta)); - std::size_t index = std::distance(absEtaEdges.begin(), binIt) - 1; - return index; + std::size_t index = std::distance(absEtaEdges.begin(), binIt); + if (index == 0) { + index = absEtaEdges.size() + 1; // positive value to check for underflow + } + return index - 1; } inline const TrackSelector::Config& TrackSelector::EtaBinnedConfig::getCuts( @@ -428,8 +440,8 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const { return track.hasReferenceSurface() && within(track.transverseMomentum(), cuts.ptMin, cuts.ptMax) && - (m_noEtaCuts || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) && - within(_eta, cuts.etaMin, cuts.etaMax))) && + (!m_isUnbinned || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) && + within(_eta, cuts.etaMin, cuts.etaMax))) && within(track.phi(), cuts.phiMin, cuts.phiMax) && within(track.loc0(), cuts.loc0Min, cuts.loc0Max) && within(track.loc1(), cuts.loc1Min, cuts.loc1Max) && @@ -452,26 +464,19 @@ inline TrackSelector::TrackSelector( "TrackSelector cut / eta bin configuration is inconsistent"}; } - m_isUnbinned = false; if (m_cfg.nEtaBins() == 1) { static const std::vector infVec = {0, inf}; - bool limitEta = m_cfg.absEtaEdges != infVec; - m_isUnbinned = !limitEta; // single bin, no eta edges given - - const Config& cuts = m_cfg.cutSets[0]; - - if (limitEta && (cuts.etaMin != -inf || cuts.etaMax != inf || - cuts.absEtaMin != 0.0 || cuts.absEtaMax != inf)) { - throw std::invalid_argument{ - "Explicit eta cuts are only valid for single eta bin"}; - } + m_isUnbinned = + m_cfg.absEtaEdges == infVec; // single bin, no eta edges given } - m_noEtaCuts = m_isUnbinned; - for (const auto& cuts : m_cfg.cutSets) { - if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 || - cuts.absEtaMax != inf) { - m_noEtaCuts = false; + if (!m_isUnbinned) { + for (const auto& cuts : m_cfg.cutSets) { + if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 || + cuts.absEtaMax != inf) { + throw std::invalid_argument{ + "Explicit eta cuts are only valid for single eta bin"}; + } } } } From fbf332b17bf6f339b5a63fa206cc9b687f439fb7 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Fri, 18 Oct 2024 05:53:07 +0200 Subject: [PATCH 6/8] refactor: explicit `hasStackTraces()` in FPE (#3727) Also fixes a typo --- .../Framework/src/Framework/Sequencer.cpp | 2 +- .../Acts/Plugins/FpeMonitoring/FpeMonitor.hpp | 6 ++-- Plugins/FpeMonitoring/src/FpeMonitor.cpp | 34 +++++++++---------- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/Examples/Framework/src/Framework/Sequencer.cpp b/Examples/Framework/src/Framework/Sequencer.cpp index 78444c58ac4..f6bb69d44c0 100644 --- a/Examples/Framework/src/Framework/Sequencer.cpp +++ b/Examples/Framework/src/Framework/Sequencer.cpp @@ -611,7 +611,7 @@ void Sequencer::fpeReport() const { auto merged = std::accumulate( fpe.begin(), fpe.end(), Acts::FpeMonitor::Result{}, [](const auto& lhs, const auto& rhs) { return lhs.merged(rhs); }); - if (!merged) { + if (!merged.hasStackTraces()) { // no FPEs to report continue; } diff --git a/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp b/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp index d1b9fabec64..022a07a3426 100644 --- a/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp +++ b/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp @@ -40,7 +40,7 @@ std::ostream &operator<<(std::ostream &os, FpeType type); class FpeMonitor { public: struct Buffer { - Buffer(std::size_t bufferSize) + explicit Buffer(std::size_t bufferSize) : m_data{std::make_unique(bufferSize)}, m_size{bufferSize} {} @@ -105,12 +105,12 @@ class FpeMonitor { Result() = default; - operator bool() const { return !m_stracktraces.empty(); } + bool hasStackTraces() const { return !m_stackTraces.empty(); } void add(Acts::FpeType type, void *stackPtr, std::size_t bufferSize); private: - std::vector m_stracktraces; + std::vector m_stackTraces; std::array m_counts{}; friend FpeMonitor; diff --git a/Plugins/FpeMonitoring/src/FpeMonitor.cpp b/Plugins/FpeMonitoring/src/FpeMonitor.cpp index db41dc1f8b2..09dce1f4b0a 100644 --- a/Plugins/FpeMonitoring/src/FpeMonitor.cpp +++ b/Plugins/FpeMonitoring/src/FpeMonitor.cpp @@ -61,10 +61,10 @@ FpeMonitor::Result FpeMonitor::Result::merged(const Result &with) const { result.m_counts[i] = m_counts[i] + with.m_counts[i]; } - std::copy(with.m_stracktraces.begin(), with.m_stracktraces.end(), - std::back_inserter(result.m_stracktraces)); - std::copy(m_stracktraces.begin(), m_stracktraces.end(), - std::back_inserter(result.m_stracktraces)); + std::copy(with.m_stackTraces.begin(), with.m_stackTraces.end(), + std::back_inserter(result.m_stackTraces)); + std::copy(m_stackTraces.begin(), m_stackTraces.end(), + std::back_inserter(result.m_stackTraces)); result.deduplicate(); @@ -76,8 +76,8 @@ void FpeMonitor::Result::merge(const Result &with) { m_counts[i] = m_counts[i] + with.m_counts[i]; } - std::copy(with.m_stracktraces.begin(), with.m_stracktraces.end(), - std::back_inserter(m_stracktraces)); + std::copy(with.m_stackTraces.begin(), with.m_stackTraces.end(), + std::back_inserter(m_stackTraces)); deduplicate(); } @@ -87,20 +87,20 @@ void FpeMonitor::Result::add(FpeType type, void *stackPtr, auto st = std::make_unique( boost::stacktrace::stacktrace::from_dump(stackPtr, bufferSize)); - auto it = std::ranges::find_if(m_stracktraces, [&](const FpeInfo &el) { + auto it = std::ranges::find_if(m_stackTraces, [&](const FpeInfo &el) { return areFpesEquivalent({el.type, *el.st}, {type, *st}); }); - if (it != m_stracktraces.end()) { + if (it != m_stackTraces.end()) { it->count += 1; } else { - m_stracktraces.push_back({1, type, std::move(st)}); + m_stackTraces.push_back({1, type, std::move(st)}); } } bool FpeMonitor::Result::contains( FpeType type, const boost::stacktrace::stacktrace &st) const { - return std::ranges::any_of(m_stracktraces, [&](const FpeInfo &el) { + return std::ranges::any_of(m_stackTraces, [&](const FpeInfo &el) { return areFpesEquivalent({el.type, *el.st}, {type, st}); }); } @@ -128,12 +128,12 @@ unsigned int FpeMonitor::Result::count(FpeType type) const { } unsigned int FpeMonitor::Result::numStackTraces() const { - return m_stracktraces.size(); + return m_stackTraces.size(); } const std::vector & FpeMonitor::Result::stackTraces() const { - return m_stracktraces; + return m_stackTraces; } bool FpeMonitor::Result::encountered(FpeType type) const { @@ -161,18 +161,18 @@ void FpeMonitor::Result::summary(std::ostream &os, std::size_t depth) const { void FpeMonitor::Result::deduplicate() { std::vector copy{}; - copy = std::move(m_stracktraces); - m_stracktraces.clear(); + copy = std::move(m_stackTraces); + m_stackTraces.clear(); for (auto &info : copy) { - auto it = std::ranges::find_if(m_stracktraces, [&info](const FpeInfo &el) { + auto it = std::ranges::find_if(m_stackTraces, [&info](const FpeInfo &el) { return areFpesEquivalent({el.type, *el.st}, {info.type, *info.st}); }); - if (it != m_stracktraces.end()) { + if (it != m_stackTraces.end()) { it->count += info.count; continue; } - m_stracktraces.push_back({info.count, info.type, std::move(info.st)}); + m_stackTraces.push_back({info.count, info.type, std::move(info.st)}); } } From c06a60ad4407f30d1121868e8d35cd945b79c11c Mon Sep 17 00:00:00 2001 From: ssdetlab <113530373+ssdetlab@users.noreply.github.com> Date: Fri, 18 Oct 2024 09:59:30 +0300 Subject: [PATCH 7/8] feat: Track parameters json converter (#3724) Adding json converter for the different track parameters types. The converter follows the pattern outlined in [nlohmann::json documentation](https://json.nlohmann.me/features/arbitrary_types/#how-can-i-use-get-for-non-default-constructiblenon-copyable-types), rather than the one present in the codebase due to the a) track parameters types not being default constructible, b) use cases connected to the grid-of-of-track-parameters conversion. To keep the converter to one implementation of the json serializer, the liberty was taken to add one more constructor to the `FreeTrackParameters`. The `GridJsonConverter` code was modified a bit to be able to handle non-default constructible types. Header for the track parameters json converter in `GridJsonConverter.hpp` is required for the grid-of-of-track-parameters conversion to be possible by the `nlohmann::json` standard. As a side note, all the types in the framework that have `to_json`/`from_json` overloaded can have grid-of-type json conversion implemented for free, if the respective headers are added to the `GridJsonConverter.hpp`. --------- Co-authored-by: kodiakhq[bot] <49736102+kodiakhq[bot]@users.noreply.github.com> --- .../EventData/GenericFreeTrackParameters.hpp | 28 ++- .../Acts/Plugins/Json/GridJsonConverter.hpp | 11 +- .../Json/TrackParametersJsonConverter.hpp | 221 ++++++++++++++++++ Tests/UnitTests/Plugins/Json/CMakeLists.txt | 1 + .../TrackParametersJsonConverterTests.cpp | 97 ++++++++ 5 files changed, 352 insertions(+), 6 deletions(-) create mode 100644 Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp create mode 100644 Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp diff --git a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp index 857d6c328b9..214aa2e1551 100644 --- a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp +++ b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp @@ -16,6 +16,7 @@ #include "Acts/EventData/detail/PrintParameters.hpp" #include "Acts/Utilities/MathHelpers.hpp" #include "Acts/Utilities/UnitVectors.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" #include #include @@ -56,6 +57,29 @@ class GenericFreeTrackParameters { m_cov(std::move(cov)), m_particleHypothesis(std::move(particleHypothesis)) {} + /// Construct from four-position, direction, absolute momentum, and charge. + /// + /// @param pos4 Track position/time four-vector + /// @param dir Track direction three-vector; normalization is ignored. + /// @param qOverP Charge over momentum + /// @param cov Free parameters covariance matrix + /// @param particleHypothesis Particle hypothesis + GenericFreeTrackParameters(const Vector4& pos4, const Vector3& dir, + Scalar qOverP, std::optional cov, + ParticleHypothesis particleHypothesis) + : m_params(FreeVector::Zero()), + m_cov(std::move(cov)), + m_particleHypothesis(std::move(particleHypothesis)) { + m_params[eFreePos0] = pos4[ePos0]; + m_params[eFreePos1] = pos4[ePos1]; + m_params[eFreePos2] = pos4[ePos2]; + m_params[eFreeTime] = pos4[eTime]; + m_params[eFreeDir0] = dir[eMom0]; + m_params[eFreeDir1] = dir[eMom1]; + m_params[eFreeDir2] = dir[eMom2]; + m_params[eFreeQOverP] = qOverP; + } + /// Construct from four-position, angles, absolute momentum, and charge. /// /// @param pos4 Track position/time four-vector @@ -136,9 +160,9 @@ class GenericFreeTrackParameters { Scalar time() const { return m_params[eFreeTime]; } /// Phi direction. - Scalar phi() const { return phi(direction()); } + Scalar phi() const { return VectorHelpers::phi(direction()); } /// Theta direction. - Scalar theta() const { return theta(direction()); } + Scalar theta() const { return VectorHelpers::theta(direction()); } /// Charge over momentum. Scalar qOverP() const { return m_params[eFreeQOverP]; } diff --git a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp index 20065d488bb..3215bd66eae 100644 --- a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp +++ b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Plugins/Json/ActsJson.hpp" +#include "Acts/Plugins/Json/TrackParametersJsonConverter.hpp" #include "Acts/Utilities/AxisFwd.hpp" #include "Acts/Utilities/GridAccessHelpers.hpp" #include "Acts/Utilities/IAxis.hpp" @@ -266,15 +267,17 @@ auto fromJson(const nlohmann::json& jGrid, if constexpr (GridType::DIM == 1u) { for (const auto& jd : jData) { std::array lbin = jd[0u]; - value_type values = jd[1u]; - grid.atLocalBins(lbin) = values; + if (!jd[1u].is_null()) { + grid.atLocalBins(lbin) = jd[1u].get(); + } } } if constexpr (GridType::DIM == 2u) { for (const auto& jd : jData) { std::array lbin = jd[0u]; - value_type values = jd[1u]; - grid.atLocalBins(lbin) = values; + if (!jd[1u].is_null()) { + grid.atLocalBins(lbin) = jd[1u].get(); + } } } return grid; diff --git a/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp new file mode 100644 index 00000000000..ebf7d5c6054 --- /dev/null +++ b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp @@ -0,0 +1,221 @@ +// 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/EventData/TrackParameters.hpp" +#include "Acts/Plugins/Json/ActsJson.hpp" +#include "Acts/Plugins/Json/SurfaceJsonConverter.hpp" + +#include + +namespace { + +// Alias to bound adl_serializer specialization +// only to track parameters +template +concept TrackParameters = Acts::FreeTrackParametersConcept || + Acts::BoundTrackParametersConcept; + +// Shorthand for bound track parameters +template +concept IsGenericBound = + std::same_as>; + +} // namespace + +namespace Acts { +NLOHMANN_JSON_SERIALIZE_ENUM(Acts::PdgParticle, + + {{Acts::PdgParticle::eInvalid, "Invalid"}, + {Acts::PdgParticle::eElectron, "Electron"}, + {Acts::PdgParticle::eAntiElectron, + "AntiElectron"}, + {Acts::PdgParticle::ePositron, "Positron"}, + {Acts::PdgParticle::eMuon, "Muon"}, + {Acts::PdgParticle::eAntiMuon, "AntiMuon"}, + {Acts::PdgParticle::eTau, "Tau"}, + {Acts::PdgParticle::eAntiTau, "AntiTau"}, + {Acts::PdgParticle::eGamma, "Gamma"}, + {Acts::PdgParticle::ePionZero, "PionZero"}, + {Acts::PdgParticle::ePionPlus, "PionPlus"}, + {Acts::PdgParticle::ePionMinus, "PionMinus"}, + {Acts::PdgParticle::eKaonPlus, "KaonPlus"}, + {Acts::PdgParticle::eKaonMinus, "KaonMinus"}, + {Acts::PdgParticle::eNeutron, "Neutron"}, + {Acts::PdgParticle::eAntiNeutron, "AntiNeutron"}, + {Acts::PdgParticle::eProton, "Proton"}, + {Acts::PdgParticle::eAntiProton, "AntiProton"}, + {Acts::PdgParticle::eLead, "Lead"}} + +) +} + +namespace nlohmann { + +/// @brief Serialize a track parameters object to json +/// +/// nlohmann::json serializer specialized for track parameters +/// as they are not default constructible. Is able to serialize +/// either bound or free track parameters given that the constructor +/// convention is followed. +/// +/// @tparam parameters_t The track parameters type +template +struct adl_serializer { + /// Covariance matrix type attached to the parameters + using CovarianceMatrix = typename parameters_t::CovarianceMatrix; + + /// @brief Serialize track parameters object to json + /// + /// @param j Json object to write to + /// @param t Track parameters object to serialize + static void to_json(nlohmann::json& j, const parameters_t& t) { + // Serialize parameters + // common to all track parameters + j["direction"] = t.direction(); + j["qOverP"] = t.qOverP(); + j["particleHypothesis"] = t.particleHypothesis().absolutePdg(); + + // Covariance is optional + j["covariance"]; + if (t.covariance().has_value()) { + // Extract covariance matrix + // parameters and serialize + auto cov = t.covariance().value(); + constexpr unsigned int size = cov.rows(); + std::array covData{}; + for (std::size_t n = 0; n < size; ++n) { + for (std::size_t m = 0; m < size; ++m) { + covData[n * size + m] = cov(n, m); + } + } + j["covariance"] = covData; + } + // Bound track parameters have + // reference surface attached + // and position takes a geometry context + if constexpr (IsGenericBound) { + Acts::GeometryContext gctx; + j["position"] = t.fourPosition(gctx); + + j["referenceSurface"] = + Acts::SurfaceJsonConverter::toJson(gctx, t.referenceSurface()); + } else { + j["position"] = t.fourPosition(); + } + } + + /// @brief Deserialize track parameters object from json + /// + /// @param j Json object to read from + /// @return Track parameters object + static parameters_t from_json(const nlohmann::json& j) { + // Extract common parameters + std::array posData = j.at("position"); + Acts::Vector4 position(posData[0], posData[1], posData[2], posData[3]); + + std::array dirData = j.at("direction"); + Acts::Vector3 direction(dirData[0], dirData[1], dirData[2]); + + Acts::ActsScalar qOverP = j.at("qOverP"); + Acts::PdgParticle absPdg = j.at("particleHypothesis"); + + // Covariance is optional + std::optional cov; + if (j.at("covariance").is_null()) { + cov = std::nullopt; + } else { + // Extract covariance matrix + // parameters and deserialize + CovarianceMatrix mat; + constexpr unsigned int size = mat.rows(); + std::array covData = j.at("covariance"); + for (std::size_t n = 0; n < size; ++n) { + for (std::size_t m = 0; m < size; ++m) { + mat(n, m) = covData[n * size + m]; + } + } + cov.emplace(std::move(mat)); + } + + // Create particle hypothesis + typename parameters_t::ParticleHypothesis particle(absPdg); + + // Bound track parameters have + // reference surface attached + // and constructor is hidden + // behind a factory method + if constexpr (IsGenericBound) { + Acts::GeometryContext gctx; + auto referenceSurface = + Acts::SurfaceJsonConverter::fromJson(j.at("referenceSurface")); + + auto res = parameters_t::create(referenceSurface, gctx, position, + direction, qOverP, cov, particle); + + if (!res.ok()) { + throw std::invalid_argument("Invalid bound track parameters"); + } + return res.value(); + } else { + return parameters_t(position, direction, qOverP, cov, particle); + } + } +}; + +/// @brief Serialize a shared pointer to track parameters object to json +/// +/// nlohmann::json serializer specialized for shared pointers to track +/// parameters as they are not default constructible. Is able to serialize +/// either bound or free track parameters given that the constructor +/// convention is followed. +/// +/// @tparam parameters_t The track parameters type +template +struct adl_serializer> { + using CovarianceMatrix = typename parameters_t::CovarianceMatrix; + static void to_json(nlohmann::json& j, + const std::shared_ptr& t) { + if (t == nullptr) { + return; + } + j = *t; + } + + static std::shared_ptr from_json(const nlohmann::json& j) { + return std::make_shared(j.get()); + } +}; + +/// @brief Serialize a unique pointer to track parameters object to json +/// +/// nlohmann::json serializer specialized for unique pointers to track +/// parameters as they are not default constructible. Is able to serialize +/// either bound or free track parameters given that the constructor +/// convention is followed. +/// +/// @tparam parameters_t The track parameters type +template +struct adl_serializer> { + using CovarianceMatrix = typename parameters_t::CovarianceMatrix; + static void to_json(nlohmann::json& j, + const std::unique_ptr& t) { + if (t == nullptr) { + return; + } + j = *t; + } + + static std::unique_ptr from_json(const nlohmann::json& j) { + return std::make_unique(j.get()); + } +}; + +} // namespace nlohmann diff --git a/Tests/UnitTests/Plugins/Json/CMakeLists.txt b/Tests/UnitTests/Plugins/Json/CMakeLists.txt index 1f8573fdf50..141f86fb3fe 100644 --- a/Tests/UnitTests/Plugins/Json/CMakeLists.txt +++ b/Tests/UnitTests/Plugins/Json/CMakeLists.txt @@ -15,3 +15,4 @@ add_unittest(UtilitiesJsonConverter UtilitiesJsonConverterTests.cpp) add_unittest(SurfaceBoundsJsonConverter SurfaceBoundsJsonConverterTests.cpp) add_unittest(SurfaceJsonConverter SurfaceJsonConverterTests.cpp) add_unittest(VolumeBoundsJsonConverter VolumeBoundsJsonConverterTests.cpp) +add_unittest(TrackParametersJsonConverter TrackParametersJsonConverterTests.cpp) diff --git a/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp new file mode 100644 index 00000000000..566464a1889 --- /dev/null +++ b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp @@ -0,0 +1,97 @@ +// 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 + +#include "Acts/Plugins/Json/TrackParametersJsonConverter.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" + +#include +#include + +#include + +BOOST_AUTO_TEST_SUITE(TrackParametersJsonIO) + +BOOST_AUTO_TEST_CASE(TrackParametersJsonIO) { + Acts::GeometryContext gctx; + + // Track parameters + Acts::Vector4 position(1., 2., 3., 4.); + Acts::ActsScalar phi = 0.1; + Acts::ActsScalar theta = 0.2; + Acts::ActsScalar qOverP = 3.0; + Acts::ParticleHypothesis particle = Acts::ParticleHypothesis::electron(); + Acts::FreeMatrix freeCov = Acts::FreeMatrix::Identity(); + Acts::BoundMatrix boundCov = Acts::BoundMatrix::Identity(); + + auto surface = Acts::Surface::makeShared( + Acts::Transform3::Identity(), + std::make_shared(10., 10.)); + surface->assignGeometryId(Acts::GeometryIdentifier(1u)); + + // Free track parameters conversion + Acts::FreeTrackParameters ftp(position, phi, theta, qOverP, freeCov, + particle); + + nlohmann::json ftpJson = ftp; + + Acts::FreeTrackParameters ftpRead = ftpJson; + + BOOST_CHECK_EQUAL(ftp.position(), ftpRead.position()); + BOOST_CHECK_EQUAL(ftp.direction(), ftpRead.direction()); + BOOST_CHECK_EQUAL(ftp.qOverP(), ftpRead.qOverP()); + BOOST_CHECK_EQUAL(ftp.covariance().value(), ftpRead.covariance().value()); + BOOST_CHECK_EQUAL(ftp.particleHypothesis(), ftpRead.particleHypothesis()); + + // Curvilinear track parameters conversion + Acts::CurvilinearTrackParameters ctp(position, phi, theta, qOverP, boundCov, + particle); + + nlohmann::json ctpJson = ctp; + + Acts::CurvilinearTrackParameters ctpRead = ctpJson; + + BOOST_CHECK_EQUAL(ctp.position(), ctpRead.position()); + BOOST_CHECK_EQUAL(ctp.direction(), ctpRead.direction()); + BOOST_CHECK_EQUAL(ctp.qOverP(), ctpRead.qOverP()); + BOOST_CHECK_EQUAL(ctp.covariance().value(), ctpRead.covariance().value()); + BOOST_CHECK_EQUAL(ctp.particleHypothesis(), ctpRead.particleHypothesis()); + + BOOST_CHECK(ctp.referenceSurface().transform(gctx).isApprox( + ctpRead.referenceSurface().transform(gctx))); + BOOST_CHECK_EQUAL(ctp.referenceSurface().geometryId(), + ctpRead.referenceSurface().geometryId()); + BOOST_CHECK_EQUAL(ctp.referenceSurface().bounds(), + ctpRead.referenceSurface().bounds()); + + // Bound track parameters conversion + Acts::BoundVector boundPosition{1., 2., 3., 4., 5., 6.}; + Acts::BoundTrackParameters btp(surface, boundPosition, boundCov, particle); + + nlohmann::json btpJson = btp; + + Acts::BoundTrackParameters btpRead = btpJson; + + BOOST_CHECK_EQUAL(btp.position(gctx), btpRead.position(gctx)); + BOOST_CHECK_EQUAL(btp.direction(), btpRead.direction()); + BOOST_CHECK_EQUAL(btp.qOverP(), btpRead.qOverP()); + BOOST_CHECK_EQUAL(btp.covariance().value(), btpRead.covariance().value()); + BOOST_CHECK_EQUAL(btp.particleHypothesis(), btpRead.particleHypothesis()); + + BOOST_CHECK(btp.referenceSurface().transform(gctx).isApprox( + btpRead.referenceSurface().transform(gctx))); + BOOST_CHECK_EQUAL(btp.referenceSurface().geometryId(), + btpRead.referenceSurface().geometryId()); + BOOST_CHECK_EQUAL(btp.referenceSurface().bounds(), + btpRead.referenceSurface().bounds()); +} + +BOOST_AUTO_TEST_SUITE_END() From 41066dc3e90520693fe52baafe494a46e62169db Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Fri, 18 Oct 2024 10:32:29 +0200 Subject: [PATCH 8/8] feat: Make it possible to turn off warnings in GSF (#3739) This adds an option that should disable the warnings appearing in the GSF in a meaningful way... Worksaround but does not solve #3738 --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 31 +++++++++++++++---- Core/src/TrackFitting/BetheHeitlerApprox.cpp | 8 +++-- .../python/acts/examples/reconstruction.py | 7 ++++- Examples/Python/src/TrackFitting.cpp | 18 +++++++---- Examples/Python/tests/root_file_hashes.txt | 8 ++--- Examples/Python/tests/test_examples.py | 2 +- .../Python/truth_tracking_gsf_refitting.py | 7 ++++- 7 files changed, 59 insertions(+), 22 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 653bef7e899..666ebf6aa23 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -125,6 +125,7 @@ class AtlasBetheHeitlerApprox { constexpr static double m_singleGaussianLimit = 0.002; double m_lowLimit = 0.10; double m_highLimit = 0.20; + bool m_clampToRange = false; public: /// Construct the Bethe-Heitler approximation description with two @@ -138,16 +139,19 @@ class AtlasBetheHeitlerApprox { /// @param highTransform whether the high data need to be transformed /// @param lowLimit the upper limit for the low data /// @param highLimit the upper limit for the high data + /// @param clampToRange whether to clamp the input x/x0 to the allowed range constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData, bool lowTransform, bool highTransform, double lowLimit = 0.1, - double highLimit = 0.2) + double highLimit = 0.2, + bool clampToRange = false) : m_lowData(lowData), m_highData(highData), m_lowTransform(lowTransform), m_highTransform(highTransform), m_lowLimit(lowLimit), - m_highLimit(highLimit) {} + m_highLimit(highLimit), + m_clampToRange(clampToRange) {} /// Returns the number of components the returned mixture will have constexpr auto numComponents() const { return NComponents; } @@ -155,7 +159,13 @@ class AtlasBetheHeitlerApprox { /// Checks if an input is valid for the parameterization /// /// @param x pathlength in terms of the radiation length - constexpr bool validXOverX0(ActsScalar x) const { return x < m_highLimit; } + constexpr bool validXOverX0(ActsScalar x) const { + if (m_clampToRange) { + return true; + } else { + return x < m_highLimit; + } + } /// Generates the mixture from the polynomials and reweights them, so /// that the sum of all weights is 1 @@ -164,6 +174,11 @@ class AtlasBetheHeitlerApprox { auto mixture(ActsScalar x) const { using Array = boost::container::static_vector; + + if (m_clampToRange) { + x = std::clamp(x, 0.0, m_highLimit); + } + // Build a polynom auto poly = [](ActsScalar xx, const std::array &coeffs) { @@ -238,9 +253,11 @@ class AtlasBetheHeitlerApprox { /// the parameterization for high x/x0 /// @param lowLimit the upper limit for the low x/x0-data /// @param highLimit the upper limit for the high x/x0-data + /// @param clampToRange forwarded to constructor static auto loadFromFiles(const std::string &low_parameters_path, const std::string &high_parameters_path, - double lowLimit = 0.1, double highLimit = 0.2) { + double lowLimit = 0.1, double highLimit = 0.2, + bool clampToRange = false) { auto read_file = [](const std::string &filepath) { std::ifstream file(filepath); @@ -284,7 +301,8 @@ class AtlasBetheHeitlerApprox { const auto [highData, highTransform] = read_file(high_parameters_path); return AtlasBetheHeitlerApprox(lowData, highData, lowTransform, - highTransform, lowLimit, highLimit); + highTransform, lowLimit, highLimit, + clampToRange); } }; @@ -292,6 +310,7 @@ class AtlasBetheHeitlerApprox { /// configuration, that are stored as static data in the source code. /// This may not be an optimal configuration, but should allow to run /// the GSF without the need to load files -AtlasBetheHeitlerApprox<6, 5> makeDefaultBetheHeitlerApprox(); +AtlasBetheHeitlerApprox<6, 5> makeDefaultBetheHeitlerApprox( + bool clampToRange = false); } // namespace Acts diff --git a/Core/src/TrackFitting/BetheHeitlerApprox.cpp b/Core/src/TrackFitting/BetheHeitlerApprox.cpp index b19b9921902..272d6e08c8d 100644 --- a/Core/src/TrackFitting/BetheHeitlerApprox.cpp +++ b/Core/src/TrackFitting/BetheHeitlerApprox.cpp @@ -8,7 +8,8 @@ #include "Acts/TrackFitting/BetheHeitlerApprox.hpp" -Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() { +Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox( + bool clampToRange) { // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par // clang-format off constexpr static AtlasBetheHeitlerApprox<6, 5>::Data cdf_cmps6_order5_data = {{ @@ -51,6 +52,7 @@ Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() { }}; // clang-format on - return AtlasBetheHeitlerApprox<6, 5>( - cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true, 0.2, 0.2); + return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data, + cdf_cmps6_order5_data, true, true, 0.2, + 0.2, clampToRange); } diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 3bb898ca998..371345a17b5 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1375,8 +1375,13 @@ def addTruthTrackingGsf( ) -> None: customLogLevel = acts.examples.defaultLogging(s, logLevel) + # NOTE we specify clampToRange as True to silence warnings in the test about + # queries to the loss distribution outside the specified range, since no dedicated + # approximation for the ODD is done yet. + bha = acts.examples.AtlasBetheHeitlerApprox.makeDefault(clampToRange=True) + gsfOptions = { - "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(), + "betheHeitlerApprox": bha, "maxComponents": 12, "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight, "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance, diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 20c6010c8ee..c8dc72f9eef 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -39,6 +39,7 @@ namespace py = pybind11; using namespace ActsExamples; using namespace Acts; +using namespace py::literals; namespace Acts::Python { @@ -106,12 +107,17 @@ void addTrackFitting(Context& ctx) { .value("KLDistance", MixtureReductionAlgorithm::KLDistance); py::class_(mex, "AtlasBetheHeitlerApprox") - .def_static("loadFromFiles", - &ActsExamples::BetheHeitlerApprox::loadFromFiles, - py::arg("lowParametersPath"), py::arg("highParametersPath"), - py::arg("lowLimit") = 0.1, py::arg("highLimit") = 0.2) - .def_static("makeDefault", - []() { return Acts::makeDefaultBetheHeitlerApprox(); }); + .def_static( + "loadFromFiles", &ActsExamples::BetheHeitlerApprox::loadFromFiles, + "lowParametersPath"_a, "highParametersPath"_a, "lowLimit"_a = 0.1, + "highLimit"_a = 0.2, "clampToRange"_a = false) + .def_static( + "makeDefault", + [](bool clampToRange) { + return Acts::makeDefaultBetheHeitlerApprox(clampToRange); + }, + "clampToRange"_a = false); + mex.def( "makeGsfFitterFunction", [](std::shared_ptr trackingGeometry, diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 39df454b3e8..54c81ea9ba1 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -74,10 +74,10 @@ test_exatrkx[cpu-torch]__performance_track_finding.root: 36b3045589c4c17c038dbc8 test_exatrkx[gpu-onnx]__performance_track_finding.root: 9090de10ffb1489d3f1993e2a3081a3038227e3e5c453e98a9a4f33ea3d6d817 test_exatrkx[gpu-torch]__performance_track_finding.root: 36b3045589c4c17c038dbc87943366f4af4440f7eea6887afb763871ac149b05 test_ML_Ambiguity_Solver__performance_ambiML.root: 284ff5c3a08c0b810938e4ac2f8ba8fe2babb17d4c202b624ed69fff731a9006 -test_refitting[odd]__trackstates_gsf_refit.root: a61fe2d80d5d10d3b3505000e8abb589d88302bf5f54b625948793418f2a7fe8 -test_refitting[odd]__tracksummary_gsf_refit.root: 082789fc1a85e578b3cf9a750723d2dcee01b5019e871c6a63e0b271f4e907b1 -test_refitting[generic]__trackstates_gsf_refit.root: 9fa7af9eff12081504c0d648f6debae64f6996e0ca610cf58187d23aa5a13251 -test_refitting[generic]__tracksummary_gsf_refit.root: 35b5ac6f208cae093fff94038d217a2e9915a5ce075da2a95718dda696f2d4a2 +test_refitting[odd]__trackstates_gsf_refit.root: e297749dc1e7eda3b8dea13defa0499986c584740d93e723a901b498b8e90c71 +test_refitting[odd]__tracksummary_gsf_refit.root: d5085882e45a0b699194dff9f40a36e9291227bf65f9aaaf9087f9242ef5ae22 +test_refitting[generic]__trackstates_gsf_refit.root: 4424fdf2f27575db825c1a59f8e53a1595946211cbd5b2c8d3a2f71cdcc77ae9 +test_refitting[generic]__tracksummary_gsf_refit.root: 562deecee4cfb97ceee72eff53d63da079e3249fb62d6bcd556e6f27d495dfd9 test_truth_tracking_kalman[generic-False-0.0]__trackstates_kf.root: 9f77962b92037cb760b1629a602b1dae61f45e659c45d9a87baa784f6190960e test_truth_tracking_kalman[generic-False-0.0]__tracksummary_kf.root: 562deecee4cfb97ceee72eff53d63da079e3249fb62d6bcd556e6f27d495dfd9 test_truth_tracking_kalman[generic-False-1000.0]__trackstates_kf.root: 56a1bd989b9c1316b9098c65fa75df9e6683e62e35ae68d8f72d27220be0fd7d diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py index 0e34a8a3c1c..d5bb593ca9f 100644 --- a/Examples/Python/tests/test_examples.py +++ b/Examples/Python/tests/test_examples.py @@ -689,7 +689,7 @@ def test_refitting(tmp_path, detector_config, assert_root_hash): field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) seq = Sequencer( - events=3, + events=10, numThreads=1, ) diff --git a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py index 797bfc59569..e8b41c9aae8 100755 --- a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py +++ b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py @@ -25,8 +25,13 @@ def runRefittingGsf( s=s, ) + # NOTE we specify clampToRange as True to silence warnings in the test about + # queries to the loss distribution outside the specified range, since no dedicated + # approximation for the ODD is done yet. + bha = acts.examples.AtlasBetheHeitlerApprox.makeDefault(clampToRange=True) + gsfOptions = { - "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(), + "betheHeitlerApprox": bha, "maxComponents": 12, "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight, "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance,