From 2eee79491fcb9e48ab0ca90731f421b664c6e0a6 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 29 Nov 2024 09:52:57 +0100 Subject: [PATCH 1/8] refactor!: Clean track EDM projector (#3605) Clean deprecated projector functions from `TrackStateProxy` and rename the new methods to something recognizable. Also deals with the ACTS internal downstream changes. ## Summary by CodeRabbit - **New Features** - Enhanced handling of projector subspace indices across various components for improved accuracy in alignment and calibration processes. - **Bug Fixes** - Updated method calls to ensure consistent use of `setProjectorSubspaceIndices`, enhancing the clarity and functionality of the track state management. - **Documentation** - Minor adjustments to comments for better clarity on the changes related to projector handling. - **Tests** - Updated test cases to reflect changes in method names and ensure proper validation of projector subspace indices. --- .../Kernel/detail/AlignmentEngine.hpp | 4 +- .../Acts/EventData/SubspaceHelpers.hpp | 30 +-- .../Acts/EventData/TrackStateProxy.hpp | 171 ++++-------------- .../Acts/EventData/TrackStateProxy.ipp | 5 - .../Acts/EventData/TrackStateProxyConcept.hpp | 6 - .../detail/MultiTrajectoryTestsCommon.hpp | 31 +--- .../Acts/EventData/detail/TestSourceLink.hpp | 22 +-- .../Acts/TrackFinding/MeasurementSelector.ipp | 2 +- .../Acts/TrackFitting/GainMatrixSmoother.hpp | 1 - .../Acts/TrackFitting/GainMatrixUpdater.hpp | 2 +- .../TrackFitting/GlobalChiSquareFitter.hpp | 2 +- .../include/Acts/TrackFitting/MbfSmoother.hpp | 7 +- .../detail/GainMatrixUpdaterImpl.hpp | 2 +- .../Acts/TrackFitting/detail/GsfUtils.hpp | 8 +- Core/include/Acts/Utilities/TrackHelpers.hpp | 7 +- Core/src/TrackFinding/MeasurementSelector.cpp | 5 +- Core/src/TrackFitting/GsfUtils.cpp | 13 +- Core/src/TrackFitting/MbfSmoother.cpp | 16 +- .../TrackFitting/src/RefittingCalibrator.cpp | 2 +- .../Framework/ML/src/NeuralCalibrator.cpp | 2 +- .../src/EventData/MeasurementCalibration.cpp | 2 +- .../src/EventData/ScalingCalibrator.cpp | 2 +- .../Io/Root/src/RootTrackStatesWriter.cpp | 8 +- Tests/Benchmarks/TrackEdmBenchmark.cpp | 2 +- .../Core/EventData/TrackTestsExtra.cpp | 2 +- .../Core/TrackFitting/FitterTestsCommon.hpp | 11 +- .../Core/TrackFitting/MbfSmootherTests.cpp | 9 +- 27 files changed, 127 insertions(+), 247 deletions(-) diff --git a/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp b/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp index 1bb21fdc8bd..b26175abcdd 100644 --- a/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp +++ b/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp @@ -208,7 +208,9 @@ TrackAlignmentState trackAlignmentState( measdim) = measCovariance; // (b) Get and fill the bound parameters to measurement projection matrix - const ActsDynamicMatrix H = state.effectiveProjector(); + const ActsDynamicMatrix H = + state.projectorSubspaceHelper().fullProjector().topLeftCorner( + measdim, eBoundSize); alignState.projectionMatrix.block(iMeasurement, iParams, measdim, eBoundSize) = H; // (c) Get and fill the residual diff --git a/Core/include/Acts/EventData/SubspaceHelpers.hpp b/Core/include/Acts/EventData/SubspaceHelpers.hpp index 029dbf372b3..f48690cf036 100644 --- a/Core/include/Acts/EventData/SubspaceHelpers.hpp +++ b/Core/include/Acts/EventData/SubspaceHelpers.hpp @@ -14,8 +14,7 @@ #include "Acts/Utilities/AlgebraHelpers.hpp" #include "Acts/Utilities/Enumerate.hpp" -#include -#include +#include #include @@ -25,30 +24,30 @@ namespace Acts { /// /// Indices must be unique and within the full size of the subspace /// -/// @tparam Container type of the container +/// @tparam index_range_t the type of the container of indices /// -/// @param container the container of indices +/// @param indexRange the range of indices /// @param fullSize the full size of the subspace /// @param subspaceSize the size of the subspace /// /// @return true if the indices are consistent -template -inline static bool checkSubspaceIndices(const Container& container, +template +inline static bool checkSubspaceIndices(const index_range_t& indexRange, std::size_t fullSize, std::size_t subspaceSize) { if (subspaceSize > fullSize) { return false; } - if (static_cast(container.size()) != subspaceSize) { + if (static_cast(indexRange.size()) != subspaceSize) { return false; } - for (auto it = container.begin(); it != container.end();) { + for (auto it = indexRange.begin(); it != indexRange.end();) { auto index = *it; if (index >= fullSize) { return false; } ++it; - if (std::find(it, container.end(), index) != container.end()) { + if (std::find(it, indexRange.end(), index) != indexRange.end()) { return false; } } @@ -69,7 +68,8 @@ inline static SerializedSubspaceIndices serializeSubspaceIndices( { SerializedSubspaceIndices result = 0; for (std::size_t i = 0; i < FullSize; ++i) { - result |= static_cast(indices[i]) << (i * 8); + result |= static_cast(indices[i] & 0xFF) + << (i * 8); } return result; } @@ -88,7 +88,7 @@ inline static SubspaceIndices deserializeSubspaceIndices( { SubspaceIndices result; for (std::size_t i = 0; i < FullSize; ++i) { - result[i] = static_cast(serialized >> (i * 8)); + result[i] = static_cast((serialized >> (i * 8)) & 0xFF); } return result; } @@ -187,8 +187,8 @@ class VariableSubspaceHelper using IndexType = index_t; using Container = boost::container::static_vector; - template - explicit VariableSubspaceHelper(const OtherContainer& indices) { + template + explicit VariableSubspaceHelper(const other_index_range_t& indices) { assert(checkSubspaceIndices(indices, kFullSize, indices.size()) && "Invalid indices"); m_indices.resize(indices.size()); @@ -236,8 +236,8 @@ class FixedSubspaceHelper using IndexType = index_t; using Container = std::array; - template - explicit FixedSubspaceHelper(const OtherContainer& indices) { + template + explicit FixedSubspaceHelper(const other_index_range_t& indices) { assert(checkSubspaceIndices(indices, kFullSize, kSubspaceSize) && "Invalid indices"); std::transform(indices.begin(), indices.end(), m_indices.begin(), diff --git a/Core/include/Acts/EventData/TrackStateProxy.hpp b/Core/include/Acts/EventData/TrackStateProxy.hpp index 3f7ddf1a028..4338c434f77 100644 --- a/Core/include/Acts/EventData/TrackStateProxy.hpp +++ b/Core/include/Acts/EventData/TrackStateProxy.hpp @@ -8,7 +8,6 @@ #pragma once -#include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/SubspaceHelpers.hpp" @@ -17,11 +16,11 @@ #include "Acts/EventData/TrackStateType.hpp" #include "Acts/EventData/Types.hpp" #include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/AlgebraHelpers.hpp" #include "Acts/Utilities/HashedString.hpp" #include "Acts/Utilities/Helpers.hpp" #include +#include #include #include @@ -127,11 +126,6 @@ struct TrackStateTraits { typename detail_lt::DynamicSizeTypes::CoefficientsMap; using EffectiveCalibratedCovariance = typename detail_lt::DynamicSizeTypes::CovarianceMap; - - constexpr static auto ProjectorFlags = Eigen::RowMajor | Eigen::AutoAlign; - using Projector = Eigen::Matrix; - using EffectiveProjector = Eigen::Matrix; }; /// Proxy object to access a single point on the trajectory. @@ -207,17 +201,6 @@ class TrackStateProxy { /// Sentinel value that indicates an invalid index static constexpr IndexType kInvalid = kTrackIndexInvalid; - /// Matrix representing the projector (measurement mapping function) for a - /// measurement. This is not a map type, but an actual matrix. This matrix - /// is always \f$M \times M\f$, even if the local measurement dimension is lower. - /// The actual \f$N\times M\f$ projector is given by the top \f$N\f$ rows. - using Projector = typename TrackStateTraits::Projector; - - /// Dynamic variant of the projector matrix - /// @warning Using this type is discouraged, as it has a runtime overhead - using EffectiveProjector = - typename TrackStateTraits::EffectiveProjector; - /// The track state container backend given as a template parameter using Trajectory = trajectory_t; @@ -615,145 +598,61 @@ class TrackStateProxy { /// /// @{ - /// Returns the projector (measurement mapping function) for this track - /// state. It is derived from the uncalibrated measurement - /// @note This function returns the overallocated projector. This means it - /// is of dimension MxM, where M is the maximum number of measurement - /// dimensions. The NxM submatrix, where N is the actual dimension of the - /// measurement, is located in the top left corner, everything else is zero. - /// @return The overallocated projector - Projector projector() const; - - /// Returns whether a projector is set - /// @return Whether it is set - bool hasProjector() const { return has(); } - - /// Returns the projector (measurement mapping function) for this track - /// state. It is derived from the uncalibrated measurement - /// @warning This function returns the effective projector. This means it - /// is of dimension \f$N\times M\f$, where \f$N\f$ is the actual dimension of the - /// measurement. - /// @return The effective projector - EffectiveProjector effectiveProjector() const { - return projector().topLeftCorner(calibratedSize(), M); - } - - /// Set the projector on this track state - /// This will convert the projector to a more compact bitset representation - /// and store it. - /// @param projector The projector in the form of a dense matrix - /// @note @p projector is assumed to only have 0s or 1s as components. - template - [[deprecated("use setProjector(span) instead")]] void setProjector( - const Eigen::MatrixBase& projector) - requires(!ReadOnly) + /// Set the projector subspace indices + /// @param subspaceIndices The projector subspace indices to set + template + void setProjectorSubspaceIndices(const index_range_t& subspaceIndices) + requires(!ReadOnly && + std::convertible_to, + std::uint8_t>) { - constexpr int rows = Eigen::MatrixBase::RowsAtCompileTime; - constexpr int cols = Eigen::MatrixBase::ColsAtCompileTime; - - static_assert(rows != -1 && cols != -1, - "Assignment of dynamic matrices is currently not supported."); - assert(has()); - - static_assert(rows <= M, "Given projector has too many rows"); - static_assert(cols <= eBoundSize, "Given projector has too many columns"); - - // set up full size projector with only zeros - typename TrackStateProxy::Projector fullProjector = - decltype(fullProjector)::Zero(); - - // assign (potentially) smaller actual projector to matrix, preserving - // zeroes outside of smaller matrix block. - fullProjector.template topLeftCorner() = projector; - - // convert to bitset before storing - ProjectorBitset projectorBitset = matrixToBitset(fullProjector).to_ulong(); - setProjectorBitset(projectorBitset); - } - - /// Directly get the projector bitset, a compressed form of a projection - /// matrix - /// @note This is mainly to copy explicitly a projector from one state - /// to another. Use the `projector` or `effectiveProjector` method if - /// you want to access the matrix. - /// @return The projector bitset - [[deprecated("use projector() instead")]] ProjectorBitset projectorBitset() - const { - return variableBoundSubspaceHelper().projectorBitset(); + assert(subspaceIndices.size() <= eBoundSize); + BoundSubspaceIndices boundSubspace{}; + std::transform(subspaceIndices.begin(), subspaceIndices.end(), + boundSubspace.begin(), + [](auto i) { return static_cast(i); }); + component() = + serializeSubspaceIndices(boundSubspace); } - /// Set the projector bitset, a compressed form of a projection matrix - /// @param proj The projector bitset - /// - /// @note This is mainly to copy explicitly a projector from one state - /// to another. If you have a projection matrix, set it with - /// `setProjector`. - [[deprecated("use setProjector(span) instead")]] void setProjectorBitset( - ProjectorBitset proj) - requires(!ReadOnly) - { - BoundMatrix projMatrix = bitsetToMatrix(proj); - BoundSubspaceIndices boundSubspace = - projectorToSubspaceIndices(projMatrix); - setBoundSubspaceIndices(boundSubspace); - } + /// Returns whether a projector is set + /// @return Whether it is set + bool hasProjector() const { return has(); } - BoundSubspaceIndices boundSubspaceIndices() const { + /// Returns the projector subspace indices + /// @return The projector subspace indices + BoundSubspaceIndices projectorSubspaceIndices() const { assert(has()); return deserializeSubspaceIndices( component()); } + /// Returns the projector subspace indices + /// @return The projector subspace indices template - SubspaceIndices subspaceIndices() const { - BoundSubspaceIndices boundSubspace = BoundSubspaceIndices(); + SubspaceIndices projectorSubspaceIndices() const { + BoundSubspaceIndices boundSubspace = projectorSubspaceIndices(); SubspaceIndices subspace; std::copy(boundSubspace.begin(), boundSubspace.begin() + measdim, subspace.begin()); return subspace; } - void setBoundSubspaceIndices(BoundSubspaceIndices boundSubspace) - requires(!ReadOnly) - { - assert(has()); - component() = - serializeSubspaceIndices(boundSubspace); - } - - template - void setSubspaceIndices(SubspaceIndices subspace) - requires(!ReadOnly && measdim <= eBoundSize) - { - assert(has()); - BoundSubspaceIndices boundSubspace{}; - std::copy(subspace.begin(), subspace.end(), boundSubspace.begin()); - setBoundSubspaceIndices(boundSubspace); - } - - template - void setSubspaceIndices(std::array subspaceIndices) - requires(!ReadOnly && measdim <= eBoundSize) - { - assert(has()); - BoundSubspaceIndices boundSubspace{}; - std::transform(subspaceIndices.begin(), subspaceIndices.end(), - boundSubspace.begin(), - [](index_t i) { return static_cast(i); }); - setBoundSubspaceIndices(boundSubspace); - } - - VariableBoundSubspaceHelper variableBoundSubspaceHelper() const { - BoundSubspaceIndices boundSubspace = boundSubspaceIndices(); + /// Creates a variable size subspace helper + /// @return The subspace helper + VariableBoundSubspaceHelper projectorSubspaceHelper() const { + BoundSubspaceIndices boundSubspace = projectorSubspaceIndices(); std::span validSubspaceIndices( boundSubspace.begin(), boundSubspace.begin() + calibratedSize()); return VariableBoundSubspaceHelper(validSubspaceIndices); } + /// Creates a fixed size subspace helper + /// @return The subspace helper template - FixedBoundSubspaceHelper fixedBoundSubspaceHelper() const { - SubspaceIndices subspace = subspaceIndices(); + FixedBoundSubspaceHelper projectorSubspaceHelper() const { + SubspaceIndices subspace = projectorSubspaceIndices(); return FixedBoundSubspaceHelper(subspace); } @@ -1028,7 +927,7 @@ class TrackStateProxy { other.template calibratedCovariance().eval()); }); - setBoundSubspaceIndices(other.boundSubspaceIndices()); + setProjectorSubspaceIndices(other.projectorSubspaceIndices()); } } else { if (ACTS_CHECK_BIT(mask, PM::Predicted) && @@ -1073,7 +972,7 @@ class TrackStateProxy { other.template calibratedCovariance().eval()); }); - setBoundSubspaceIndices(other.boundSubspaceIndices()); + setProjectorSubspaceIndices(other.projectorSubspaceIndices()); } } diff --git a/Core/include/Acts/EventData/TrackStateProxy.ipp b/Core/include/Acts/EventData/TrackStateProxy.ipp index e36be5f3c25..16abe8fb3aa 100644 --- a/Core/include/Acts/EventData/TrackStateProxy.ipp +++ b/Core/include/Acts/EventData/TrackStateProxy.ipp @@ -61,11 +61,6 @@ inline auto TrackStateProxy::covariance() const } } -template -inline auto TrackStateProxy::projector() const -> Projector { - return variableBoundSubspaceHelper().fullProjector(); -} - template inline auto TrackStateProxy::getUncalibratedSourceLink() const -> SourceLink { diff --git a/Core/include/Acts/EventData/TrackStateProxyConcept.hpp b/Core/include/Acts/EventData/TrackStateProxyConcept.hpp index 61e98adbc48..4b8a713f864 100644 --- a/Core/include/Acts/EventData/TrackStateProxyConcept.hpp +++ b/Core/include/Acts/EventData/TrackStateProxyConcept.hpp @@ -111,12 +111,6 @@ concept TrackStateProxyConcept = { cv.hasProjector() } -> std::same_as; { v.hasProjector() } -> std::same_as; - { cv.effectiveProjector() } -> std::same_as; - { v.effectiveProjector() } -> std::same_as; - - { cv.projectorBitset() } -> std::same_as; - { v.projectorBitset() } -> std::same_as; - { cv.getUncalibratedSourceLink() } -> std::same_as; { v.getUncalibratedSourceLink() } -> std::same_as; diff --git a/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp b/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp index 368138fe38d..7926ca3ba77 100644 --- a/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp +++ b/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp @@ -471,8 +471,8 @@ class MultiTrajectoryTestsCommon { pc.sourceLink.sourceId); // Explicitly unset to avoid error below ts.unset(TrackStatePropMask::Calibrated); - testSourceLinkCalibratorReturn( - gctx, cctx, SourceLink{ttsb.sourceLink}, ts); + testSourceLinkCalibrator(gctx, cctx, + SourceLink{ttsb.sourceLink}, ts); BOOST_CHECK_EQUAL( ts.getUncalibratedSourceLink().template get().sourceId, ttsb.sourceLink.sourceId); @@ -549,21 +549,6 @@ class MultiTrajectoryTestsCommon { BOOST_CHECK_EQUAL(ts.template calibratedCovariance(), expCov); }); } - - BOOST_CHECK(ts.hasProjector()); - ActsMatrix fullProj; - fullProj.setZero(); - { - Acts::GeometryContext gctx; - Acts::CalibrationContext cctx; - // create a temporary measurement to extract the projector matrix - testSourceLinkCalibratorReturn( - gctx, cctx, SourceLink{pc.sourceLink}, ts); - fullProj = ts.projector(); - } - BOOST_CHECK_EQUAL(ts.effectiveProjector(), - fullProj.topLeftCorner(nMeasurements, eBoundSize)); - BOOST_CHECK_EQUAL(ts.projector(), fullProj); } void testTrackStateProxyAllocations(std::default_random_engine& rng) { @@ -782,7 +767,8 @@ class MultiTrajectoryTestsCommon { }); BOOST_CHECK_NE(ts1.calibratedSize(), ts2.calibratedSize()); - BOOST_CHECK_NE(ts1.projector(), ts2.projector()); + BOOST_CHECK(ts1.projectorSubspaceIndices() != + ts2.projectorSubspaceIndices()); BOOST_CHECK_NE(ts1.jacobian(), ts2.jacobian()); BOOST_CHECK_NE(ts1.chi2(), ts2.chi2()); @@ -813,7 +799,8 @@ class MultiTrajectoryTestsCommon { }); BOOST_CHECK_EQUAL(ts1.calibratedSize(), ts2.calibratedSize()); - BOOST_CHECK_EQUAL(ts1.projector(), ts2.projector()); + BOOST_CHECK(ts1.projectorSubspaceIndices() == + ts2.projectorSubspaceIndices()); BOOST_CHECK_EQUAL(ts1.jacobian(), ts2.jacobian()); BOOST_CHECK_EQUAL(ts1.chi2(), ts2.chi2()); @@ -840,7 +827,8 @@ class MultiTrajectoryTestsCommon { }); BOOST_CHECK_NE(ts1.calibratedSize(), ts2.calibratedSize()); - BOOST_CHECK_NE(ts1.projector(), ts2.projector()); + BOOST_CHECK(ts1.projectorSubspaceIndices() != + ts2.projectorSubspaceIndices()); BOOST_CHECK_NE(ts1.jacobian(), ts2.jacobian()); BOOST_CHECK_NE(ts1.chi2(), ts2.chi2()); @@ -864,7 +852,8 @@ class MultiTrajectoryTestsCommon { }); BOOST_CHECK_EQUAL(ts1.calibratedSize(), ts2.calibratedSize()); - BOOST_CHECK_EQUAL(ts1.projector(), ts2.projector()); + BOOST_CHECK(ts1.projectorSubspaceIndices() == + ts2.projectorSubspaceIndices()); BOOST_CHECK_EQUAL(ts1.jacobian(), ts2.jacobian()); BOOST_CHECK_EQUAL(ts1.chi2(), ts2.chi2()); // always copied diff --git a/Core/include/Acts/EventData/detail/TestSourceLink.hpp b/Core/include/Acts/EventData/detail/TestSourceLink.hpp index 4a98b1a28a9..9a6350e7df4 100644 --- a/Core/include/Acts/EventData/detail/TestSourceLink.hpp +++ b/Core/include/Acts/EventData/detail/TestSourceLink.hpp @@ -18,9 +18,7 @@ #include "Acts/Geometry/TrackingGeometry.hpp" #include "Acts/Utilities/CalibrationContext.hpp" -#include #include -#include #include #include #include @@ -126,7 +124,7 @@ inline std::ostream& operator<<(std::ostream& os, /// @param gctx Unused /// @param trackState TrackState to calibrated template -void testSourceLinkCalibratorReturn( +void testSourceLinkCalibrator( const GeometryContext& /*gctx*/, const CalibrationContext& /*cctx*/, const SourceLink& sourceLink, typename trajectory_t::TrackStateProxy trackState) { @@ -139,30 +137,18 @@ void testSourceLinkCalibratorReturn( trackState.allocateCalibrated(2); trackState.template calibrated<2>() = sl.parameters; trackState.template calibratedCovariance<2>() = sl.covariance; - trackState.setSubspaceIndices(std::array{sl.indices[0], sl.indices[1]}); + trackState.setProjectorSubspaceIndices( + 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.setSubspaceIndices(std::array{sl.indices[0]}); + trackState.setProjectorSubspaceIndices(std::array{sl.indices[0]}); } else { throw std::runtime_error( "Tried to extract measurement from invalid TestSourceLink"); } } -/// Extract the measurement from a TestSourceLink. -/// -/// @param gctx Unused -/// @param trackState TrackState to calibrated -template -void testSourceLinkCalibrator( - const GeometryContext& gctx, const CalibrationContext& cctx, - const SourceLink& sourceLink, - typename trajectory_t::TrackStateProxy trackState) { - testSourceLinkCalibratorReturn(gctx, cctx, sourceLink, - trackState); -} - } // namespace Acts::detail::Test diff --git a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp index 30d09c002be..5e13ffd904a 100644 --- a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp +++ b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp @@ -63,7 +63,7 @@ MeasurementSelector::select( trackState.effectiveCalibrated().data(), trackState.effectiveCalibratedCovariance().data(), trackState.predicted(), trackState.predictedCovariance(), - trackState.boundSubspaceIndices(), trackState.calibratedSize()); + trackState.projectorSubspaceIndices(), trackState.calibratedSize()); trackState.chi2() = chi2; if (chi2 < minChi2) { diff --git a/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp b/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp index cf5b372f487..5adc8340cfe 100644 --- a/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp +++ b/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp @@ -10,7 +10,6 @@ #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/TrackFitting/KalmanFitterError.hpp" #include "Acts/Utilities/Delegate.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" diff --git a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp index 0ac45a79b1d..9d8ae5a0d03 100644 --- a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp +++ b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp @@ -80,7 +80,7 @@ class GainMatrixUpdater { // shape later trackState.effectiveCalibrated().data(), trackState.effectiveCalibratedCovariance().data(), - trackState.boundSubspaceIndices(), + trackState.projectorSubspaceIndices(), trackState.predicted(), trackState.predictedCovariance(), trackState.filtered(), diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index eaa9020053d..dc0bc221e9a 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -423,7 +423,7 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, trackState.template calibrated(); const ActsMatrix projector = - trackState.projector().template topLeftCorner(); + trackState.template projectorSubspaceHelper().projector(); const Eigen::MatrixXd projJacobian = projector * extendedJacobian; diff --git a/Core/include/Acts/TrackFitting/MbfSmoother.hpp b/Core/include/Acts/TrackFitting/MbfSmoother.hpp index 02066a642e7..be7a8c12311 100644 --- a/Core/include/Acts/TrackFitting/MbfSmoother.hpp +++ b/Core/include/Acts/TrackFitting/MbfSmoother.hpp @@ -87,9 +87,6 @@ class MbfSmoother { /// Internal track state representation for the smoother. /// @note This allows us to move parts of the implementation into the .cpp struct InternalTrackState final { - using Projector = - typename TrackStateTraits::Projector; using Jacobian = typename TrackStateTraits::Covariance; @@ -105,14 +102,14 @@ class MbfSmoother { // This is used to build a covariance matrix view in the .cpp file const double* calibrated{nullptr}; const double* calibratedCovariance{nullptr}; - Projector projector; + BoundSubspaceIndices projector; template explicit Measurement(TrackStateProxy ts) : calibratedSize(ts.calibratedSize()), calibrated(ts.effectiveCalibrated().data()), calibratedCovariance(ts.effectiveCalibratedCovariance().data()), - projector(ts.projector()) {} + projector(ts.projectorSubspaceIndices()) {} }; Jacobian jacobian; diff --git a/Core/include/Acts/TrackFitting/detail/GainMatrixUpdaterImpl.hpp b/Core/include/Acts/TrackFitting/detail/GainMatrixUpdaterImpl.hpp index 0dce3c2ba21..0437ca551a3 100644 --- a/Core/include/Acts/TrackFitting/detail/GainMatrixUpdaterImpl.hpp +++ b/Core/include/Acts/TrackFitting/detail/GainMatrixUpdaterImpl.hpp @@ -35,7 +35,7 @@ std::tuple GainMatrixUpdater::visitMeasurementImpl( ACTS_VERBOSE("Calibrated measurement: " << calibrated.transpose()); ACTS_VERBOSE("Calibrated measurement covariance:\n" << calibratedCovariance); - std::span validSubspaceIndices( + std::span validSubspaceIndices( trackState.projector.begin(), trackState.projector.begin() + kMeasurementSize); FixedBoundSubspaceHelper subspaceHelper( diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp index 98605052bf6..574a69e8f71 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp @@ -13,6 +13,7 @@ #include "Acts/EventData/MultiComponentTrackParameters.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/Types.hpp" #include "Acts/Utilities/Logger.hpp" #include @@ -154,9 +155,7 @@ double calculateDeterminant( const double *fullCalibratedCovariance, TrackStateTraits::Covariance predictedCovariance, - TrackStateTraits::Projector - projector, - unsigned int calibratedSize); + BoundSubspaceIndices projector, unsigned int calibratedSize); /// Reweight the components according to `R. Frühwirth, "Track fitting /// with non-Gaussian noise"`. See also the implementation in Athena at @@ -190,7 +189,8 @@ void computePosteriorWeights( .template calibratedCovariance< MultiTrajectoryTraits::MeasurementSizeMax>() .data(), - state.predictedCovariance(), state.projector(), state.calibratedSize()); + state.predictedCovariance(), state.projectorSubspaceIndices(), + state.calibratedSize()); const auto factor = std::sqrt(1. / detR) * safeExp(-0.5 * chi2); diff --git a/Core/include/Acts/Utilities/TrackHelpers.hpp b/Core/include/Acts/Utilities/TrackHelpers.hpp index 2e86a62910d..2e40c6b4988 100644 --- a/Core/include/Acts/Utilities/TrackHelpers.hpp +++ b/Core/include/Acts/Utilities/TrackHelpers.hpp @@ -693,8 +693,11 @@ std::pair calculateUnbiasedParametersCovariance( return visit_measurement( trackState.calibratedSize(), [&](std::integral_constant) { - auto H = trackState.projector() - .template topLeftCorner(); + FixedBoundSubspaceHelper subspaceHelper = + trackState.template projectorSubspaceHelper(); + + // TODO use subspace helper for projection instead + auto H = subspaceHelper.projector(); auto s = trackState.smoothed(); auto C = trackState.smoothedCovariance(); auto m = trackState.template calibrated(); diff --git a/Core/src/TrackFinding/MeasurementSelector.cpp b/Core/src/TrackFinding/MeasurementSelector.cpp index 2a3c0e2d390..676cf8074ff 100644 --- a/Core/src/TrackFinding/MeasurementSelector.cpp +++ b/Core/src/TrackFinding/MeasurementSelector.cpp @@ -9,7 +9,6 @@ #include "Acts/TrackFinding/MeasurementSelector.hpp" #include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/MeasurementHelpers.hpp" #include "Acts/EventData/SubspaceHelpers.hpp" #include "Acts/EventData/Types.hpp" @@ -19,6 +18,8 @@ #include #include #include +#include +#include namespace Acts { @@ -99,7 +100,7 @@ double MeasurementSelector::calculateChi2( using ParametersVector = ActsVector; - std::span validSubspaceIndices( + std::span validSubspaceIndices( projector.begin(), projector.begin() + kMeasurementSize); FixedBoundSubspaceHelper subspaceHelper( validSubspaceIndices); diff --git a/Core/src/TrackFitting/GsfUtils.cpp b/Core/src/TrackFitting/GsfUtils.cpp index b70544e47c6..8ff52f2c1b6 100644 --- a/Core/src/TrackFitting/GsfUtils.cpp +++ b/Core/src/TrackFitting/GsfUtils.cpp @@ -9,8 +9,12 @@ #include "Acts/TrackFitting/detail/GsfUtils.hpp" #include "Acts/EventData/MeasurementHelpers.hpp" +#include "Acts/EventData/SubspaceHelpers.hpp" +#include "Acts/EventData/Types.hpp" #include +#include +#include namespace Acts::detail { @@ -19,17 +23,20 @@ using TrackStateTraits = double calculateDeterminant(const double* fullCalibratedCovariance, TrackStateTraits::Covariance predictedCovariance, - TrackStateTraits::Projector projector, + BoundSubspaceIndices projector, unsigned int calibratedSize) { return visit_measurement(calibratedSize, [&](auto N) { constexpr std::size_t kMeasurementSize = decltype(N)::value; + std::span validSubspaceIndices( + projector.begin(), projector.begin() + kMeasurementSize); + FixedBoundSubspaceHelper subspaceHelper( + validSubspaceIndices); typename Acts::TrackStateTraits< kMeasurementSize, true>::CalibratedCovariance calibratedCovariance{ fullCalibratedCovariance}; - const auto H = - projector.template topLeftCorner().eval(); + const auto H = subspaceHelper.projector(); return (H * predictedCovariance * H.transpose() + calibratedCovariance) .determinant(); diff --git a/Core/src/TrackFitting/MbfSmoother.cpp b/Core/src/TrackFitting/MbfSmoother.cpp index e260a578e38..9b300989e90 100644 --- a/Core/src/TrackFitting/MbfSmoother.cpp +++ b/Core/src/TrackFitting/MbfSmoother.cpp @@ -10,6 +10,8 @@ #include "Acts/EventData/TrackParameterHelpers.hpp" +#include + namespace Acts { void MbfSmoother::calculateSmoothed(InternalTrackState& ts, @@ -42,9 +44,13 @@ void MbfSmoother::visitMeasurement(const InternalTrackState& ts, visit_measurement(measurement.calibratedSize, [&](auto N) -> void { constexpr std::size_t kMeasurementSize = decltype(N)::value; + std::span validSubspaceIndices( + measurement.projector.begin(), + measurement.projector.begin() + kMeasurementSize); + FixedBoundSubspaceHelper subspaceHelper( + validSubspaceIndices); - using MeasurementMatrix = - Eigen::Matrix; + using ProjectorMatrix = Eigen::Matrix; using CovarianceMatrix = Eigen::Matrix; using KalmanGainMatrix = @@ -55,10 +61,8 @@ void MbfSmoother::visitMeasurement(const InternalTrackState& ts, typename TrackStateTraits::CalibratedCovariance calibratedCovariance{measurement.calibratedCovariance}; - // Measurement matrix - const MeasurementMatrix H = - measurement.projector - .template topLeftCorner(); + // Projector matrix + const ProjectorMatrix H = subspaceHelper.projector(); // Residual covariance const CovarianceMatrix S = diff --git a/Examples/Algorithms/TrackFitting/src/RefittingCalibrator.cpp b/Examples/Algorithms/TrackFitting/src/RefittingCalibrator.cpp index eaf2f74bffa..154dc5d7285 100644 --- a/Examples/Algorithms/TrackFitting/src/RefittingCalibrator.cpp +++ b/Examples/Algorithms/TrackFitting/src/RefittingCalibrator.cpp @@ -35,7 +35,7 @@ void RefittingCalibrator::calibrate(const Acts::GeometryContext& /*gctx*/, sl.state.template calibratedCovariance().eval()); }); - trackState.setBoundSubspaceIndices(sl.state.boundSubspaceIndices()); + trackState.setProjectorSubspaceIndices(sl.state.projectorSubspaceIndices()); } } // namespace ActsExamples diff --git a/Examples/Framework/ML/src/NeuralCalibrator.cpp b/Examples/Framework/ML/src/NeuralCalibrator.cpp index 83e4d11354e..83c7a2bc09a 100644 --- a/Examples/Framework/ML/src/NeuralCalibrator.cpp +++ b/Examples/Framework/ML/src/NeuralCalibrator.cpp @@ -190,6 +190,6 @@ void ActsExamples::NeuralCalibrator::calibrate( calibratedCovariance(boundLoc1, boundLoc1) = output[iVar0 + 1]; trackState.allocateCalibrated(calibratedParameters, calibratedCovariance); - trackState.setSubspaceIndices(fixedMeasurement.subspaceIndices()); + trackState.setProjectorSubspaceIndices(fixedMeasurement.subspaceIndices()); }); } diff --git a/Examples/Framework/src/EventData/MeasurementCalibration.cpp b/Examples/Framework/src/EventData/MeasurementCalibration.cpp index 190d37153bb..d7b8c96da70 100644 --- a/Examples/Framework/src/EventData/MeasurementCalibration.cpp +++ b/Examples/Framework/src/EventData/MeasurementCalibration.cpp @@ -43,7 +43,7 @@ void ActsExamples::PassThroughCalibrator::calibrate( trackState.allocateCalibrated(fixedMeasurement.parameters().eval(), fixedMeasurement.covariance().eval()); - trackState.setSubspaceIndices(fixedMeasurement.subspaceIndices()); + trackState.setProjectorSubspaceIndices(fixedMeasurement.subspaceIndices()); }); } diff --git a/Examples/Framework/src/EventData/ScalingCalibrator.cpp b/Examples/Framework/src/EventData/ScalingCalibrator.cpp index 2a956902513..50e693e8b15 100644 --- a/Examples/Framework/src/EventData/ScalingCalibrator.cpp +++ b/Examples/Framework/src/EventData/ScalingCalibrator.cpp @@ -179,6 +179,6 @@ void ActsExamples::ScalingCalibrator::calibrate( calibratedCovariance(boundLoc1, boundLoc1) *= ct.y_scale; trackState.allocateCalibrated(calibratedParameters, calibratedCovariance); - trackState.setSubspaceIndices(fixedMeasurement.subspaceIndices()); + trackState.setProjectorSubspaceIndices(fixedMeasurement.subspaceIndices()); }); } diff --git a/Examples/Io/Root/src/RootTrackStatesWriter.cpp b/Examples/Io/Root/src/RootTrackStatesWriter.cpp index f90fa96ab0e..37d7cdf2fc9 100644 --- a/Examples/Io/Root/src/RootTrackStatesWriter.cpp +++ b/Examples/Io/Root/src/RootTrackStatesWriter.cpp @@ -452,8 +452,8 @@ ProcessCode RootTrackStatesWriter::writeT(const AlgorithmContext& ctx, m_t_eT.push_back(static_cast(truthParams[Acts::eBoundTime])); // expand the local measurements into the full bound space - Acts::BoundVector meas = state.effectiveProjector().transpose() * - state.effectiveCalibrated(); + Acts::BoundVector meas = state.projectorSubspaceHelper().expandVector( + state.effectiveCalibrated()); // extract local and global position Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]); Acts::Vector3 global = @@ -633,7 +633,9 @@ ProcessCode RootTrackStatesWriter::writeT(const AlgorithmContext& ctx, if (ipar == ePredicted) { // local hit residual info - auto H = state.effectiveProjector(); + auto H = + state.projectorSubspaceHelper().fullProjector().topLeftCorner( + state.calibratedSize(), Acts::eBoundSize); auto V = state.effectiveCalibratedCovariance(); auto resCov = V + H * covariance * H.transpose(); Acts::ActsDynamicVector res = diff --git a/Tests/Benchmarks/TrackEdmBenchmark.cpp b/Tests/Benchmarks/TrackEdmBenchmark.cpp index a85b0482be2..ce14eef8e8a 100644 --- a/Tests/Benchmarks/TrackEdmBenchmark.cpp +++ b/Tests/Benchmarks/TrackEdmBenchmark.cpp @@ -149,7 +149,7 @@ int main(int /*argc*/, char** /*argv[]*/) { std::array indices{0}; std::iota(indices.begin(), indices.end(), 0); - trackState.setBoundSubspaceIndices(indices); + trackState.setProjectorSubspaceIndices(indices); }); trackState.typeFlags().set(TrackStateFlag::MeasurementFlag); diff --git a/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp b/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp index f2f27ed8522..ecf49d7cd5b 100644 --- a/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp +++ b/Tests/UnitTests/Core/EventData/TrackTestsExtra.cpp @@ -457,7 +457,7 @@ BOOST_AUTO_TEST_CASE(CopyTrackProxyCalibrated) { auto track1 = tc.makeTrack(); auto ts = track1.appendTrackState(TrackStatePropMask::Calibrated); ts.allocateCalibrated(kMeasurementSize); - ts.setSubspaceIndices(BoundSubspaceIndices{}); + ts.setProjectorSubspaceIndices(BoundSubspaceIndices{}); auto tsCopy = track1.appendTrackState(TrackStatePropMask::Calibrated); tsCopy.copyFrom(ts, TrackStatePropMask::Calibrated, false); diff --git a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp index 269370d4fc4..dca8ac045ca 100644 --- a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp +++ b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp @@ -10,6 +10,7 @@ #include +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/ProxyAccessor.hpp" @@ -56,9 +57,13 @@ struct TestOutlierFinder { if (!state.hasCalibrated() || !state.hasPredicted()) { return false; } - auto residuals = (state.effectiveCalibrated() - - state.effectiveProjector() * state.predicted()) - .eval(); + auto subspaceHelper = state.projectorSubspaceHelper(); + auto projector = + subspaceHelper.fullProjector() + .topLeftCorner(state.calibratedSize(), Acts::eBoundSize) + .eval(); + auto residuals = + (state.effectiveCalibrated() - projector * state.predicted()).eval(); auto distance = residuals.norm(); return (distanceMax <= distance); } diff --git a/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp b/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp index 21d0a347e56..e87d386b5b5 100644 --- a/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp @@ -8,9 +8,7 @@ #include -#include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/TrackStatePropMask.hpp" #include "Acts/EventData/TrackStateType.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" @@ -19,7 +17,6 @@ #include "Acts/TrackFitting/MbfSmoother.hpp" #include "Acts/Utilities/Result.hpp" -#include #include #include @@ -59,7 +56,7 @@ BOOST_AUTO_TEST_CASE(Smooth) { ts.allocateCalibrated(2); ts.calibrated<2>() << 0.351, 0.473; ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8; - ts.setSubspaceIndices<2>(projector); + ts.setProjectorSubspaceIndices(projector); ts.filtered() << 0.301, 0.503, std::numbers::pi / 2., 0., 1 / 100., 0.; ts.filteredCovariance() = covTrk; @@ -76,7 +73,7 @@ BOOST_AUTO_TEST_CASE(Smooth) { ts.allocateCalibrated(2); ts.calibrated<2>() << 0.351, 0.473; ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8; - ts.setSubspaceIndices<2>(projector); + ts.setProjectorSubspaceIndices(projector); ts.filtered() << 0.27, 0.53, std::numbers::pi / 2., 0., 1 / 100., 0.; ts.filteredCovariance() = covTrk; @@ -93,7 +90,7 @@ BOOST_AUTO_TEST_CASE(Smooth) { ts.allocateCalibrated(2); ts.calibrated<2>() << 0.351, 0.473; ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8; - ts.setSubspaceIndices<2>(projector); + ts.setProjectorSubspaceIndices(projector); ts.filtered() << 0.33, 0.43, std::numbers::pi / 2., 0., 1 / 100., 0.; ts.filteredCovariance() = covTrk; From 9c9444e35b0fd71374ce0fa4c9ff3c55854dc976 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 29 Nov 2024 11:17:57 +0100 Subject: [PATCH 2/8] refactor: Replace `ParticleSmearing` with `TrackParameterSmearing` in Examples (#3784) - extract track parameters from particles - replace `ParticleSmearing` with `TrackParameterSmearing` ## Summary by CodeRabbit - **New Features** - Introduced `TrackParameterSmearing` and `ParticleTrackParamExtractor` classes for enhanced track parameter processing. - Updated smearing configurations and parameters to improve clarity and usability. - Enhanced Python bindings for various generator classes, allowing for customizable event generation. - **Bug Fixes** - Adjusted logging levels for new algorithms to ensure appropriate verbosity. - **Documentation** - Updated comments and documentation to reflect changes in parameter names and functionalities. - **Chores** - Removed deprecated `ParticleSmearing` algorithm and associated files from the project. --- .../workflows/physmon_trackfinding_1muon.py | 18 +- .../Generators/EventGenerator.cpp | 21 ++- .../Generators/EventGenerator.hpp | 4 +- .../TruthTracking/ParticleSmearing.cpp | 158 ---------------- .../ParticleTrackParamExtractor.cpp | 72 +++++++ .../ParticleTrackParamExtractor.hpp | 48 +++++ .../TruthTracking/TrackParameterSmearing.cpp | 175 ++++++++++++++++++ ...mearing.hpp => TrackParameterSmearing.hpp} | 49 ++--- .../Algorithms/TruthTracking/CMakeLists.txt | 3 +- .../python/acts/examples/reconstruction.py | 67 ++++--- Examples/Python/src/Generators.cpp | 2 - Examples/Python/src/TruthTracking.cpp | 17 +- Examples/Python/tests/conftest.py | 19 +- Examples/Python/tests/root_file_hashes.txt | 12 +- Examples/Python/tests/test_algorithms.py | 4 +- Examples/Python/tests/test_propagation.py | 25 +-- Examples/Scripts/Optimization/ckf.py | 18 +- Examples/Scripts/Python/ckf_tracks.py | 18 +- Examples/Scripts/Python/full_chain_test.py | 4 +- .../Scripts/Python/material_validation.py | 16 +- .../Scripts/Python/material_validation_itk.py | 22 +-- Examples/Scripts/Python/propagation.py | 17 +- Examples/Scripts/Python/vertex_fitting.py | 13 +- 23 files changed, 471 insertions(+), 331 deletions(-) delete mode 100644 Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp create mode 100644 Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.cpp create mode 100644 Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp create mode 100644 Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.cpp rename Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/{ParticleSmearing.hpp => TrackParameterSmearing.hpp} (61%) diff --git a/CI/physmon/workflows/physmon_trackfinding_1muon.py b/CI/physmon/workflows/physmon_trackfinding_1muon.py index aaa4bc9f71d..90305112099 100755 --- a/CI/physmon/workflows/physmon_trackfinding_1muon.py +++ b/CI/physmon/workflows/physmon_trackfinding_1muon.py @@ -18,7 +18,7 @@ from acts.examples.reconstruction import ( addSeeding, - ParticleSmearingSigmas, + TrackSmearingSigmas, SeedFinderConfigArg, SeedFinderOptionsArg, SeedingAlgorithm, @@ -91,15 +91,15 @@ def run_ckf_tracking(label, seeding): s, setup.trackingGeometry, setup.field, - ParticleSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared + TrackSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared # zero eveything so the CKF has a chance to find the measurements - d0=0, - d0PtA=0, - d0PtB=0, - z0=0, - z0PtA=0, - z0PtB=0, - t0=0, + loc0=0, + loc0PtA=0, + loc0PtB=0, + loc1=0, + loc1PtA=0, + loc1PtB=0, + time=0, phi=0, theta=0, ptRel=0, diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.cpp b/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.cpp index a11d62c928f..9d5c1b873d3 100644 --- a/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.cpp +++ b/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.cpp @@ -8,17 +8,23 @@ #include "ActsExamples/Generators/EventGenerator.hpp" +#include "Acts/Surfaces/PerigeeSurface.hpp" #include "ActsExamples/EventData/SimVertex.hpp" +#include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" #include "ActsFatras/EventData/Barcode.hpp" #include "ActsFatras/EventData/Particle.hpp" #include +#include +#include #include #include +#include -ActsExamples::EventGenerator::EventGenerator(const Config& cfg, - Acts::Logging::Level lvl) +namespace ActsExamples { + +EventGenerator::EventGenerator(const Config& cfg, Acts::Logging::Level lvl) : m_cfg(cfg), m_logger(Acts::getDefaultLogger("EventGenerator", lvl)) { if (m_cfg.outputParticles.empty()) { throw std::invalid_argument("Missing output particles collection"); @@ -37,17 +43,15 @@ ActsExamples::EventGenerator::EventGenerator(const Config& cfg, m_outputVertices.initialize(m_cfg.outputVertices); } -std::string ActsExamples::EventGenerator::name() const { +std::string EventGenerator::name() const { return "EventGenerator"; } -std::pair -ActsExamples::EventGenerator::availableEvents() const { +std::pair EventGenerator::availableEvents() const { return {0u, std::numeric_limits::max()}; } -ActsExamples::ProcessCode ActsExamples::EventGenerator::read( - const AlgorithmContext& ctx) { +ProcessCode EventGenerator::read(const AlgorithmContext& ctx) { SimParticleContainer particles; SimVertexContainer vertices; @@ -120,5 +124,8 @@ ActsExamples::ProcessCode ActsExamples::EventGenerator::read( // move generated event to the store m_outputParticles(ctx, std::move(particles)); m_outputVertices(ctx, std::move(vertices)); + return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.hpp b/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.hpp index 46720b9fff1..7a76ac03d71 100644 --- a/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.hpp +++ b/Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.hpp @@ -18,7 +18,6 @@ #include "ActsExamples/Framework/RandomNumbers.hpp" #include -#include #include #include #include @@ -93,8 +92,9 @@ class EventGenerator final : public ActsExamples::IReader { struct Config { /// Name of the output particles collection. std::string outputParticles; - /// Name of the vertex collection. + /// Name of the output vertex collection. std::string outputVertices; + /// List of generators that should be used to generate the event. std::vector generators; /// The random number service. diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp deleted file mode 100644 index 4c09134ef20..00000000000 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp +++ /dev/null @@ -1,158 +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/. - -#include "ActsExamples/TruthTracking/ParticleSmearing.hpp" - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" -#include "Acts/Surfaces/PerigeeSurface.hpp" -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/detail/periodic.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/EventData/Track.hpp" -#include "ActsExamples/Framework/AlgorithmContext.hpp" -#include "ActsExamples/Framework/RandomNumbers.hpp" - -#include -#include -#include -#include -#include - -ActsExamples::ParticleSmearing::ParticleSmearing(const Config& config, - Acts::Logging::Level level) - : IAlgorithm("ParticleSmearing", level), m_cfg(config) { - if (m_cfg.inputParticles.empty()) { - throw std::invalid_argument("Missing input truth particles collection"); - } - if (m_cfg.outputTrackParameters.empty()) { - throw std::invalid_argument("Missing output tracks parameters collection"); - } - if (m_cfg.randomNumbers == nullptr) { - throw std::invalid_argument("Missing random numbers tool"); - } - - if (m_cfg.particleHypothesis) { - ACTS_INFO("Override truth particle hypothesis with " - << *m_cfg.particleHypothesis); - } - - m_inputParticles.initialize(m_cfg.inputParticles); - m_outputTrackParameters.initialize(m_cfg.outputTrackParameters); -} - -ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute( - const AlgorithmContext& ctx) const { - // setup input and output containers - const auto& particles = m_inputParticles(ctx); - TrackParametersContainer parameters; - parameters.reserve(particles.size()); - - // setup random number generator and standard gaussian - auto rng = m_cfg.randomNumbers->spawnGenerator(ctx); - std::normal_distribution stdNormal(0.0, 1.0); - - for (auto&& [vtxId, vtxParticles] : groupBySecondaryVertex(particles)) { - // a group contains at least one particle by construction. assume that all - // particles within the group originate from the same position and use it to - // as the reference position for the perigee frame. - auto perigee = Acts::Surface::makeShared( - vtxParticles.begin()->position()); - - for (const auto& particle : vtxParticles) { - const auto time = particle.time(); - const auto phi = Acts::VectorHelpers::phi(particle.direction()); - const auto theta = Acts::VectorHelpers::theta(particle.direction()); - const auto pt = particle.transverseMomentum(); - const auto p = particle.absoluteMomentum(); - const auto qOverP = particle.qOverP(); - const auto particleHypothesis = - m_cfg.particleHypothesis.value_or(particle.hypothesis()); - - // compute momentum-dependent resolutions - const double sigmaD0 = - m_cfg.sigmaD0 + - m_cfg.sigmaD0PtA * std::exp(-1.0 * std::abs(m_cfg.sigmaD0PtB) * pt); - const double sigmaZ0 = - m_cfg.sigmaZ0 + - m_cfg.sigmaZ0PtA * std::exp(-1.0 * std::abs(m_cfg.sigmaZ0PtB) * pt); - // shortcuts for other resolutions - const double sigmaT0 = m_cfg.sigmaT0; - const double sigmaPhi = m_cfg.sigmaPhi; - const double sigmaTheta = m_cfg.sigmaTheta; - const double sigmaQOverP = - std::sqrt(std::pow(m_cfg.sigmaPtRel * qOverP, 2) + - std::pow(sigmaTheta * (qOverP * std::tan(theta)), 2)); - - Acts::BoundVector params = Acts::BoundVector::Zero(); - // smear the position/time - // note that we smear d0 and z0 in the perigee frame - params[Acts::eBoundLoc0] = sigmaD0 * stdNormal(rng); - params[Acts::eBoundLoc1] = sigmaZ0 * stdNormal(rng); - params[Acts::eBoundTime] = time + sigmaT0 * stdNormal(rng); - // smear direction angles phi,theta ensuring correct bounds - const auto [newPhi, newTheta] = Acts::detail::normalizePhiTheta( - phi + sigmaPhi * stdNormal(rng), theta + sigmaTheta * stdNormal(rng)); - params[Acts::eBoundPhi] = newPhi; - params[Acts::eBoundTheta] = newTheta; - // compute smeared q/p - params[Acts::eBoundQOverP] = qOverP + sigmaQOverP * stdNormal(rng); - - ACTS_VERBOSE("Smearing particle (pos, time, phi, theta, q/p):"); - ACTS_VERBOSE(" from: " << particle.position().transpose() << ", " << time - << ", " << phi << ", " << theta << ", " << qOverP); - ACTS_VERBOSE(" to: " << perigee - ->localToGlobal( - ctx.geoContext, - Acts::Vector2{params[Acts::eBoundLoc0], - params[Acts::eBoundLoc1]}, - particle.direction() * p) - .transpose() - << ", " << params[Acts::eBoundTime] << ", " - << params[Acts::eBoundPhi] << ", " - << params[Acts::eBoundTheta] << ", " - << params[Acts::eBoundQOverP]); - - // build the track covariance matrix using the smearing sigmas - Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero(); - if (m_cfg.initialSigmas) { - // use the initial sigmas if set - - Acts::EstimateTrackParamCovarianceConfig config{ - .initialSigmas = - Eigen::Map{ - m_cfg.initialSigmas->data()}, - .initialSigmaPtRel = m_cfg.initialSigmaPtRel, - .initialVarInflation = Eigen::Map{ - m_cfg.initialVarInflation.data()}}; - - cov = Acts::estimateTrackParamCovariance(config, params, false); - } else { - // otherwise use the smearing sigmas - - Acts::BoundVector sigmas = Acts::BoundVector( - {sigmaD0, sigmaZ0, sigmaPhi, sigmaTheta, sigmaQOverP, sigmaT0}); - - for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) { - double sigma = sigmas[i]; - double variance = sigma * sigma; - - // Inflate the initial covariance - variance *= m_cfg.initialVarInflation[i]; - - cov(i, i) = variance; - } - } - parameters.emplace_back(perigee, params, cov, particleHypothesis); - } - } - - m_outputTrackParameters(ctx, std::move(parameters)); - return ProcessCode::SUCCESS; -} diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.cpp new file mode 100644 index 00000000000..d480ca8f88f --- /dev/null +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.cpp @@ -0,0 +1,72 @@ +// 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 "ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp" + +#include "Acts/Surfaces/PerigeeSurface.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/Framework/AlgorithmContext.hpp" + +#include +#include + +namespace ActsExamples { + +ParticleTrackParamExtractor::ParticleTrackParamExtractor( + const Config& config, Acts::Logging::Level level) + : IAlgorithm("ParticleTrackParamExtractor", level), m_cfg(config) { + if (m_cfg.inputParticles.empty()) { + throw std::invalid_argument("Missing input particles collection"); + } + if (m_cfg.outputTrackParameters.empty()) { + throw std::invalid_argument("Missing output track parameters collection"); + } + + m_inputParticles.initialize(m_cfg.inputParticles); + m_outputTrackParameters.initialize(m_cfg.outputTrackParameters); +} + +ActsExamples::ProcessCode ParticleTrackParamExtractor::execute( + const AlgorithmContext& ctx) const { + const SimParticleContainer& particles = m_inputParticles(ctx); + + std::unordered_map> + perigeeSurfaces; + + for (auto&& [vtxId, vtxParticles] : groupBySecondaryVertex(particles)) { + // a group contains at least one particle by construction. assume that all + // particles within the group originate from the same position and use it + // to as the reference position for the perigee frame. + auto perigee = Acts::Surface::makeShared( + vtxParticles.begin()->position()); + perigeeSurfaces[vtxId] = perigee; + } + + // create track parameters from the particles + TrackParametersContainer trackParameters; + + for (const auto& particle : particles) { + const auto vtxId = particle.particleId().vertexId(); + const auto particleHypothesis = particle.hypothesis(); + const auto phi = Acts::VectorHelpers::phi(particle.direction()); + const auto theta = Acts::VectorHelpers::theta(particle.direction()); + const auto qOverP = particle.qOverP(); + const auto time = particle.time(); + + trackParameters.emplace_back( + perigeeSurfaces.at(vtxId), + Acts::BoundVector{0, 0, phi, theta, qOverP, time}, std::nullopt, + particleHypothesis); + } + + m_outputTrackParameters(ctx, std::move(trackParameters)); + + return ProcessCode::SUCCESS; +} + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp new file mode 100644 index 00000000000..60882fa9ace --- /dev/null +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp @@ -0,0 +1,48 @@ +// 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/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" + +#include + +namespace ActsExamples { +struct AlgorithmContext; + +/// Extract track parameters from particles. +class ParticleTrackParamExtractor final : public IAlgorithm { + public: + struct Config { + /// The input particles collection. + std::string inputParticles; + /// The output track parameters collection. + std::string outputTrackParameters; + }; + + ParticleTrackParamExtractor(const Config& config, Acts::Logging::Level level); + + ProcessCode execute(const AlgorithmContext& ctx) const final; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + Config m_cfg; + + ReadDataHandle m_inputParticles{this, "InputParticles"}; + WriteDataHandle m_outputTrackParameters{ + this, "OutputTrackParameters"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.cpp new file mode 100644 index 00000000000..f1789749f77 --- /dev/null +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.cpp @@ -0,0 +1,175 @@ +// 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 "ActsExamples/TruthTracking/TrackParameterSmearing.hpp" + +#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Framework/RandomNumbers.hpp" + +#include +#include +#include + +namespace ActsExamples { + +TrackParameterSmearing::TrackParameterSmearing(const Config& config, + Acts::Logging::Level level) + : IAlgorithm("TrackParameterSmearing", level), m_cfg(config) { + if (m_cfg.inputTrackParameters.empty()) { + throw std::invalid_argument("Missing input track parameters collection"); + } + if (m_cfg.outputTrackParameters.empty()) { + throw std::invalid_argument("Missing output track parameters collection"); + } + if (m_cfg.randomNumbers == nullptr) { + throw std::invalid_argument("Missing random numbers tool"); + } + + if (m_cfg.particleHypothesis) { + ACTS_INFO("Override truth particle hypothesis with " + << *m_cfg.particleHypothesis); + } + + m_inputTrackParameters.initialize(m_cfg.inputTrackParameters); + m_outputTrackParameters.initialize(m_cfg.outputTrackParameters); + + ACTS_DEBUG("smearing track param loc0 " << m_cfg.sigmaLoc0 << " A " + << m_cfg.sigmaLoc0PtA << " B " + << m_cfg.sigmaLoc0PtB); + ACTS_DEBUG("smearing track param loc1 " << m_cfg.sigmaLoc1 << " A " + << m_cfg.sigmaLoc1PtA << " B " + << m_cfg.sigmaLoc1PtB); + ACTS_DEBUG("smearing track param time " << m_cfg.sigmaTime); + ACTS_DEBUG("smearing track param phi " << m_cfg.sigmaPhi); + ACTS_DEBUG("smearing track param theta " << m_cfg.sigmaTheta); + ACTS_DEBUG("smearing track param q/p " << m_cfg.sigmaPtRel); + ACTS_DEBUG( + "initial sigmas " + << Acts::BoundVector( + m_cfg.initialSigmas.value_or(std::array()).data()) + .transpose()); + ACTS_DEBUG("initial sigma pt rel " << m_cfg.initialSigmaPtRel); + ACTS_DEBUG( + "initial var inflation " + << Acts::BoundVector(m_cfg.initialVarInflation.data()).transpose()); + if (m_cfg.particleHypothesis) { + ACTS_DEBUG("particle hypothesis " << *m_cfg.particleHypothesis); + } else { + ACTS_DEBUG("particle hypothesis truth"); + } +} + +ProcessCode TrackParameterSmearing::execute(const AlgorithmContext& ctx) const { + // setup input and output containers + const auto& inputTrackParametersContainer = m_inputTrackParameters(ctx); + + ACTS_VERBOSE("Smearing " << inputTrackParametersContainer.size() + << " track parameters"); + + TrackParametersContainer outputTrackParametersContainer; + outputTrackParametersContainer.reserve(inputTrackParametersContainer.size()); + + // setup random number generator and standard gaussian + auto rng = m_cfg.randomNumbers->spawnGenerator(ctx); + std::normal_distribution stdNormal(0.0, 1.0); + + for (const auto& inputTrackParameters : inputTrackParametersContainer) { + const auto position = inputTrackParameters.localPosition(); + const auto time = inputTrackParameters.time(); + const auto phi = inputTrackParameters.phi(); + const auto theta = inputTrackParameters.theta(); + const auto pt = inputTrackParameters.transverseMomentum(); + const auto qOverP = inputTrackParameters.qOverP(); + const auto particleHypothesis = m_cfg.particleHypothesis.value_or( + inputTrackParameters.particleHypothesis()); + + // compute momentum-dependent resolutions + const double sigmaLoc0 = + m_cfg.sigmaLoc0 + + m_cfg.sigmaLoc0PtA * std::exp(-1.0 * std::abs(m_cfg.sigmaLoc0PtB) * pt); + const double sigmaLoc1 = + m_cfg.sigmaLoc1 + + m_cfg.sigmaLoc1PtA * std::exp(-1.0 * std::abs(m_cfg.sigmaLoc1PtB) * pt); + // shortcuts for other resolutions + const double sigmaTime = m_cfg.sigmaTime; + const double sigmaPhi = m_cfg.sigmaPhi; + const double sigmaTheta = m_cfg.sigmaTheta; + const double sigmaQOverP = + std::sqrt(std::pow(m_cfg.sigmaPtRel * qOverP, 2) + + std::pow(sigmaTheta * (qOverP * std::tan(theta)), 2)); + + Acts::BoundVector params = Acts::BoundVector::Zero(); + // smear the position/time + // note that we smear d0 and z0 in the perigee frame + params[Acts::eBoundLoc0] = position[0] + sigmaLoc0 * stdNormal(rng); + params[Acts::eBoundLoc1] = position[1] + sigmaLoc1 * stdNormal(rng); + params[Acts::eBoundTime] = time + sigmaTime * stdNormal(rng); + // smear direction angles phi,theta ensuring correct bounds + const auto [newPhi, newTheta] = Acts::detail::normalizePhiTheta( + phi + sigmaPhi * stdNormal(rng), theta + sigmaTheta * stdNormal(rng)); + params[Acts::eBoundPhi] = newPhi; + params[Acts::eBoundTheta] = newTheta; + // compute smeared q/p + params[Acts::eBoundQOverP] = qOverP + sigmaQOverP * stdNormal(rng); + + // build the track covariance matrix using the smearing sigmas + Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero(); + if (m_cfg.initialSigmas) { + // use the initial sigmas if set + + Acts::EstimateTrackParamCovarianceConfig config{ + .initialSigmas = + Eigen::Map{m_cfg.initialSigmas->data()}, + .initialSigmaPtRel = m_cfg.initialSigmaPtRel, + .initialVarInflation = Eigen::Map{ + m_cfg.initialVarInflation.data()}}; + + cov = Acts::estimateTrackParamCovariance(config, params, false); + } else { + // otherwise use the smearing sigmas + + Acts::BoundVector sigmas = Acts::BoundVector( + {sigmaLoc0, sigmaLoc1, sigmaPhi, sigmaTheta, sigmaQOverP, sigmaTime}); + + for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) { + double sigma = sigmas[i]; + double variance = sigma * sigma; + + // Inflate the initial covariance + variance *= m_cfg.initialVarInflation[i]; + + cov(i, i) = variance; + } + } + + const auto& outputTrackParameters = + outputTrackParametersContainer.emplace_back( + inputTrackParameters.referenceSurface().shared_from_this(), params, + cov, particleHypothesis); + + ACTS_VERBOSE("Smearing particle (pos, time, phi, theta, q/p):"); + ACTS_VERBOSE( + " from: " << inputTrackParameters.position(ctx.geoContext).transpose() + << ", " << time << ", " << phi << ", " << theta << ", " + << qOverP); + ACTS_VERBOSE( + " to: " << outputTrackParameters.position(ctx.geoContext).transpose() + << ", " << params[Acts::eBoundTime] << ", " + << params[Acts::eBoundPhi] << ", " + << params[Acts::eBoundTheta] << ", " + << params[Acts::eBoundQOverP]); + } + + m_outputTrackParameters(ctx, std::move(outputTrackParametersContainer)); + + return ProcessCode::SUCCESS; +} + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.hpp similarity index 61% rename from Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp rename to Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.hpp index 3b60f12b099..53bab3b62e5 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.hpp @@ -11,7 +11,6 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IAlgorithm.hpp" @@ -19,7 +18,6 @@ #include "ActsExamples/Framework/RandomNumbers.hpp" #include -#include #include #include #include @@ -28,36 +26,38 @@ namespace ActsExamples { class RandomNumbers; struct AlgorithmContext; -/// Create track states by smearing truth particle information. +/// @brief Smear track parameters. /// -/// Particles are smeared in the perigee frame anchored at their true vertex -/// position. The `d0` and `z0` parameters are always defined within that -/// perigee frame and not globally. The generated bound parameters are stored in -/// the same order as the input particles. -class ParticleSmearing final : public IAlgorithm { +/// Track parameters are smeared in the local frame. The `loc0` and `loc1` +/// parameters are always defined within that local frame and not globally. The +/// generated bound parameters are stored in the same order as the input +/// parameters. +class TrackParameterSmearing final : public IAlgorithm { public: struct Config { - /// Input truth particles collection. - std::string inputParticles; - /// Output smeared tracks parameters collection. + /// Input track parameters collection. + std::string inputTrackParameters; + /// Output smeared track parameters collection. std::string outputTrackParameters; /// Random numbers service. std::shared_ptr randomNumbers = nullptr; // Smearing parameters - /// Constant term of the d0 resolution. - double sigmaD0 = 20 * Acts::UnitConstants::um; - /// Pt-dependent d0 resolution of the form sigma_d0 = A*exp(-1.*abs(B)*pt). - double sigmaD0PtA = 30 * Acts::UnitConstants::um; - double sigmaD0PtB = 0.3 / Acts::UnitConstants::GeV; - /// Constant term of the z0 resolution. - double sigmaZ0 = 20 * Acts::UnitConstants::um; - /// Pt-dependent z0 resolution of the form sigma_z0 = A*exp(-1.*abs(B)*pt). - double sigmaZ0PtA = 30 * Acts::UnitConstants::um; - double sigmaZ0PtB = 0.3 / Acts::UnitConstants::GeV; + /// Constant term of the loc0 resolution. + double sigmaLoc0 = 20 * Acts::UnitConstants::um; + /// Pt-dependent loc0 resolution of the form sigma_loc0 = + /// A*exp(-1.*abs(B)*pt). + double sigmaLoc0PtA = 30 * Acts::UnitConstants::um; + double sigmaLoc0PtB = 0.3 / Acts::UnitConstants::GeV; + /// Constant term of the loc1 resolution. + double sigmaLoc1 = 20 * Acts::UnitConstants::um; + /// Pt-dependent loc1 resolution of the form sigma_loc1 = + /// A*exp(-1.*abs(B)*pt). + double sigmaLoc1PtA = 30 * Acts::UnitConstants::um; + double sigmaLoc1PtB = 0.3 / Acts::UnitConstants::GeV; /// Time resolution. - double sigmaT0 = 1 * Acts::UnitConstants::ns; + double sigmaTime = 1 * Acts::UnitConstants::ns; /// Phi angular resolution. double sigmaPhi = 1 * Acts::UnitConstants::degree; /// Theta angular resolution. @@ -77,7 +77,7 @@ class ParticleSmearing final : public IAlgorithm { std::optional particleHypothesis = std::nullopt; }; - ParticleSmearing(const Config& config, Acts::Logging::Level level); + TrackParameterSmearing(const Config& config, Acts::Logging::Level level); ProcessCode execute(const AlgorithmContext& ctx) const override; @@ -87,7 +87,8 @@ class ParticleSmearing final : public IAlgorithm { private: Config m_cfg; - ReadDataHandle m_inputParticles{this, "InputParticles"}; + ReadDataHandle m_inputTrackParameters{ + this, "InputTrackParameters"}; WriteDataHandle m_outputTrackParameters{ this, "OutputTrackParameters"}; diff --git a/Examples/Algorithms/TruthTracking/CMakeLists.txt b/Examples/Algorithms/TruthTracking/CMakeLists.txt index 280e1104ce2..9f576efa1b4 100644 --- a/Examples/Algorithms/TruthTracking/CMakeLists.txt +++ b/Examples/Algorithms/TruthTracking/CMakeLists.txt @@ -2,7 +2,8 @@ add_library( ActsExamplesTruthTracking SHARED ActsExamples/TruthTracking/ParticleSelector.cpp - ActsExamples/TruthTracking/ParticleSmearing.cpp + ActsExamples/TruthTracking/ParticleTrackParamExtractor.cpp + ActsExamples/TruthTracking/TrackParameterSmearing.cpp ActsExamples/TruthTracking/TrackParameterSelector.cpp ActsExamples/TruthTracking/TrackModifier.cpp ActsExamples/TruthTracking/TrackTruthMatcher.cpp diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 0d1131d45cf..a365572cabc 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -13,9 +13,20 @@ "Default TruthSmeared TruthEstimated Orthogonal HoughTransform Gbts Hashing", ) -ParticleSmearingSigmas = namedtuple( - "ParticleSmearingSigmas", - ["d0", "d0PtA", "d0PtB", "z0", "z0PtA", "z0PtB", "t0", "phi", "theta", "ptRel"], +TrackSmearingSigmas = namedtuple( + "TrackSmearingSigmas", + [ + "loc0", + "loc0PtA", + "loc0PtB", + "loc1", + "loc1PtA", + "loc1PtB", + "time", + "phi", + "theta", + "ptRel", + ], defaults=[None] * 10, ) @@ -256,7 +267,7 @@ class VertexFinder(Enum): @acts.examples.NamedTypeArgs( seedingAlgorithm=SeedingAlgorithm, - particleSmearingSigmas=ParticleSmearingSigmas, + trackSmearingSigmas=TrackSmearingSigmas, seedFinderConfigArg=SeedFinderConfigArg, seedFinderOptionsArg=SeedFinderOptionsArg, seedFilterConfigArg=SeedFilterConfigArg, @@ -275,7 +286,7 @@ def addSeeding( layerMappingConfigFile: Optional[Union[Path, str]] = None, connector_inputConfigFile: Optional[Union[Path, str]] = None, seedingAlgorithm: SeedingAlgorithm = SeedingAlgorithm.Default, - particleSmearingSigmas: ParticleSmearingSigmas = ParticleSmearingSigmas(), + trackSmearingSigmas: TrackSmearingSigmas = TrackSmearingSigmas(), initialSigmas: Optional[list] = None, initialSigmaPtRel: Optional[float] = None, initialVarInflation: Optional[list] = None, @@ -313,15 +324,15 @@ def addSeeding( Json file for space point geometry selection. Not required for SeedingAlgorithm.TruthSmeared. seedingAlgorithm : SeedingAlgorithm, Default seeding algorithm to use: one of Default (no truth information used), TruthSmeared, TruthEstimated - particleSmearingSigmas : ParticleSmearingSigmas(d0, d0PtA, d0PtB, z0, z0PtA, z0PtB, t0, phi, theta, ptRel) - ParticleSmearing configuration. - Defaults specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp + trackSmearingSigmas : TrackSmearingSigmas(loc0, loc0PtA, loc0PtB, loc1, loc1PtA, loc1PtB, time, phi, theta, ptRel) + TrackSmearing configuration. + Defaults specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.hpp initialSigmas : list Sets the initial covariance matrix diagonal. This is ignored in case of TruthSmearing. Defaults specified in Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp initialVarInflation : list List of 6 scale factors to inflate the initial covariance matrix - Defaults (all 1) specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp + Defaults (all 1) specified in Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TrackParameterSmearing.hpp seedFinderConfigArg : SeedFinderConfigArg(maxSeedsPerSpM, cotThetaMax, sigmaScattering, radLengthPerSeed, minPt, impactMax, deltaPhiMax, interactionPointCut, deltaZMax, maxPtScattering, zBinEdges, zBinsCustomLooping, rRangeMiddleSP, useVariableMiddleSPRange, binSizeR, seedConfirmation, centralSeedConfirmationRange, forwardSeedConfirmationRange, deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z) SeedFinderConfig settings. deltaR, deltaRBottomSP, deltaRTopSP, deltaRMiddleSPRange, collisionRegion, r, z. Defaults specified in Core/include/Acts/Seeding/SeedFinderConfig.hpp @@ -366,7 +377,7 @@ def addSeeding( s=s, rnd=rnd, selectedParticles=selectedParticles, - particleSmearingSigmas=particleSmearingSigmas, + trackSmearingSigmas=trackSmearingSigmas, initialSigmas=initialSigmas, initialSigmaPtRel=initialSigmaPtRel, initialVarInflation=initialVarInflation, @@ -520,7 +531,7 @@ def addTruthSmearedSeeding( s: acts.examples.Sequencer, rnd: Optional[acts.examples.RandomNumbers], selectedParticles: str, - particleSmearingSigmas: ParticleSmearingSigmas, + trackSmearingSigmas: TrackSmearingSigmas, initialSigmas: Optional[List[float]], initialSigmaPtRel: Optional[float], initialVarInflation: Optional[List[float]], @@ -532,31 +543,39 @@ def addTruthSmearedSeeding( """ rnd = rnd or acts.examples.RandomNumbers(seed=42) - # Run particle smearing - ptclSmear = acts.examples.ParticleSmearing( + + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( level=logLevel, inputParticles=selectedParticles, + outputTrackParameters="trueparameters", + ) + s.addAlgorithm(trkParamExtractor) + + # Smearing track parameters + trkSmear = acts.examples.TrackParameterSmearing( + level=logLevel, + inputTrackParameters=trkParamExtractor.config.outputTrackParameters, outputTrackParameters="estimatedparameters", randomNumbers=rnd, # gaussian sigmas to smear particle parameters **acts.examples.defaultKWArgs( - sigmaD0=particleSmearingSigmas.d0, - sigmaD0PtA=particleSmearingSigmas.d0PtA, - sigmaD0PtB=particleSmearingSigmas.d0PtB, - sigmaZ0=particleSmearingSigmas.z0, - sigmaZ0PtA=particleSmearingSigmas.z0PtA, - sigmaZ0PtB=particleSmearingSigmas.z0PtB, - sigmaT0=particleSmearingSigmas.t0, - sigmaPhi=particleSmearingSigmas.phi, - sigmaTheta=particleSmearingSigmas.theta, - sigmaPtRel=particleSmearingSigmas.ptRel, + sigmaLoc0=trackSmearingSigmas.loc0, + sigmaLoc0PtA=trackSmearingSigmas.loc0PtA, + sigmaLoc0PtB=trackSmearingSigmas.loc0PtB, + sigmaLoc1=trackSmearingSigmas.loc1, + sigmaLoc1PtA=trackSmearingSigmas.loc1PtA, + sigmaLoc1PtB=trackSmearingSigmas.loc1PtB, + sigmaTime=trackSmearingSigmas.time, + sigmaPhi=trackSmearingSigmas.phi, + sigmaTheta=trackSmearingSigmas.theta, + sigmaPtRel=trackSmearingSigmas.ptRel, initialSigmas=initialSigmas, initialSigmaPtRel=initialSigmaPtRel, initialVarInflation=initialVarInflation, particleHypothesis=particleHypothesis, ), ) - s.addAlgorithm(ptclSmear) + s.addAlgorithm(trkSmear) truthTrkFndAlg = acts.examples.TruthTrackFinder( level=logLevel, diff --git a/Examples/Python/src/Generators.cpp b/Examples/Python/src/Generators.cpp index 454227ddbf6..bc24b8e43f7 100644 --- a/Examples/Python/src/Generators.cpp +++ b/Examples/Python/src/Generators.cpp @@ -16,8 +16,6 @@ #include "ActsExamples/Generators/ParametricParticleGenerator.hpp" #include "ActsExamples/Generators/VertexGenerators.hpp" -#include -#include #include #include #include diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index 1502ce23d40..319ee24e2ee 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -10,9 +10,10 @@ #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/TruthTracking/HitSelector.hpp" #include "ActsExamples/TruthTracking/ParticleSelector.hpp" -#include "ActsExamples/TruthTracking/ParticleSmearing.hpp" +#include "ActsExamples/TruthTracking/ParticleTrackParamExtractor.hpp" #include "ActsExamples/TruthTracking/TrackModifier.hpp" #include "ActsExamples/TruthTracking/TrackParameterSelector.hpp" +#include "ActsExamples/TruthTracking/TrackParameterSmearing.hpp" #include "ActsExamples/TruthTracking/TrackTruthMatcher.hpp" #include "ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp" #include "ActsExamples/TruthTracking/TruthTrackFinder.hpp" @@ -41,12 +42,16 @@ void addTruthTracking(Context& ctx) { ActsExamples::TruthTrackFinder, mex, "TruthTrackFinder", inputParticles, inputMeasurementParticlesMap, outputProtoTracks); + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::ParticleTrackParamExtractor, mex, + "ParticleTrackParamExtractor", inputParticles, + outputTrackParameters); + ACTS_PYTHON_DECLARE_ALGORITHM( - ActsExamples::ParticleSmearing, mex, "ParticleSmearing", inputParticles, - outputTrackParameters, sigmaD0, sigmaD0PtA, sigmaD0PtB, sigmaZ0, - sigmaZ0PtA, sigmaZ0PtB, sigmaT0, sigmaPhi, sigmaTheta, sigmaPtRel, - initialSigmas, initialSigmaPtRel, initialVarInflation, particleHypothesis, - randomNumbers); + ActsExamples::TrackParameterSmearing, mex, "TrackParameterSmearing", + inputTrackParameters, outputTrackParameters, sigmaLoc0, sigmaLoc0PtA, + sigmaLoc0PtB, sigmaLoc1, sigmaLoc1PtA, sigmaLoc1PtB, sigmaTime, sigmaPhi, + sigmaTheta, sigmaPtRel, initialSigmas, initialSigmaPtRel, + initialVarInflation, particleHypothesis, randomNumbers); { using Alg = ActsExamples::ParticleSelector; diff --git a/Examples/Python/tests/conftest.py b/Examples/Python/tests/conftest.py index 71c0899fe3f..4bd8a430730 100644 --- a/Examples/Python/tests/conftest.py +++ b/Examples/Python/tests/conftest.py @@ -209,19 +209,12 @@ def _basic_prop_seq_factory(geo, s=None): rnd=rng, ) - # Run particle smearing - trackParametersGenerator = acts.examples.ParticleSmearing( - level=acts.logging.INFO, + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( + level=acts.logging.WARNING, inputParticles="particles_input", - outputTrackParameters="start_parameters", - randomNumbers=rng, - sigmaD0=0.0, - sigmaZ0=0.0, - sigmaPhi=0.0, - sigmaTheta=0.0, - sigmaPtRel=0.0, + outputTrackParameters="params_particles_input", ) - s.addAlgorithm(trackParametersGenerator) + s.addAlgorithm(trkParamExtractor) nav = acts.Navigator(trackingGeometry=geo) stepper = acts.StraightLineStepper() @@ -232,11 +225,11 @@ def _basic_prop_seq_factory(geo, s=None): level=acts.logging.WARNING, propagatorImpl=prop, sterileLogger=False, - inputTrackParameters="start_parameters", + inputTrackParameters="params_particles_input", outputSummaryCollection="propagation_summary", ) - s.addAlgorithm(alg) + return s, alg return _basic_prop_seq_factory diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index ab5e04706bb..6eeaaebff44 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -16,7 +16,7 @@ test_itk_seeding__estimatedparams.root: fc042037f12a434f2236df7d225b8ca24209b691 test_itk_seeding__performance_seeding.root: 78ebda54cd0f026ba4b7f316724ffd946de56a932735914baf1b7bba9505c29d test_itk_seeding__particles.root: 907ff693262c0db14b12c74b16586cb20d79caf5f03f93b178943e41ed35a1b6 test_itk_seeding__particles_simulation.root: ef0246069aa697019f28a8b270a68de95312cae5f2f2c74848566c3ce4f70363 -test_propagation__propagation_summary.root: 400043dfeed9eda16512ec5fd4a1389f94b98565083f06f439c7b7b59e46c544 +test_propagation__propagation_summary.root: de0c105ab0de0096241855fe3db46d7c5d054d897625ede4df276472a7e43c57 test_material_recording__geant4_material_tracks.root: c022b9362249b29f57a07926b20644e3ab4ab8ebcf03f773fbf46c446fc1a0a1 test_truth_tracking_gsf[generic]__trackstates_gsf.root: 4df2c69d5dd7d5446a547651e4e962daf17924f5c8617165a93a3223c8ba18fd test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 8c01d139cb865afa1959c62dbca76f3a1fb8b684c57ea4c2968baa6ffedadb6f @@ -24,9 +24,9 @@ test_truth_tracking_gsf[odd]__trackstates_gsf.root: c7397e53ea093f2432943ae263fc test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 4562341f12a61ea0d5e25872b6bf466b79a73781dc95fc18ef9c6515f0a47916 test_particle_gun__particles.root: 669d0304eb8bcf244aa627809a117944e5e3b994fdfcfb8710f2b9a8f9a62d3b test_material_mapping__material-map_tracks.root: 938b1a855369e9304401cb10d2751df3fd7acf32a2573f2057eb1691cd94edf3 -test_material_mapping__propagation-material.root: 5eacd0cb804d381171c8fb65d7b2d36e1d57db9f4cb2b58c0f24703479218cbf +test_material_mapping__propagation-material.root: e2b4eade0d8124c03c89e01bf6ff5029dd12e3c9efc0a19c22a12c5cd2800e77 test_volume_material_mapping__material-map-volume_tracks.root: 98e212d32ca054fa3d01af4167c1f49755a139d43b82c57908197f5985e0a4ff -test_volume_material_mapping__propagation-volume-material.root: 3e9d38cc541a1956b2f33be320d457559bb230311130a8531bf09371c272f913 +test_volume_material_mapping__propagation-volume-material.root: 42bb2fd9c50d44210c914e487903bbdef65a0491f0a85a7159c0e078ac983a80 test_digitization_example[smeared]__measurements.root: 2b583b886b76f94786c6c2832aa84d176059cbbc35882fb34b4fd1f22dd4c006 test_digitization_example[geometric]__measurements.root: 85efa861d14207fd7d1798dab093edcd0a2685bc189b4fc9a96b07d1001013a0 test_digitization_example_input[smeared]__particles.root: 669d0304eb8bcf244aa627809a117944e5e3b994fdfcfb8710f2b9a8f9a62d3b @@ -56,9 +56,9 @@ test_vertex_fitting_reading[AMVF-False-100]__performance_vertexing.root: 009e4b1 test_vertex_fitting_reading[AMVF-True-100]__performance_vertexing.root: 2d0dc1e02bfd1f7eaae26ef8ac657ce0291f70c7e4efddd35d171d31988a631e test_bfield_writing__solenoid.root: 7be51f0ed9cb99f59ae0271ba79cdb84635e6ee3d2109ea8a4b521875029c21d test_bfield_writing__solenoid2.root: 2db149336c9cd749dc50025076b49f9bc0586d53792b87a0fdd7f21a649a01a5 -test_root_prop_step_writer[configPosConstructor]__prop_steps.root: a783d525eebb4737e6e6bcf20d63d6f35520f4bcf23301ae8c4657309cbccdda -test_root_prop_step_writer[configKwConstructor]__prop_steps.root: a783d525eebb4737e6e6bcf20d63d6f35520f4bcf23301ae8c4657309cbccdda -test_root_prop_step_writer[kwargsConstructor]__prop_steps.root: a783d525eebb4737e6e6bcf20d63d6f35520f4bcf23301ae8c4657309cbccdda +test_root_prop_step_writer[configPosConstructor]__prop_steps.root: 200ece7cde60eb0dd8835cf4830835e9058d897ba5759099d96914b06a1df092 +test_root_prop_step_writer[configKwConstructor]__prop_steps.root: 200ece7cde60eb0dd8835cf4830835e9058d897ba5759099d96914b06a1df092 +test_root_prop_step_writer[kwargsConstructor]__prop_steps.root: 200ece7cde60eb0dd8835cf4830835e9058d897ba5759099d96914b06a1df092 test_root_particle_writer[configPosConstructor]__particles.root: e5d723e138b4e121c6e74a6dba072072f622995e117a40b8e63755ac784baad6 test_root_particle_writer[configKwConstructor]__particles.root: e5d723e138b4e121c6e74a6dba072072f622995e117a40b8e63755ac784baad6 test_root_particle_writer[kwargsConstructor]__particles.root: e5d723e138b4e121c6e74a6dba072072f622995e117a40b8e63755ac784baad6 diff --git a/Examples/Python/tests/test_algorithms.py b/Examples/Python/tests/test_algorithms.py index 4cca7199291..6311dc6460a 100644 --- a/Examples/Python/tests/test_algorithms.py +++ b/Examples/Python/tests/test_algorithms.py @@ -14,7 +14,7 @@ TruthTrackFinder, ParticleSelector, TruthVertexFinder, - ParticleSmearing, + TrackParameterSmearing, TrackSelectorAlgorithm, TrackFittingAlgorithm, SurfaceSortingAlgorithm, @@ -44,7 +44,7 @@ TruthTrackFinder, ParticleSelector, TruthVertexFinder, - ParticleSmearing, + TrackParameterSmearing, TrackSelectorAlgorithm, TrackFittingAlgorithm, SurfaceSortingAlgorithm, diff --git a/Examples/Python/tests/test_propagation.py b/Examples/Python/tests/test_propagation.py index 6583098f427..adfc55ba8c2 100644 --- a/Examples/Python/tests/test_propagation.py +++ b/Examples/Python/tests/test_propagation.py @@ -56,34 +56,27 @@ def test_steppers(conf_const, trk_geo): rnd=rnd, ) - # Run particle smearing - trackParametersGenerator = acts.examples.ParticleSmearing( - level=acts.logging.INFO, - inputParticles="particles_input", - outputTrackParameters="start_parameters", - randomNumbers=rnd, - sigmaD0=0.0, - sigmaZ0=0.0, - sigmaPhi=0.0, - sigmaTheta=0.0, - sigmaPtRel=0.0, - ) - seq.addAlgorithm(trackParametersGenerator) - prop = acts.examples.ConcretePropagator( acts.Propagator(stepper=s, navigator=nav) ) + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( + level=acts.logging.WARNING, + inputParticles="particles_input", + outputTrackParameters="params_particles_input", + ) + seq.addAlgorithm(trkParamExtractor) + alg = conf_const( acts.examples.PropagationAlgorithm, level=acts.logging.WARNING, propagatorImpl=prop, - inputTrackParameters="start_parameters", + inputTrackParameters="params_particles_input", outputSummaryCollection="propagation_summary", sterileLogger=False, ) - seq.addAlgorithm(alg) + chkAlg = AssertCollectionExistsAlg( "propagation_summary", "chk_alg", level=acts.logging.WARNING ) diff --git a/Examples/Scripts/Optimization/ckf.py b/Examples/Scripts/Optimization/ckf.py index ba18291735c..e0cd5fd9870 100755 --- a/Examples/Scripts/Optimization/ckf.py +++ b/Examples/Scripts/Optimization/ckf.py @@ -123,7 +123,7 @@ def runCKFTracks( from acts.examples.reconstruction import ( addSeeding, - ParticleSmearingSigmas, + TrackSmearingSigmas, SeedFinderConfigArg, SeedFinderOptionsArg, SeedingAlgorithm, @@ -188,15 +188,15 @@ def runCKFTracks( s, trackingGeometry, field, - ParticleSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared + TrackSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared # zero eveything so the CKF has a chance to find the measurements - d0=0, - d0PtA=0, - d0PtB=0, - z0=0, - z0PtA=0, - z0PtB=0, - t0=0, + loc0=0, + loc0PtA=0, + loc0PtB=0, + loc1=0, + loc1PtA=0, + loc1PtB=0, + time=0, phi=0, theta=0, ptRel=0, diff --git a/Examples/Scripts/Python/ckf_tracks.py b/Examples/Scripts/Python/ckf_tracks.py index 97f56d17ff0..75ec8a7c316 100755 --- a/Examples/Scripts/Python/ckf_tracks.py +++ b/Examples/Scripts/Python/ckf_tracks.py @@ -34,7 +34,7 @@ def runCKFTracks( from acts.examples.reconstruction import ( addSeeding, - ParticleSmearingSigmas, + TrackSmearingSigmas, SeedFinderConfigArg, SeedFinderOptionsArg, SeedingAlgorithm, @@ -99,15 +99,15 @@ def runCKFTracks( s, trackingGeometry, field, - ParticleSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared + TrackSmearingSigmas( # only used by SeedingAlgorithm.TruthSmeared # zero eveything so the CKF has a chance to find the measurements - d0=0, - d0PtA=0, - d0PtB=0, - z0=0, - z0PtA=0, - z0PtB=0, - t0=0, + loc0=0, + loc0PtA=0, + loc0PtB=0, + loc1=0, + loc1PtA=0, + loc1PtB=0, + time=0, phi=0, theta=0, ptRel=0, diff --git a/Examples/Scripts/Python/full_chain_test.py b/Examples/Scripts/Python/full_chain_test.py index f0010bd498e..595dbd43a4e 100755 --- a/Examples/Scripts/Python/full_chain_test.py +++ b/Examples/Scripts/Python/full_chain_test.py @@ -521,7 +521,7 @@ def full_chain(args): from acts.examples.reconstruction import ( addSeeding, - ParticleSmearingSigmas, + TrackSmearingSigmas, addCKFTracks, CkfConfig, SeedingAlgorithm, @@ -547,7 +547,7 @@ def full_chain(args): seedingAlgorithm=args.seeding_algorithm, **( dict( - particleSmearingSigmas=ParticleSmearingSigmas(ptRel=0.01), + trackSmearingSigmas=TrackSmearingSigmas(ptRel=0.01), rnd=rnd, ) if args.seeding_algorithm == SeedingAlgorithm.TruthSmeared diff --git a/Examples/Scripts/Python/material_validation.py b/Examples/Scripts/Python/material_validation.py index ebc26041cd2..94c0f21b94a 100755 --- a/Examples/Scripts/Python/material_validation.py +++ b/Examples/Scripts/Python/material_validation.py @@ -40,30 +40,22 @@ def runMaterialValidation( rnd=rnd, ) - # Run particle smearing - trackParametersGenerator = acts.examples.ParticleSmearing( + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( level=acts.logging.INFO, inputParticles="particles_input", - outputTrackParameters="start_parameters", - randomNumbers=rnd, - sigmaD0=0.0, - sigmaZ0=0.0, - sigmaPhi=0.0, - sigmaTheta=0.0, - sigmaPtRel=0.0, + outputTrackParameters="params_particles_input", ) - s.addAlgorithm(trackParametersGenerator) + s.addAlgorithm(trkParamExtractor) alg = acts.examples.PropagationAlgorithm( propagatorImpl=prop, level=acts.logging.INFO, sterileLogger=True, recordMaterialInteractions=True, - inputTrackParameters="start_parameters", + inputTrackParameters="params_particles_input", outputSummaryCollection="propagation_summary", outputMaterialCollection="material_tracks", ) - s.addAlgorithm(alg) s.addWriter( diff --git a/Examples/Scripts/Python/material_validation_itk.py b/Examples/Scripts/Python/material_validation_itk.py index 759f75a5ea2..4d5f095ef10 100755 --- a/Examples/Scripts/Python/material_validation_itk.py +++ b/Examples/Scripts/Python/material_validation_itk.py @@ -10,6 +10,8 @@ def runMaterialValidation( + nevents, + ntracks, trackingGeometry, decorators, field, @@ -18,7 +20,8 @@ def runMaterialValidation( dumpPropagationSteps=False, s=None, ): - s = s or Sequencer(events=1000, numThreads=-1) + # Create a sequencer + s = s or Sequencer(events=nevents, numThreads=-1) rnd = acts.examples.RandomNumbers(seed=42) @@ -44,31 +47,22 @@ def runMaterialValidation( rnd=rnd, ) - # Run particle smearing - trackParametersGenerator = acts.examples.ParticleSmearing( + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( level=acts.logging.INFO, inputParticles="particles_input", - outputTrackParameters="start_parameters", - randomNumbers=rnd, - sigmaD0=0.0, - sigmaZ0=0.0, - sigmaPhi=0.0, - sigmaTheta=0.0, - sigmaPRel=0.0, - addCovariances=False, + outputTrackParameters="params_particles_input", ) - s.addAlgorithm(trackParametersGenerator) + s.addAlgorithm(trkParamExtractor) alg = acts.examples.PropagationAlgorithm( propagatorImpl=prop, level=acts.logging.INFO, sterileLogger=False, recordMaterialInteractions=True, - inputTrackParameters="start_parameters", + inputTrackParameters="params_particles_input", outputPropagationSteps="propagation_steps", outputMaterialTracks="material-tracks", ) - s.addAlgorithm(alg) s.addWriter( diff --git a/Examples/Scripts/Python/propagation.py b/Examples/Scripts/Python/propagation.py index ea54919a89f..59be3385ca0 100755 --- a/Examples/Scripts/Python/propagation.py +++ b/Examples/Scripts/Python/propagation.py @@ -31,19 +31,12 @@ def runPropagation(trackingGeometry, field, outputDir, s=None, decorators=[]): rnd=rnd, ) - # Run particle smearing - trackParametersGenerator = acts.examples.ParticleSmearing( - level=acts.logging.INFO, + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( + level=acts.logging.WARNING, inputParticles="particles_input", - outputTrackParameters="start_parameters", - randomNumbers=rnd, - sigmaD0=0.0, - sigmaZ0=0.0, - sigmaPhi=0.0, - sigmaTheta=0.0, - sigmaPtRel=0.0, + outputTrackParameters="params_particles_input", ) - s.addAlgorithm(trackParametersGenerator) + s.addAlgorithm(trkParamExtractor) nav = acts.Navigator(trackingGeometry=trackingGeometry) @@ -57,7 +50,7 @@ def runPropagation(trackingGeometry, field, outputDir, s=None, decorators=[]): propagatorImpl=propagator, level=acts.logging.INFO, sterileLogger=True, - inputTrackParameters="start_parameters", + inputTrackParameters="params_particles_input", outputSummaryCollection="propagation_summary", ) s.addAlgorithm(propagationAlgorithm) diff --git a/Examples/Scripts/Python/vertex_fitting.py b/Examples/Scripts/Python/vertex_fitting.py index e8797a5ed12..c6088f618a6 100755 --- a/Examples/Scripts/Python/vertex_fitting.py +++ b/Examples/Scripts/Python/vertex_fitting.py @@ -6,7 +6,7 @@ from acts.examples import ( Sequencer, ParticleSelector, - ParticleSmearing, + TrackParameterSmearing, TrackParameterSelector, ) from acts.examples.simulation import addPythia8 @@ -65,9 +65,16 @@ def runVertexFitting( if inputTrackSummary is None or inputParticlePath is None: logger.info("Using smeared particles") - ptclSmearing = ParticleSmearing( - level=acts.logging.INFO, + trkParamExtractor = acts.examples.ParticleTrackParamExtractor( + level=acts.logging.WARNING, inputParticles=selectedParticles, + outputTrackParameters="params_particles_input", + ) + s.addAlgorithm(trkParamExtractor) + + ptclSmearing = TrackParameterSmearing( + level=acts.logging.INFO, + inputTrackParameters="params_particles_input", outputTrackParameters=trackParameters, randomNumbers=rnd, ) From eb6a5f11df765261ac9f1b75e0eeb8f3e855eb59 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 29 Nov 2024 18:10:50 +0100 Subject: [PATCH 3/8] feat: compatibility for EDM4hep < `0.99` & >= `0.99` (#3921) This updates our code so it compiles with both EDM4hep versions below `0.99` and above/equal to `0.99`. Breaking changes are: - `TrackerHit` becomes `TrackerHit3D`. - `SimTrackerHit::getMCParticle` becomes `SimTrackerHit::getParticle`, with the original name being deprecated. I'm aliasing the new name `TrackerHit3D` and `TrackerHit3DCollection` to the old types, so that when we drop support for < `0.99` versions, they're already correct. ## Summary by CodeRabbit - **New Features** - Enhanced compatibility with different versions of the EDM4hep library for particle handling. - Introduced utility functions for retrieving and setting particle information in `SimTrackerHit` objects. - Added type aliases to streamline the handling of tracker hit collections. - **Bug Fixes** - Improved error handling in the `EDM4hepReader` for missing particles. - **Refactor** - Updated function signatures and internal logic to utilize new utility functions for better maintainability. - Streamlined the `EDM4hepMeasurementWriter` and `EDM4hepMeasurementReader` to accommodate new data structures. - Adjusted header inclusions and conditional compilation for better flexibility. - **Documentation** - Added comments to indicate changes in parameter handling for clarity. --- .../Io/EDM4hep/EDM4hepMeasurementWriter.hpp | 3 - .../ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp | 21 +++--- .../EDM4hep/src/EDM4hepMeasurementReader.cpp | 5 +- .../EDM4hep/src/EDM4hepMeasurementWriter.cpp | 5 +- Examples/Io/EDM4hep/src/EDM4hepReader.cpp | 6 +- Examples/Io/EDM4hep/src/EDM4hepUtil.cpp | 73 ++++++------------- .../Acts/Plugins/EDM4hep/EDM4hepUtil.hpp | 12 ++- .../EDM4hep/TrackerHitCompatibility.hpp | 22 ++++++ Plugins/EDM4hep/src/EDM4hepUtil.cpp | 34 ++++++++- 9 files changed, 105 insertions(+), 76 deletions(-) create mode 100644 Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp index d559ee9a5e4..04cf8805df4 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp @@ -16,9 +16,6 @@ #include -#include -#include - namespace ActsExamples { /// Write out a measurement cluster collection to EDM4hep. diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp index bdce08bc876..4973763a087 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp @@ -8,6 +8,7 @@ #pragma once #include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/SimParticle.hpp" @@ -16,15 +17,13 @@ #include -#include "edm4hep/MCParticle.h" -#include "edm4hep/MutableMCParticle.h" -#include "edm4hep/MutableSimTrackerHit.h" -#include "edm4hep/MutableTrack.h" -#include "edm4hep/MutableTrackerHitPlane.h" -#include "edm4hep/SimTrackerHit.h" -#include "edm4hep/TrackerHit.h" -#include "edm4hep/TrackerHitCollection.h" -#include "edm4hep/TrackerHitPlane.h" +#include +#include +#include +#include +#include +#include +#include namespace ActsExamples::EDM4hepUtil { @@ -91,7 +90,7 @@ void writeSimHit(const ActsFatras::Hit& from, edm4hep::MutableSimTrackerHit to, /// - local 2D coordinates and time are read from position VariableBoundMeasurementProxy readMeasurement( MeasurementContainer& container, const edm4hep::TrackerHitPlane& from, - const edm4hep::TrackerHitCollection* fromClusters, Cluster* toCluster, + const edm4hep::TrackerHit3DCollection* fromClusters, Cluster* toCluster, const MapGeometryIdFrom& geometryMapper); /// Writes a measurement cluster to EDM4hep. @@ -107,7 +106,7 @@ VariableBoundMeasurementProxy readMeasurement( void writeMeasurement(const ConstVariableBoundMeasurementProxy& from, edm4hep::MutableTrackerHitPlane to, const Cluster* fromCluster, - edm4hep::TrackerHitCollection& toClusters, + edm4hep::TrackerHit3DCollection& toClusters, const MapGeometryIdTo& geometryMapper); /// Writes a trajectory to EDM4hep. diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp index 237ff19e3b6..5299da121ca 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp @@ -9,6 +9,7 @@ #include "ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp" #include "Acts/Definitions/Units.hpp" +#include "Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp" #include "Acts/Plugins/Podio/PodioUtil.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" @@ -18,8 +19,6 @@ #include #include -#include -#include #include #include @@ -60,7 +59,7 @@ ProcessCode EDM4hepMeasurementReader::read(const AlgorithmContext& ctx) { const auto& trackerHitPlaneCollection = frame.get("ActsTrackerHitsPlane"); const auto& trackerHitRawCollection = - frame.get("ActsTrackerHitsRaw"); + frame.get("ActsTrackerHitsRaw"); for (const auto& trackerHitPlane : trackerHitPlaneCollection) { Cluster cluster; diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp index 0f9e9422d28..23e291cc2dd 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp @@ -9,6 +9,7 @@ #include "ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp" #include "Acts/Definitions/Units.hpp" +#include "Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Framework/WhiteBoard.hpp" @@ -16,6 +17,8 @@ #include +#include +#include #include namespace ActsExamples { @@ -44,7 +47,7 @@ ActsExamples::ProcessCode EDM4hepMeasurementWriter::writeT( podio::Frame frame; edm4hep::TrackerHitPlaneCollection hitsPlane; - edm4hep::TrackerHitCollection hits; + edm4hep::TrackerHit3DCollection hits; if (!m_cfg.inputClusters.empty()) { ACTS_VERBOSE("Fetch clusters for writing: " << m_cfg.inputClusters); diff --git a/Examples/Io/EDM4hep/src/EDM4hepReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepReader.cpp index 4561bc61db0..0b82a20181b 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepReader.cpp @@ -10,6 +10,7 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" +#include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp" #include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" #include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" @@ -305,8 +306,9 @@ ProcessCode EDM4hepReader::read(const AlgorithmContext& ctx) { auto simHit = EDM4hepUtil::readSimHit( hit, [&](const auto& inParticle) { - ACTS_VERBOSE("SimHit has source particle: " - << hit.getMCParticle().getObjectID().index); + ACTS_VERBOSE( + "SimHit has source particle: " + << Acts::EDM4hepUtil::getParticle(hit).getObjectID().index); auto it = edm4hepParticleMap.find(inParticle.getObjectID().index); if (it == edm4hepParticleMap.end()) { ACTS_ERROR( diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp index 38e26038779..adb433891d9 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp @@ -70,9 +70,10 @@ void EDM4hepUtil::writeParticle(const SimParticle& from, ActsFatras::Hit EDM4hepUtil::readSimHit( const edm4hep::SimTrackerHit& from, const MapParticleIdFrom& particleMapper, const MapGeometryIdFrom& geometryMapper) { - ActsFatras::Barcode particleId = particleMapper(from.getMCParticle()); + auto particle = Acts::EDM4hepUtil::getParticle(from); + ActsFatras::Barcode particleId = particleMapper(particle); - const auto mass = from.getMCParticle().getMass() * 1_GeV; + const auto mass = particle.getMass() * 1_GeV; const Acts::Vector3 momentum{ from.getMomentum().x * 1_GeV, from.getMomentum().y * 1_GeV, @@ -116,7 +117,7 @@ void EDM4hepUtil::writeSimHit(const ActsFatras::Hit& from, const auto delta4 = from.momentum4After() - momentum4Before; if (particleMapper) { - to.setMCParticle(particleMapper(from.particleId())); + Acts::EDM4hepUtil::setParticle(to, particleMapper(from.particleId())); } if (geometryMapper) { @@ -146,8 +147,8 @@ void EDM4hepUtil::writeSimHit(const ActsFatras::Hit& from, VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement( MeasurementContainer& container, const edm4hep::TrackerHitPlane& from, - const edm4hep::TrackerHitCollection* fromClusters, Cluster* toCluster, - const MapGeometryIdFrom& geometryMapper) { + const edm4hep::TrackerHit3DCollection* /*fromClusters*/, + Cluster* /*toCluster*/, const MapGeometryIdFrom& geometryMapper) { // no need for digitization as we only want to identify the sensor Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID()); @@ -171,32 +172,15 @@ VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement( auto to = createMeasurement(container, geometryId, dParameters); - if (fromClusters != nullptr) { - for (const auto objectId : from.getRawHits()) { - const auto& c = fromClusters->at(objectId.index); - - // TODO get EDM4hep fixed - // misusing some fields to store ACTS specific information - // don't ask ... - ActsFatras::Segmentizer::Bin2D bin{ - static_cast(c.getType()), - static_cast(c.getQuality())}; - ActsFatras::Segmentizer::Segment2D path2D{ - {Acts::Vector2::Zero(), Acts::Vector2::Zero()}}; - double activation = c.getTime(); - ActsFatras::Segmentizer::ChannelSegment cell{bin, path2D, activation}; - - toCluster->channels.push_back(cell); - } - } + // @TODO: Figure out if cell information is accessible return to; } void EDM4hepUtil::writeMeasurement( const ConstVariableBoundMeasurementProxy& from, - edm4hep::MutableTrackerHitPlane to, const Cluster* fromCluster, - edm4hep::TrackerHitCollection& toClusters, + edm4hep::MutableTrackerHitPlane to, const Cluster* /*fromCluster*/, + edm4hep::TrackerHit3DCollection& /*toClusters*/, const MapGeometryIdTo& geometryMapper) { Acts::GeometryIdentifier geoId = from.geometryId(); @@ -224,21 +208,7 @@ void EDM4hepUtil::writeMeasurement( 0, }); - if (fromCluster != nullptr) { - for (const auto& c : fromCluster->channels) { - auto toChannel = toClusters.create(); - to.addToRawHits(toChannel.getObjectID()); - - // TODO digitization channel - - // TODO get EDM4hep fixed - // misusing some fields to store ACTS specific information - // don't ask ... - toChannel.setType(c.bin[0]); - toChannel.setQuality(c.bin[1]); - toChannel.setTime(c.activation); - } - } + // @TODO: Check if we can write cell info } void EDM4hepUtil::writeTrajectory( @@ -296,17 +266,18 @@ void EDM4hepUtil::writeTrajectory( trackState.referencePoint.z = center.z(); if (converted.covariance) { - const auto& c = converted.covariance.value(); - - trackState.covMatrix = { - static_cast(c(0, 0)), static_cast(c(1, 0)), - static_cast(c(1, 1)), static_cast(c(2, 0)), - static_cast(c(2, 1)), static_cast(c(2, 2)), - static_cast(c(3, 0)), static_cast(c(3, 1)), - static_cast(c(3, 2)), static_cast(c(3, 3)), - static_cast(c(4, 0)), static_cast(c(4, 1)), - static_cast(c(4, 2)), static_cast(c(4, 3)), - static_cast(c(4, 4))}; + auto c = [&](std::size_t row, std::size_t col) { + return static_cast(converted.covariance.value()(row, col)); + }; + + // clang-format off + trackState.covMatrix = {c(0, 0), + c(1, 0), c(1, 1), + c(2, 0), c(2, 1), c(2, 2), + c(3, 0), c(3, 1), c(3, 2), c(3, 3), + c(4, 0), c(4, 1), c(4, 2), c(4, 3), c(4, 4), + c(5, 0), c(5, 1), c(5, 2), c(5, 3), c(5, 4), c(5, 5)}; + // clang-format on } to.addToTrackStates(trackState); diff --git a/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp index c5fda03bc82..e4836862769 100644 --- a/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp +++ b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp @@ -26,11 +26,13 @@ #include #include +#include +#include +#include +#include #include #include -#include "edm4hep/MutableTrack.h" - namespace Acts::EDM4hepUtil { static constexpr std::int32_t EDM4HEP_ACTS_POSITION_TYPE = 42; @@ -61,6 +63,12 @@ BoundTrackParameters convertTrackParametersFromEdm4hep( } // namespace detail +// Compatibility with EDM4hep < 0.99 and >= 0.99 +edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit); + +void setParticle(edm4hep::MutableSimTrackerHit& hit, + const edm4hep::MCParticle& particle); + template void writeTrack(const Acts::GeometryContext& gctx, track_proxy_t track, edm4hep::MutableTrack to, double Bz, diff --git a/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp new file mode 100644 index 00000000000..11115e7cbe5 --- /dev/null +++ b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/TrackerHitCompatibility.hpp @@ -0,0 +1,22 @@ +// 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 + +// Compatibility with EDM4hep < 0.99 and >= 0.99 +#if __has_include() +#include "edm4hep/TrackerHit3D.h" +#include "edm4hep/TrackerHit3DCollection.h" +#else +#include "edm4hep/TrackerHit.h" +#include "edm4hep/TrackerHitCollection.h" +namespace edm4hep { +using TrackerHit3DCollection = edm4hep::TrackerHitCollection; +using TrackerHit3D = edm4hep::TrackerHit; +} // namespace edm4hep +#endif diff --git a/Plugins/EDM4hep/src/EDM4hepUtil.cpp b/Plugins/EDM4hep/src/EDM4hepUtil.cpp index 515063da679..60d466d0427 100644 --- a/Plugins/EDM4hep/src/EDM4hepUtil.cpp +++ b/Plugins/EDM4hep/src/EDM4hepUtil.cpp @@ -18,9 +18,14 @@ #include -#include "edm4hep/TrackState.h" +#include +#include +#include +#include +#include -namespace Acts::EDM4hepUtil::detail { +namespace Acts::EDM4hepUtil { +namespace detail { ActsSquareMatrix<6> jacobianToEdm4hep(double theta, double qOverP, double Bz) { // Calculate jacobian from our internal parametrization (d0, z0, phi, theta, @@ -185,4 +190,27 @@ BoundTrackParameters convertTrackParametersFromEdm4hep( return {params.surface, targetPars, cov, params.particleHypothesis}; } -} // namespace Acts::EDM4hepUtil::detail +} // namespace detail + +#if EDM4HEP_VERSION_MAJOR >= 1 || \ + (EDM4HEP_VERSION_MAJOR == 0 && EDM4HEP_VERSION_MINOR == 99) +edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) { + return hit.getParticle(); +} + +void setParticle(edm4hep::MutableSimTrackerHit& hit, + const edm4hep::MCParticle& particle) { + hit.setParticle(particle); +} +#else +edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) { + return hit.getMCParticle(); +} + +void setParticle(edm4hep::MutableSimTrackerHit& hit, + const edm4hep::MCParticle& particle) { + hit.setMCParticle(particle); +} +#endif + +} // namespace Acts::EDM4hepUtil From 6a836db2a62deec29f0f1913f557a59088bb2978 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 29 Nov 2024 21:52:08 +0100 Subject: [PATCH 4/8] ci: Add LCG 106a job (#3918) This adds jobs for LCG 106a to the GitLab CI. I'm bumping ODD to the [`v4.0.3`](https://gitlab.cern.ch/acts/OpenDataDetector/-/tags/v4.0.3) where I cherry picked a tiny code fix that was breaking the build here. Blocked by: - #3921 ## Summary by CodeRabbit - **New Features** - Introduced a new CI/CD job `lcg_106a` for enhanced build configurations. - **Improvements** - Simplified the `lcg_105` job by directly setting `LCG_PLATFORM` for the `alma9` OS. - **Updates** - Updated the commit reference for the `OpenDataDetector` subproject to the latest version. --- .gitlab-ci.yml | 20 +++++++++++++++----- thirdparty/OpenDataDetector | 2 +- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b0a213e607d..c897c947010 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -382,11 +382,7 @@ linux_ubuntu_2204_clang: # Figure out LCG platform name based on version number and OS - > if [ "$OS" = "alma9" ]; then - if [ "$LCG_VERSION" -ge "104" ]; then - export LCG_PLATFORM="el9" - else - export LCG_PLATFORM="centos9" - fi + export LCG_PLATFORM="el9" else export LCG_PLATFORM="$OS" fi @@ -431,3 +427,17 @@ lcg_105: COMPILER: - gcc13 - clang16 + +lcg_106a: + extends: .lcg_base_job + + variables: + LCG_VERSION: "106a" + + parallel: + matrix: + - OS: [alma9] + COMPILER: + - gcc13 + - gcc14 + - clang16 diff --git a/thirdparty/OpenDataDetector b/thirdparty/OpenDataDetector index 0bf11ff33bd..4c5e7413a50 160000 --- a/thirdparty/OpenDataDetector +++ b/thirdparty/OpenDataDetector @@ -1 +1 @@ -Subproject commit 0bf11ff33bd10d3070e6ecbc4bea0a7c42d86fa0 +Subproject commit 4c5e7413a505a7743cbec9544901e8c55172e513 From 8895da1c915b56c9ea2498e32bf86dadb4f8b018 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 29 Nov 2024 23:19:49 +0100 Subject: [PATCH 5/8] fix: Remove `final` on template derived classes (#3923) This seems to cause issues when using `dynamic_cast` on AppleClang 16 and llvm 18 on macOS 15 (at least). Reproducer: https://gist.github.com/paulgessinger/2c1d2abdeb322072c507878ab5833728 ## Summary by CodeRabbit - **New Features** - Enhanced extensibility of the `GridPortalLinkT`, `Axis`, `Affine3Transformed`, `GlobalSubspace`, and `LocalSubspace` classes by allowing subclassing. - Updated constructors and methods in `GridPortalLinkT` to improve error handling and functionality. - **Bug Fixes** - Improved validation in the `resolveVolume` method to ensure positions are within valid bounds. - **Documentation** - Clarified class hierarchies and relationships due to changes in access specifiers. --- Core/include/Acts/Geometry/GridPortalLink.hpp | 2 +- Core/include/Acts/Utilities/Axis.hpp | 4 ++-- Core/include/Acts/Utilities/GridAccessHelpers.hpp | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Core/include/Acts/Geometry/GridPortalLink.hpp b/Core/include/Acts/Geometry/GridPortalLink.hpp index 411e6bb9962..8bf1d06ee08 100644 --- a/Core/include/Acts/Geometry/GridPortalLink.hpp +++ b/Core/include/Acts/Geometry/GridPortalLink.hpp @@ -396,7 +396,7 @@ class GridPortalLink : public PortalLinkBase { /// @tparam Axes The axis types of the grid template requires(sizeof...(Axes) <= 2) -class GridPortalLinkT final : public GridPortalLink { +class GridPortalLinkT : public GridPortalLink { public: /// The internal grid type using GridType = Grid; diff --git a/Core/include/Acts/Utilities/Axis.hpp b/Core/include/Acts/Utilities/Axis.hpp index 5b55d27cd77..758a8c2ec3c 100644 --- a/Core/include/Acts/Utilities/Axis.hpp +++ b/Core/include/Acts/Utilities/Axis.hpp @@ -97,7 +97,7 @@ class NeighborHoodIndices { /// This class provides some basic functionality for calculating bin indices /// for a given equidistant binning. template -class Axis final : public IAxis { +class Axis : public IAxis { public: static constexpr AxisType type = AxisType::Equidistant; @@ -417,7 +417,7 @@ class Axis final : public IAxis { /// This class provides some basic functionality for calculating bin indices /// for a given binning with variable bin sizes. template -class Axis final : public IAxis { +class Axis : public IAxis { public: static constexpr AxisType type = AxisType::Variable; diff --git a/Core/include/Acts/Utilities/GridAccessHelpers.hpp b/Core/include/Acts/Utilities/GridAccessHelpers.hpp index 085ea318f26..fbaf7f3acbd 100644 --- a/Core/include/Acts/Utilities/GridAccessHelpers.hpp +++ b/Core/include/Acts/Utilities/GridAccessHelpers.hpp @@ -110,7 +110,7 @@ class IBoundToGridLocal { }; template -class Affine3Transformed final : public IGlobalToGridLocal { +class Affine3Transformed : public IGlobalToGridLocal { public: using grid_local_t = typename global_to_grid_local_t::grid_local_t; @@ -142,7 +142,7 @@ class Affine3Transformed final : public IGlobalToGridLocal { /// position /// @tparam ...Args template -class GlobalSubspace final : public IGlobalToGridLocal { +class GlobalSubspace : public IGlobalToGridLocal { public: using grid_local_t = std::array; @@ -179,7 +179,7 @@ class GlobalSubspace final : public IGlobalToGridLocal { // The bound to grid local transformation, if only access of a subspace // is requested template -class LocalSubspace final : public IBoundToGridLocal { +class LocalSubspace : public IBoundToGridLocal { public: using grid_local_t = std::array; From bf997eb640537c9d91f0ee7f54e95216a680a98a Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Mon, 2 Dec 2024 13:36:44 +0100 Subject: [PATCH 6/8] refactor: modernise map utils (#3615) ## Summary by CodeRabbit - **New Features** - Introduced a new function to compute minimum and maximum values along with unique bin counts for input vectors. - **Refactor** - Enhanced the handling of magnetic field data by improving type usage and logic in multiple functions. - Streamlined calculations for minimum and maximum values using the new utility function. - Simplified the material mapping functions to improve maintainability and reduce complexity. - **Bug Fixes** - Adjusted error handling and control flow for grid value settings based on magnetic field data. --- .../Acts/MagneticField/BFieldMapUtils.hpp | 15 +- Core/include/Acts/Utilities/Helpers.hpp | 45 +++ Core/src/MagneticField/BFieldMapUtils.cpp | 269 +++++++----------- Core/src/Material/MaterialMapUtils.cpp | 98 ++----- 4 files changed, 175 insertions(+), 252 deletions(-) diff --git a/Core/include/Acts/MagneticField/BFieldMapUtils.hpp b/Core/include/Acts/MagneticField/BFieldMapUtils.hpp index 09ac0035d9b..4aa1934699b 100644 --- a/Core/include/Acts/MagneticField/BFieldMapUtils.hpp +++ b/Core/include/Acts/MagneticField/BFieldMapUtils.hpp @@ -75,7 +75,7 @@ fieldMapRZ(const std::function binsRZ, std::array nBinsRZ)>& localToGlobalBin, std::vector rPos, std::vector zPos, - std::vector bField, + const std::vector& bField, double lengthUnit = UnitConstants::mm, double BFieldUnit = UnitConstants::T, bool firstQuadrant = false); @@ -137,7 +137,7 @@ fieldMapXYZ( std::array nBinsXYZ)>& localToGlobalBin, std::vector xPos, std::vector yPos, - std::vector zPos, std::vector bField, + std::vector zPos, const std::vector& bField, double lengthUnit = UnitConstants::mm, double BFieldUnit = UnitConstants::T, bool firstOctant = false); @@ -145,17 +145,18 @@ fieldMapXYZ( /// creates a field mapper by sampling grid points from the analytical /// solenoid field. /// -/// @param rlim pair of r bounds -/// @param zlim pair of z bounds -/// @param nbins pair of bin counts +/// @param rLim pair of r bounds +/// @param zLim pair of z bounds +/// @param nBins pair of bin counts /// @param field the solenoid field instance /// /// @return A field map instance for use in interpolation. Acts::InterpolatedBFieldMap< Acts::Grid, Acts::Axis>> -solenoidFieldMap(std::pair rlim, std::pair zlim, - std::pair nbins, +solenoidFieldMap(const std::pair& rLim, + const std::pair& zLim, + const std::pair& nBins, const SolenoidBField& field); } // namespace Acts diff --git a/Core/include/Acts/Utilities/Helpers.hpp b/Core/include/Acts/Utilities/Helpers.hpp index 44bf23d9157..56d9b36481f 100644 --- a/Core/include/Acts/Utilities/Helpers.hpp +++ b/Core/include/Acts/Utilities/Helpers.hpp @@ -221,4 +221,49 @@ struct overloaded : Ts... { template overloaded(Ts...) -> overloaded; +namespace detail { + +/// Computes the minimum, maximum, and bin count for a given vector of values. +/// +/// This function processes a vector of doubles to compute: +/// - The minimum value (@c xMin) +/// - The maximum value (@c xMax), adjusted to include an additional bin +/// - The bin count (@c xBinCount) based on the number of unique values +/// +/// The computation is performed as follows: +/// 1. Sorts the input vector using @c std::ranges::sort to prepare for uniqueness. +/// 2. Determines the number of unique values using @c std::unique and calculates the bin count. +/// 3. Calculates the minimum and maximum using @c std::ranges::minmax. +/// 4. Adjusts the maximum to include an additional bin by adding the bin step +/// size. +/// +/// @param xPos A reference to a vector of doubles. +/// @return A tuple containing: +/// - The minimum value (double) +/// - The adjusted maximum value (double) +/// - The bin count (std::size_t) +/// +/// @note The vector xPos will be modified during the call. +inline auto getMinMaxAndBinCount(std::vector& xPos) { + // sort the values for unique() + std::ranges::sort(xPos); + + // get the number of bins over unique values + auto it = std::unique(xPos.begin(), xPos.end()); + const std::size_t xBinCount = std::distance(xPos.begin(), it); + + // get the minimum and maximum + auto [xMin, xMax] = std::ranges::minmax(xPos); + + // calculate maxima (add one last bin, because bin value always corresponds to + // left boundary) + const double stepX = (xMax - xMin) / static_cast(xBinCount - 1); + xMax += stepX; + + // Return all values as a tuple + return std::make_tuple(xMin, xMax, xBinCount); +} + +} // namespace detail + } // namespace Acts diff --git a/Core/src/MagneticField/BFieldMapUtils.cpp b/Core/src/MagneticField/BFieldMapUtils.cpp index b5b92191589..f5b1b410ffe 100644 --- a/Core/src/MagneticField/BFieldMapUtils.cpp +++ b/Core/src/MagneticField/BFieldMapUtils.cpp @@ -12,12 +12,14 @@ #include "Acts/MagneticField/SolenoidBField.hpp" #include "Acts/Utilities/Axis.hpp" #include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Utilities/VectorHelpers.hpp" #include "Acts/Utilities/detail/grid_helper.hpp" #include #include +#include #include #include #include @@ -35,100 +37,73 @@ Acts::fieldMapRZ( std::array nBinsRZ)>& localToGlobalBin, std::vector rPos, std::vector zPos, - std::vector bField, double lengthUnit, double BFieldUnit, + const std::vector& bField, double lengthUnit, double BFieldUnit, bool firstQuadrant) { // [1] Create Grid - // sort the values - std::ranges::sort(rPos); - std::ranges::sort(zPos); - // Get unique values - rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - rPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsR = rPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum. We just sorted the vectors, so these are just - // the first and last elements. - double rMin = rPos[0]; - double zMin = zPos[0]; - double rMax = rPos[nBinsR - 1]; - double zMax = zPos[nBinsZ - 1]; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; - zMax += stepZ; + const auto [rMin, rMax, rBinCount] = detail::getMinMaxAndBinCount(rPos); + auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos); + + const std::size_t nBinsR = rBinCount; + std::size_t nBinsZ = zBinCount; + if (firstQuadrant) { zMin = -zPos[nBinsZ - 1]; - nBinsZ = static_cast(2. * nBinsZ - 1); + nBinsZ = 2 * nBinsZ - 1; } // Create the axis for the grid - Acts::Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); - Acts::Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); + Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); + Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); // Create the grid - Grid grid(Type, std::move(rAxis), std::move(zAxis)); + Grid grid(Type, std::move(rAxis), std::move(zAxis)); using Grid_t = decltype(grid); // [2] Set the bField values + const std::array nIndices = {{rBinCount, zBinCount}}; for (std::size_t i = 1; i <= nBinsR; ++i) { for (std::size_t j = 1; j <= nBinsZ; ++j) { - std::array nIndices = {{rPos.size(), zPos.size()}}; Grid_t::index_t indices = {{i, j}}; + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one if (firstQuadrant) { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one - std::size_t n = - std::abs(static_cast(j) - static_cast(zPos.size())); - Grid_t::index_t indicesFirstQuadrant = {{i - 1, n}}; + std::size_t n = std::abs(static_cast(j) - + static_cast(zBinCount)); grid.atLocalBins(indices) = - bField.at(localToGlobalBin(indicesFirstQuadrant, nIndices)) * - BFieldUnit; + bField.at(localToGlobalBin({{i - 1, n}}, nIndices)) * BFieldUnit; } else { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one grid.atLocalBins(indices) = bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices)) * BFieldUnit; } } } - grid.setExteriorBins(Acts::Vector2::Zero()); + grid.setExteriorBins(Vector2::Zero()); - // [3] Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { - return Acts::Vector2(perp(pos), pos.z()); + // [3] Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { + return Vector2(perp(pos), pos.z()); }; - // [4] Create the transformation for the bfield - // map (Br,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector2& field, - const Acts::Vector3& pos) { - double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y(); - double cos_phi = 0, sin_phi = 0; - if (r_sin_theta_2 > std::numeric_limits::min()) { - double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2); - cos_phi = pos.x() * inv_r_sin_theta; - sin_phi = pos.y() * inv_r_sin_theta; - } else { - cos_phi = 1.; - sin_phi = 0.; + // [4] Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector2& field, const Vector3& pos) { + const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y(); + double cosPhi = 1.; + double sinPhi = 0.; + + if (rSinTheta2 > std::numeric_limits::min()) { + const double invRsinTheta = 1. / std::sqrt(rSinTheta2); + cosPhi = pos.x() * invRsinTheta; + sinPhi = pos.y() * invRsinTheta; } - return Acts::Vector3(field.x() * cos_phi, field.x() * sin_phi, field.y()); + + return Vector3(field.x() * cosPhi, field.x() * sinPhi, field.y()); }; - // [5] Create the mapper & BField Service - // create field mapping - return Acts::InterpolatedBFieldMap( + // [5] Create the mapper & BField Service create field mapping + return InterpolatedBFieldMap( {transformPos, transformBField, std::move(grid)}); } @@ -141,88 +116,58 @@ Acts::fieldMapXYZ( std::array nBinsXYZ)>& localToGlobalBin, std::vector xPos, std::vector yPos, - std::vector zPos, std::vector bField, + std::vector zPos, const std::vector& bField, double lengthUnit, double BFieldUnit, bool firstOctant) { // [1] Create Grid - // Sort the values - std::ranges::sort(xPos); - std::ranges::sort(yPos); - std::ranges::sort(zPos); - // Get unique values - xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end()); - yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - xPos.shrink_to_fit(); - yPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsX = xPos.size(); - std::size_t nBinsY = yPos.size(); - std::size_t nBinsZ = zPos.size(); + auto [xMin, xMax, xBinCount] = detail::getMinMaxAndBinCount(xPos); + auto [yMin, yMax, yBinCount] = detail::getMinMaxAndBinCount(yPos); + auto [zMin, zMax, zBinCount] = detail::getMinMaxAndBinCount(zPos); - // Create the axis for the grid - // get minima and maximia. We just sorted the vectors, so these are just the - // first and last elements. - double xMin = xPos[0]; - double yMin = yPos[0]; - double zMin = zPos[0]; - // get maxima - double xMax = xPos[nBinsX - 1]; - double yMax = yPos[nBinsY - 1]; - double zMax = zPos[nBinsZ - 1]; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::abs(yMax - yMin) / (nBinsY - 1); - double stepX = std::abs(xMax - xMin) / (nBinsX - 1); - xMax += stepX; - yMax += stepY; - zMax += stepZ; + std::size_t nBinsX = xBinCount; + std::size_t nBinsY = yBinCount; + std::size_t nBinsZ = zBinCount; - // If only the first octant is given if (firstOctant) { xMin = -xPos[nBinsX - 1]; - yMin = -yPos[nBinsY - 1]; - zMin = -zPos[nBinsZ - 1]; nBinsX = 2 * nBinsX - 1; + yMin = -yPos[nBinsY - 1]; nBinsY = 2 * nBinsY - 1; + zMin = -zPos[nBinsZ - 1]; nBinsZ = 2 * nBinsZ - 1; } - Acts::Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); - Acts::Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); - Acts::Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); + + Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); + Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); + Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); // Create the grid Grid grid(Type, std::move(xAxis), std::move(yAxis), std::move(zAxis)); using Grid_t = decltype(grid); // [2] Set the bField values + const std::array nIndices = { + {xBinCount, yBinCount, zBinCount}}; + + auto calcAbsDiff = [](std::size_t val, std::size_t binCount) { + return std::abs(static_cast(val) - + static_cast(binCount)); + }; + for (std::size_t i = 1; i <= nBinsX; ++i) { for (std::size_t j = 1; j <= nBinsY; ++j) { for (std::size_t k = 1; k <= nBinsZ; ++k) { Grid_t::index_t indices = {{i, j, k}}; - std::array nIndices = { - {xPos.size(), yPos.size(), zPos.size()}}; + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one if (firstOctant) { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one - std::size_t m = - std::abs(static_cast(i) - (static_cast(xPos.size()))); - std::size_t n = - std::abs(static_cast(j) - (static_cast(yPos.size()))); - std::size_t l = - std::abs(static_cast(k) - (static_cast(zPos.size()))); - Grid_t::index_t indicesFirstOctant = {{m, n, l}}; + const std::size_t l = calcAbsDiff(i, xBinCount); + const std::size_t m = calcAbsDiff(j, yBinCount); + const std::size_t n = calcAbsDiff(k, zBinCount); grid.atLocalBins(indices) = - bField.at(localToGlobalBin(indicesFirstOctant, nIndices)) * - BFieldUnit; - + bField.at(localToGlobalBin({{l, m, n}}, nIndices)) * BFieldUnit; } else { - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one grid.atLocalBins(indices) = bField.at(localToGlobalBin({{i - 1, j - 1, k - 1}}, nIndices)) * BFieldUnit; @@ -230,74 +175,67 @@ Acts::fieldMapXYZ( } } } - grid.setExteriorBins(Acts::Vector3::Zero()); + grid.setExteriorBins(Vector3::Zero()); - // [3] Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { return pos; }; + // [3] Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { return pos; }; - // [4] Create the transformation for the bfield - // map (Bx,By,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector3& field, - const Acts::Vector3& /*pos*/) { return field; }; + // [4] Create the transformation for the BField map (Bx,By,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector3& field, const Vector3& /*pos*/) { + return field; + }; - // [5] Create the mapper & BField Service - // create field mapping - return Acts::InterpolatedBFieldMap( + // [5] Create the mapper & BField Service create field mapping + return InterpolatedBFieldMap( {transformPos, transformBField, std::move(grid)}); } Acts::InterpolatedBFieldMap< Acts::Grid, Acts::Axis>> -Acts::solenoidFieldMap(std::pair rlim, - std::pair zlim, - std::pair nbins, +Acts::solenoidFieldMap(const std::pair& rLim, + const std::pair& zLim, + const std::pair& nBins, const SolenoidBField& field) { - auto [rMin, rMax] = rlim; - auto [zMin, zMax] = zlim; - const auto [nBinsR, nBinsZ] = nbins; + auto [rMin, rMax] = rLim; + auto [zMin, zMax] = zLim; + const auto [nBinsR, nBinsZ] = nBins; double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; zMax += stepZ; // Create the axis for the grid - Acts::Axis rAxis(rMin, rMax, nBinsR); - Acts::Axis zAxis(zMin, zMax, nBinsZ); + Axis rAxis(rMin, rMax, nBinsR); + Axis zAxis(zMin, zMax, nBinsZ); // Create the grid - Grid grid(Type, std::move(rAxis), std::move(zAxis)); + Grid grid(Type, std::move(rAxis), std::move(zAxis)); using Grid_t = decltype(grid); - // Create the transformation for the position - // map (x,y,z) -> (r,z) - auto transformPos = [](const Acts::Vector3& pos) { - return Acts::Vector2(perp(pos), pos.z()); + // Create the transformation for the position map (x,y,z) -> (r,z) + auto transformPos = [](const Vector3& pos) { + return Vector2(perp(pos), pos.z()); }; - // Create the transformation for the bfield - // map (Br,Bz) -> (Bx,By,Bz) - auto transformBField = [](const Acts::Vector2& bfield, - const Acts::Vector3& pos) { - double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y(); - double cos_phi = 0, sin_phi = 0; - if (r_sin_theta_2 > std::numeric_limits::min()) { - double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2); - cos_phi = pos.x() * inv_r_sin_theta; - sin_phi = pos.y() * inv_r_sin_theta; - } else { - cos_phi = 1.; - sin_phi = 0.; + // Create the transformation for the bField map (Br,Bz) -> (Bx,By,Bz) + auto transformBField = [](const Vector2& bField, const Vector3& pos) { + const double rSinTheta2 = pos.x() * pos.x() + pos.y() * pos.y(); + double cosPhi = 1.; + double sinPhi = 0.; + + if (rSinTheta2 > std::numeric_limits::min()) { + const double invRsinTheta = 1. / std::sqrt(rSinTheta2); + cosPhi = pos.x() * invRsinTheta; + sinPhi = pos.y() * invRsinTheta; } - return Acts::Vector3(bfield.x() * cos_phi, bfield.x() * sin_phi, - bfield.y()); + + return Vector3(bField.x() * cosPhi, bField.x() * sinPhi, bField.y()); }; - // iterate over all bins, set their value to the solenoid value - // at their lower left position + // iterate over all bins, set their value to the solenoid value at their lower + // left position for (std::size_t i = 0; i <= nBinsR + 1; i++) { for (std::size_t j = 0; j <= nBinsZ + 1; j++) { Grid_t::index_t index({i, j}); @@ -314,9 +252,8 @@ Acts::solenoidFieldMap(std::pair rlim, } } - // Create the mapper & BField Service - // create field mapping - Acts::InterpolatedBFieldMap map( + // Create the mapper & BField Service create field mapping + InterpolatedBFieldMap map( {transformPos, transformBField, std::move(grid)}); return map; } diff --git a/Core/src/Material/MaterialMapUtils.cpp b/Core/src/Material/MaterialMapUtils.cpp index b51b1a3e223..8f4d23cd71e 100644 --- a/Core/src/Material/MaterialMapUtils.cpp +++ b/Core/src/Material/MaterialMapUtils.cpp @@ -12,6 +12,7 @@ #include "Acts/Material/Material.hpp" #include "Acts/Utilities/Axis.hpp" #include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Helpers.hpp" #include #include @@ -43,31 +44,8 @@ auto Acts::materialMapperRZ( } // [2] Create Grid - // sort the values - std::ranges::sort(rPos); - std::ranges::sort(zPos); - // Get unique values - rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - rPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsR = rPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum - auto minMaxR = std::minmax_element(rPos.begin(), rPos.end()); - auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end()); - double rMin = *minMaxR.first; - double zMin = *minMaxZ.first; - double rMax = *minMaxR.second; - double zMax = *minMaxZ.second; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepR = std::abs(rMax - rMin) / (nBinsR - 1); - rMax += stepR; - zMax += stepZ; + const auto [rMin, rMax, nBinsR] = detail::getMinMaxAndBinCount(rPos); + const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos); // Create the axis for the grid Axis rAxis(rMin * lengthUnit, rMax * lengthUnit, nBinsR); @@ -79,13 +57,13 @@ auto Acts::materialMapperRZ( using Grid_t = decltype(grid); // [3] Set the material values + const std::array nIndices = {{nBinsR, nBinsZ}}; for (std::size_t i = 1; i <= nBinsR; ++i) { for (std::size_t j = 1; j <= nBinsZ; ++j) { - std::array nIndices = {{rPos.size(), zPos.size()}}; Grid_t::index_t indices = {{i, j}}; - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one grid.atLocalBins(indices) = materialVector.at( materialVectorToGridMapper({{i - 1, j - 1}}, nIndices)); } @@ -95,14 +73,12 @@ auto Acts::materialMapperRZ( 0., 0., 0.; grid.setExteriorBins(vec); - // [4] Create the transformation for the position - // map (x,y,z) -> (r,z) + // [4] Create the transformation for the position map (x,y,z) -> (r,z) auto transformPos = [](const Vector3& pos) { return Vector2(perp(pos), pos.z()); }; - // [5] Create the mapper & BField Service - // create material mapping + // [5] Create the mapper & BField Service create material mapping return MaterialMapper(transformPos, std::move(grid)); } @@ -125,44 +101,11 @@ auto Acts::materialMapperXYZ( } // [2] Create Grid - // Sort the values - std::ranges::sort(xPos); - std::ranges::sort(yPos); - std::ranges::sort(zPos); - // Get unique values - xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end()); - yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end()); - zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end()); - xPos.shrink_to_fit(); - yPos.shrink_to_fit(); - zPos.shrink_to_fit(); - // get the number of bins - std::size_t nBinsX = xPos.size(); - std::size_t nBinsY = yPos.size(); - std::size_t nBinsZ = zPos.size(); - - // get the minimum and maximum - auto minMaxX = std::minmax_element(xPos.begin(), xPos.end()); - auto minMaxY = std::minmax_element(yPos.begin(), yPos.end()); - auto minMaxZ = std::minmax_element(zPos.begin(), zPos.end()); - // Create the axis for the grid - // get minima - double xMin = *minMaxX.first; - double yMin = *minMaxY.first; - double zMin = *minMaxZ.first; - // get maxima - double xMax = *minMaxX.second; - double yMax = *minMaxY.second; - double zMax = *minMaxZ.second; - // calculate maxima (add one last bin, because bin value always corresponds to - // left boundary) - double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1); - double stepY = std::abs(yMax - yMin) / (nBinsY - 1); - double stepX = std::abs(xMax - xMin) / (nBinsX - 1); - xMax += stepX; - yMax += stepY; - zMax += stepZ; + const auto [xMin, xMax, nBinsX] = detail::getMinMaxAndBinCount(xPos); + const auto [yMin, yMax, nBinsY] = detail::getMinMaxAndBinCount(yPos); + const auto [zMin, zMax, nBinsZ] = detail::getMinMaxAndBinCount(zPos); + // Create the axis for the grid Axis xAxis(xMin * lengthUnit, xMax * lengthUnit, nBinsX); Axis yAxis(yMin * lengthUnit, yMax * lengthUnit, nBinsY); Axis zAxis(zMin * lengthUnit, zMax * lengthUnit, nBinsZ); @@ -172,15 +115,14 @@ auto Acts::materialMapperXYZ( using Grid_t = decltype(grid); // [3] Set the bField values + const std::array nIndices = {{nBinsX, nBinsY, nBinsZ}}; for (std::size_t i = 1; i <= nBinsX; ++i) { for (std::size_t j = 1; j <= nBinsY; ++j) { for (std::size_t k = 1; k <= nBinsZ; ++k) { Grid_t::index_t indices = {{i, j, k}}; - std::array nIndices = { - {xPos.size(), yPos.size(), zPos.size()}}; - // std::vectors begin with 0 and we do not want the user needing to - // take underflow or overflow bins in account this is why we need to - // subtract by one + // std::vectors begin with 0 and we do not want the user needing to take + // underflow or overflow bins in account this is why we need to subtract + // by one grid.atLocalBins(indices) = materialVector.at( materialVectorToGridMapper({{i - 1, j - 1, k - 1}}, nIndices)); } @@ -191,11 +133,9 @@ auto Acts::materialMapperXYZ( 0., 0., 0.; grid.setExteriorBins(vec); - // [4] Create the transformation for the position - // map (x,y,z) -> (r,z) + // [4] Create the transformation for the position map (x,y,z) -> (r,z) auto transformPos = [](const Vector3& pos) { return pos; }; - // [5] Create the mapper & BField Service - // create material mapping + // [5] Create the mapper & BField Service create material mapping return MaterialMapper(transformPos, std::move(grid)); } From 4a506f8298d13376c5c85335e3783c1e7aed43c8 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Mon, 2 Dec 2024 15:51:03 +0100 Subject: [PATCH 7/8] refactor: replace `long unsigned int` with `std::size_t` (#3930) ## Summary by CodeRabbit - **New Features** - Enhanced error handling for input space points in the GbtsSeedingAlgorithm. - Improved seed finding logic based on internal ROI descriptors. - **Bug Fixes** - Updated handling of mapping logic to ensure clarity and prevent missing keys. - **Refactor** - Changed variable types from `long unsigned int` to `std::size_t` across multiple components for improved type safety. - **Tests** - Updated loop index types in various test cases to enhance type safety and align with best practices. --- Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp | 2 +- .../include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp | 2 +- Tests/UnitTests/Plugins/GeoModel/GeoBoxToVolumeTest.cpp | 2 +- Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp | 4 ++-- Tests/UnitTests/Plugins/GeoModel/GeoPolyConverterTests.cpp | 4 ++-- Tests/UnitTests/Plugins/GeoModel/GeoTubeToVolumeTest.cpp | 2 +- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp index 677ba2b4c10..c126225aea1 100644 --- a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp @@ -355,7 +355,7 @@ ActsExamples::GbtsSeedingAlgorithm::LayerNumbering() const { } }); - for (long unsigned int i = 0; i < input_vector.size(); i++) { + for (std::size_t i = 0; i < input_vector.size(); i++) { input_vector[i].m_refCoord = input_vector[i].m_refCoord / count_vector[i]; } diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp index 2a09d8492b2..f7873bd0e2a 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp @@ -172,7 +172,7 @@ class RootAthenaDumpReader : public IReader { /// Vector of {eventNr, entryMin, entryMax} std::vector> m_eventMap; std::shared_ptr m_inputchain; - long unsigned int m_events; + std::size_t m_events; bool m_haveStripFeatures = true; static constexpr unsigned int maxCL = 1500000; diff --git a/Tests/UnitTests/Plugins/GeoModel/GeoBoxToVolumeTest.cpp b/Tests/UnitTests/Plugins/GeoModel/GeoBoxToVolumeTest.cpp index b905fda65ee..426b7cc07c3 100644 --- a/Tests/UnitTests/Plugins/GeoModel/GeoBoxToVolumeTest.cpp +++ b/Tests/UnitTests/Plugins/GeoModel/GeoBoxToVolumeTest.cpp @@ -58,7 +58,7 @@ BOOST_AUTO_TEST_CASE(GeoBoxToSensitiveConversion) { const auto* bounds = dynamic_cast(&volumeBox->volumeBounds()); std::vector convHls = bounds->values(); - for (long unsigned int i = 0; i < hls.size(); i++) { + for (std::size_t i = 0; i < hls.size(); i++) { BOOST_CHECK(hls[i] == convHls[i]); } } diff --git a/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp b/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp index 033f2c9cb4a..ae3c4ec0420 100644 --- a/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp +++ b/Tests/UnitTests/Plugins/GeoModel/GeoDetectorObjectTest.cpp @@ -43,7 +43,7 @@ void test(const Acts::GeoModelDetectorObjectFactory::Cache& cache, GeoModelDetObj::GeoDims geoDims) { for (const auto& box : cache.boundingBoxes) { const Acts::VolumeBounds& bounds = box->volumeBounds(); - for (long unsigned int i = 0; i < geoDims.boxO.size(); i++) { + for (std::size_t i = 0; i < geoDims.boxO.size(); i++) { BOOST_CHECK(geoDims.boxO[i] == bounds.values()[i]); } std::vector surfaces = box->surfaces(); @@ -75,7 +75,7 @@ void test(const Acts::GeoModelDetectorObjectFactory::Cache& cache, dynamic_cast(&sbounds); std::vector trapVerts = trapBounds->vertices(); - for (long unsigned int i = 0; i < trapVerts.size(); i++) { + for (std::size_t i = 0; i < trapVerts.size(); i++) { BOOST_CHECK(trapVerts[i][0] == geoDims.trapVerts[i][0]); BOOST_CHECK(trapVerts[i][1] == geoDims.trapVerts[i][1]); } diff --git a/Tests/UnitTests/Plugins/GeoModel/GeoPolyConverterTests.cpp b/Tests/UnitTests/Plugins/GeoModel/GeoPolyConverterTests.cpp index b273b0a2364..a9ed6fa6f56 100644 --- a/Tests/UnitTests/Plugins/GeoModel/GeoPolyConverterTests.cpp +++ b/Tests/UnitTests/Plugins/GeoModel/GeoPolyConverterTests.cpp @@ -88,7 +88,7 @@ BOOST_AUTO_TEST_CASE(GeoModelDetectorObjectFactory) { const auto* polyBounds = dynamic_cast(&polySurface->bounds()); std::vector convPolyVerts = polyBounds->vertices(); - for (long unsigned int i = 0; i < polyVerts.size(); i++) { + for (std::size_t i = 0; i < polyVerts.size(); i++) { BOOST_CHECK(polyVerts[i][0] == convPolyVerts[i][0]); BOOST_CHECK(polyVerts[i][1] == convPolyVerts[i][1]); } @@ -96,7 +96,7 @@ BOOST_AUTO_TEST_CASE(GeoModelDetectorObjectFactory) { const auto* trapBounds = dynamic_cast(&trapSurface->bounds()); std::vector convTrapVerts = trapBounds->vertices(); - for (long unsigned int i = 0; i < trapVerts.size(); i++) { + for (std::size_t i = 0; i < trapVerts.size(); i++) { BOOST_CHECK(trapVerts[i][0] == convTrapVerts[i][0]); BOOST_CHECK(trapVerts[i][1] == convTrapVerts[i][1]); } diff --git a/Tests/UnitTests/Plugins/GeoModel/GeoTubeToVolumeTest.cpp b/Tests/UnitTests/Plugins/GeoModel/GeoTubeToVolumeTest.cpp index 64d368e7f2b..fd58df2f7ff 100644 --- a/Tests/UnitTests/Plugins/GeoModel/GeoTubeToVolumeTest.cpp +++ b/Tests/UnitTests/Plugins/GeoModel/GeoTubeToVolumeTest.cpp @@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE(GeoBoxToSensitiveConversion) { const auto* bounds = dynamic_cast( &volumeTube->volumeBounds()); std::vector convDims = bounds->values(); - for (long unsigned int i = 0; i < dims.size(); i++) { + for (std::size_t i = 0; i < dims.size(); i++) { BOOST_CHECK(dims[i] == convDims[i]); } } From 6cbf203aa3703d8484e501a3c119e1c3dca32c8b Mon Sep 17 00:00:00 2001 From: Stephen Nicholas Swatman Date: Mon, 2 Dec 2024 18:45:58 +0100 Subject: [PATCH 8/8] fix: Explicitly disable paging in pre-commit hook (#3928) The pre-commit hook that checks for leftover conflict markers introduced in #3708 calls git diff. While this works well in general, I am running into an edge case where git is failing to determine that it is not being run in a TTY, causing it to fire up a pager which, in turn, causes pre-commit to hang as it waits for the pager to close. This commit fixes this issue by adding the `--no-pager` flag to the git command. --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 68c78341c5d..4158678a3e0 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -92,7 +92,7 @@ repos: - id: leftover_conflict_markers name: Leftover conflict markers language: system - entry: git diff --staged --check + entry: git --no-pager diff --staged --check - repo: local hooks: