From 3a8a0e365db45af5f251c70f008cb71becf64c65 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 8 Aug 2024 10:07:33 +0200 Subject: [PATCH] feat: Variable size measurement for Examples (#3107) This is an attempt to implement a variable size measurement as an alternative to the fixed size measurement. Currently our Examples framework is using a variant over fixed size measurements to cope with measurements of different dimensions. This obfuscates the code a bit since accessing the actual measurement requires a visitor pattern to acquire the dimensional typed measurement. Here I propose to switch this out for a variable size measurement which allocates the same amount of memory but stores the dimension as a direct member instead of pushing this detail to `std::variant`. blocked by - https://github.com/acts-project/acts/pull/3192 - https://github.com/acts-project/acts/pull/3331 --- .../Acts/EventData/MultiTrajectory.hpp | 16 ++ .../Acts/Utilities/detail/Subspace.hpp | 52 ++-- .../Digitization/src/MeasurementCreation.cpp | 12 +- .../TrackFinding/src/HoughTransformSeeder.cpp | 30 +-- .../TrackFinding/src/SpacePointMaker.cpp | 10 +- .../Utilities/src/PrototracksToTracks.cpp | 11 +- .../Framework/ML/src/NeuralCalibrator.cpp | 162 +++++------ .../ActsExamples/EventData/Measurement.hpp | 254 ++++++++---------- .../src/EventData/MeasurementCalibration.cpp | 22 +- .../src/EventData/ScalingCalibrator.cpp | 56 ++-- Examples/Io/Csv/src/CsvMeasurementWriter.cpp | 97 ++++--- Examples/Io/EDM4hep/src/EDM4hepUtil.cpp | 90 +++---- .../Io/Root/src/RootMeasurementWriter.cpp | 86 +++--- .../Core/Utilities/SubspaceTests.cpp | 12 +- .../Examples/EventData/MeasurementTests.cpp | 197 +++----------- .../Io/Csv/MeasurementReaderWriterTests.cpp | 22 +- 16 files changed, 461 insertions(+), 668 deletions(-) diff --git a/Core/include/Acts/EventData/MultiTrajectory.hpp b/Core/include/Acts/EventData/MultiTrajectory.hpp index b914556b77b..00eea920dcc 100644 --- a/Core/include/Acts/EventData/MultiTrajectory.hpp +++ b/Core/include/Acts/EventData/MultiTrajectory.hpp @@ -499,6 +499,10 @@ class MultiTrajectory { return self().template calibratedCovariance_impl(istate); } + /// Retrieve a calibrated measurement covariance proxy instance for a + /// measurement at a given index + /// @param istate The track state + /// @return Mutable proxy typename TrackStateProxy::EffectiveCalibrated effectiveCalibrated( IndexType istate) requires(!ReadOnly) { // This abuses an incorrectly sized vector / matrix to access the @@ -508,6 +512,10 @@ class MultiTrajectory { calibrated(istate).data(), calibratedSize(istate)}; } + /// Retrieve a calibrated measurement covariance proxy instance for a + /// measurement at a given index + /// @param istate The track state + /// @return Const proxy typename ConstTrackStateProxy::EffectiveCalibrated effectiveCalibrated( IndexType istate) const { // This abuses an incorrectly sized vector / matrix to access the @@ -517,6 +525,10 @@ class MultiTrajectory { calibrated(istate).data(), calibratedSize(istate)}; } + /// Retrieve a calibrated measurement covariance proxy instance for a + /// measurement at a given index + /// @param istate The track state + /// @return Mutable proxy typename TrackStateProxy::EffectiveCalibratedCovariance effectiveCalibratedCovariance(IndexType istate) requires(!ReadOnly) { // This abuses an incorrectly sized vector / matrix to access the @@ -527,6 +539,10 @@ class MultiTrajectory { calibratedSize(istate)}; } + /// Retrieve a calibrated measurement covariance proxy instance for a + /// measurement at a given index + /// @param istate The track state + /// @return Const proxy typename ConstTrackStateProxy::EffectiveCalibratedCovariance effectiveCalibratedCovariance(IndexType istate) const { // This abuses an incorrectly sized vector / matrix to access the diff --git a/Core/include/Acts/Utilities/detail/Subspace.hpp b/Core/include/Acts/Utilities/detail/Subspace.hpp index 556a4bc4b2f..fe355b44edb 100644 --- a/Core/include/Acts/Utilities/detail/Subspace.hpp +++ b/Core/include/Acts/Utilities/detail/Subspace.hpp @@ -119,6 +119,15 @@ class FixedSizeSubspace { /// @return Axis index in the full space constexpr std::size_t operator[](std::size_t i) const { return m_axes[i]; } + std::size_t indexOf(std::size_t axis) const { + for (std::size_t i = 0; i < kSize; ++i) { + if (m_axes[i] == axis) { + return i; + } + } + return kSize; + } + /// Axis indices that comprise the subspace. /// /// The specific container and index type should be considered an @@ -266,6 +275,15 @@ class VariableSizeSubspace { return m_axes[i]; } + std::size_t indexOf(std::size_t axis) const { + for (std::size_t i = 0; i < m_size; ++i) { + if (m_axes[i] == axis) { + return i; + } + } + return m_size; + } + /// Check if the given axis index in the full space is part of the subspace. constexpr bool contains(std::size_t index) const { bool isContained = false; @@ -302,40 +320,6 @@ class VariableSizeSubspace { } return expn; } - - std::uint64_t projectorBits() const { - std::uint64_t result = 0; - - for (std::size_t i = 0; i < m_size; ++i) { - for (std::size_t j = 0; j < kFullSize; ++j) { - // the bit order is defined in `Acts/Utilities/AlgebraHelpers.hpp` - // in `matrixToBitset` - std::size_t index = m_size * kFullSize - 1 - (i + j * m_size); - if (m_axes[i] == j) { - result |= (1ull << index); - } - } - } - - return result; - } - - std::uint64_t fullProjectorBits() const { - std::uint64_t result = 0; - - for (std::size_t i = 0; i < kFullSize; ++i) { - for (std::size_t j = 0; j < kFullSize; ++j) { - // the bit order is defined in `Acts/Utilities/AlgebraHelpers.hpp` - // in `matrixToBitset` - std::size_t index = kFullSize * kFullSize - 1 - (i + j * kFullSize); - if (i < m_size && m_axes[i] == j) { - result |= (1ull << index); - } - } - } - - return result; - } }; /// @} diff --git a/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp b/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp index 318174589e0..7338b431f8c 100644 --- a/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp +++ b/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp @@ -21,23 +21,19 @@ ActsExamples::Measurement ActsExamples::createMeasurement( switch (dParams.indices.size()) { case 1u: { auto [indices, par, cov] = measurementConstituents<1>(dParams); - return FixedSizeMeasurement(std::move(sl), indices, - par, cov); + return ActsExamples::Measurement(std::move(sl), indices, par, cov); } case 2u: { auto [indices, par, cov] = measurementConstituents<2>(dParams); - return FixedSizeMeasurement(std::move(sl), indices, - par, cov); + return ActsExamples::Measurement(std::move(sl), indices, par, cov); }; case 3u: { auto [indices, par, cov] = measurementConstituents<3>(dParams); - return FixedSizeMeasurement(std::move(sl), indices, - par, cov); + return ActsExamples::Measurement(std::move(sl), indices, par, cov); }; case 4u: { auto [indices, par, cov] = measurementConstituents<4>(dParams); - return FixedSizeMeasurement(std::move(sl), indices, - par, cov); + return ActsExamples::Measurement(std::move(sl), indices, par, cov); }; default: std::string errorMsg = "Invalid/mismatching measurement dimension: " + diff --git a/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp b/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp index 9333d7b11a8..f458b33795a 100644 --- a/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp +++ b/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp @@ -525,20 +525,18 @@ void ActsExamples::HoughTransformSeeder::addMeasurements( // are transformed to the bound space where we do know their location. // if the local parameters are not measured, this results in a // zero location, which is a reasonable default fall-back. - auto [localPos, localCov] = std::visit( - [](const auto& meas) { - auto expander = meas.expander(); - Acts::BoundVector par = expander * meas.parameters(); - Acts::BoundSquareMatrix cov = - expander * meas.covariance() * expander.transpose(); - // extract local position - Acts::Vector2 lpar(par[Acts::eBoundLoc0], par[Acts::eBoundLoc1]); - // extract local position covariance. - Acts::SquareMatrix2 lcov = - cov.block<2, 2>(Acts::eBoundLoc0, Acts::eBoundLoc0); - return std::make_pair(lpar, lcov); - }, - measurements[sourceLink.index()]); + const auto& measurement = measurements[sourceLink.index()]; + + assert(measurement.contains(Acts::eBoundLoc0) && + "Measurement does not contain the required bound loc0"); + assert(measurement.contains(Acts::eBoundLoc1) && + "Measurement does not contain the required bound loc1"); + + auto boundLoc0 = measurement.subspace().indexOf(Acts::eBoundLoc0); + auto boundLoc1 = measurement.subspace().indexOf(Acts::eBoundLoc1); + + Acts::Vector2 localPos{measurement.effectiveParameters()[boundLoc0], + measurement.effectiveParameters()[boundLoc1]}; // transform local position to global coordinates Acts::Vector3 globalFakeMom(1, 1, 1); @@ -551,10 +549,10 @@ void ActsExamples::HoughTransformSeeder::addMeasurements( if (hitlayer.ok()) { std::vector index; index.push_back(sourceLink.index()); - auto meas = std::shared_ptr( + auto houghMeas = std::shared_ptr( new HoughMeasurementStruct(hitlayer.value(), phi, r, z, index, HoughHitType::MEASUREMENT)); - houghMeasurementStructs.push_back(meas); + houghMeasurementStructs.push_back(houghMeas); } } } diff --git a/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp b/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp index 4cf3205f2b0..49802b57087 100644 --- a/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp +++ b/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp @@ -128,15 +128,7 @@ ActsExamples::ProcessCode ActsExamples::SpacePointMaker::execute( const auto islink = slink.get(); const auto& meas = measurements[islink.index()]; - return std::visit( - [](const auto& measurement) { - auto expander = measurement.expander(); - Acts::BoundVector par = expander * measurement.parameters(); - Acts::BoundSquareMatrix cov = - expander * measurement.covariance() * expander.transpose(); - return std::make_pair(par, cov); - }, - meas); + return std::make_pair(meas.fullParameters(), meas.fullCovariance()); }; SimSpacePointContainer spacePoints; diff --git a/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp b/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp index d3d0156ec51..9506bc9b1e8 100644 --- a/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp +++ b/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp @@ -32,14 +32,9 @@ ProcessCode PrototracksToTracks::execute(const AlgorithmContext& ctx) const { TrackContainer tracks(trackContainer, mtj); boost::container::flat_map slMap; - for (const auto& varm : m_inputMeasurements(ctx)) { - std::visit( - [&](const auto& m) { - const auto idx = - m.sourceLink().template get().index(); - slMap.insert(std::pair{idx, m.sourceLink()}); - }, - varm); + for (const auto& m : m_inputMeasurements(ctx)) { + const auto idx = m.sourceLink().template get().index(); + slMap.insert(std::pair{idx, m.sourceLink()}); } const auto& prototracks = m_inputProtoTracks(ctx); diff --git a/Examples/Framework/ML/src/NeuralCalibrator.cpp b/Examples/Framework/ML/src/NeuralCalibrator.cpp index 29be5a07ad9..6dcbe1cc7cf 100644 --- a/Examples/Framework/ML/src/NeuralCalibrator.cpp +++ b/Examples/Framework/ML/src/NeuralCalibrator.cpp @@ -8,6 +8,8 @@ #include "ActsExamples/EventData/NeuralCalibrator.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/MeasurementHelpers.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/Utilities/CalibrationContext.hpp" #include "Acts/Utilities/UnitVectors.hpp" @@ -102,81 +104,89 @@ void ActsExamples::NeuralCalibrator::calibrate( input[iInput++] = idxSourceLink.geometryId().layer(); const Acts::Surface& referenceSurface = trackState.referenceSurface(); + auto trackParameters = trackState.parameters(); + + const auto& measurement = measurements[idxSourceLink.index()]; + + assert(measurement.contains(Acts::eBoundLoc0) && + "Measurement does not contain the required bound loc0"); + assert(measurement.contains(Acts::eBoundLoc1) && + "Measurement does not contain the required bound loc1"); + + auto boundLoc0 = measurement.subspace().indexOf(Acts::eBoundLoc0); + auto boundLoc1 = measurement.subspace().indexOf(Acts::eBoundLoc1); + + Acts::Vector2 localPosition{measurement.effectiveParameters()[boundLoc0], + measurement.effectiveParameters()[boundLoc1]}; + Acts::Vector2 localCovariance{ + measurement.effectiveCovariance()(boundLoc0, boundLoc0), + measurement.effectiveCovariance()(boundLoc1, boundLoc1)}; + + Acts::Vector3 dir = Acts::makeDirectionFromPhiTheta( + trackParameters[Acts::eBoundPhi], trackParameters[Acts::eBoundTheta]); + Acts::Vector3 globalPosition = + referenceSurface.localToGlobal(gctx, localPosition, dir); + + // Rotation matrix. When applied to global coordinates, they + // are rotated into the local reference frame of the + // surface. Note that this such a rotation can be found by + // inverting a matrix whose columns correspond to the + // coordinate axes of the local coordinate system. + Acts::RotationMatrix3 rot = + referenceSurface.referenceFrame(gctx, globalPosition, dir).inverse(); + std::pair angles = + Acts::VectorHelpers::incidentAngles(dir, rot); + + input[iInput++] = angles.first; + input[iInput++] = angles.second; + input[iInput++] = localPosition[0]; + input[iInput++] = localPosition[1]; + input[iInput++] = localCovariance[0]; + input[iInput++] = localCovariance[1]; + if (iInput != m_nInputs) { + throw std::runtime_error("Expected input size of " + + std::to_string(m_nInputs) + + ", got: " + std::to_string(iInput)); + } + + // Input is a single row, hence .front() + std::vector output = m_model.runONNXInference(inputBatch).front(); + // Assuming 2-D measurements, the expected params structure is: + // [ 0, nComponent[ --> priors + // [ nComponent, 3*nComponent[ --> means + // [3*nComponent, 5*nComponent[ --> variances + std::size_t nParams = 5 * m_nComponents; + if (output.size() != nParams) { + throw std::runtime_error("Got output vector of size " + + std::to_string(output.size()) + + ", expected size " + std::to_string(nParams)); + } - std::visit( - [&](const auto& measurement) { - auto E = measurement.expander(); - auto P = measurement.projector(); - Acts::ActsVector fpar = E * measurement.parameters(); - Acts::ActsSquareMatrix fcov = - E * measurement.covariance() * E.transpose(); - - Acts::Vector3 dir = Acts::makeDirectionFromPhiTheta( - fpar[Acts::eBoundPhi], fpar[Acts::eBoundTheta]); - Acts::Vector3 globalPosition = referenceSurface.localToGlobal( - gctx, fpar.segment<2>(Acts::eBoundLoc0), dir); - - // Rotation matrix. When applied to global coordinates, they - // are rotated into the local reference frame of the - // surface. Note that this such a rotation can be found by - // inverting a matrix whose columns correspond to the - // coordinate axes of the local coordinate system. - Acts::RotationMatrix3 rot = - referenceSurface.referenceFrame(gctx, globalPosition, dir) - .inverse(); - std::pair angles = - Acts::VectorHelpers::incidentAngles(dir, rot); - - input[iInput++] = angles.first; - input[iInput++] = angles.second; - input[iInput++] = fpar[Acts::eBoundLoc0]; - input[iInput++] = fpar[Acts::eBoundLoc1]; - input[iInput++] = fcov(Acts::eBoundLoc0, Acts::eBoundLoc0); - input[iInput++] = fcov(Acts::eBoundLoc1, Acts::eBoundLoc1); - if (iInput != m_nInputs) { - throw std::runtime_error("Expected input size of " + - std::to_string(m_nInputs) + - ", got: " + std::to_string(iInput)); - } - - // Input is a single row, hence .front() - std::vector output = - m_model.runONNXInference(inputBatch).front(); - // Assuming 2-D measurements, the expected params structure is: - // [ 0, nComponent[ --> priors - // [ nComponent, 3*nComponent[ --> means - // [3*nComponent, 5*nComponent[ --> variances - std::size_t nParams = 5 * m_nComponents; - if (output.size() != nParams) { - throw std::runtime_error( - "Got output vector of size " + std::to_string(output.size()) + - ", expected size " + std::to_string(nParams)); - } - - // Most probable value computation of mixture density - std::size_t iMax = 0; - if (m_nComponents > 1) { - iMax = std::distance( - output.begin(), - std::max_element(output.begin(), output.begin() + m_nComponents)); - } - std::size_t iLoc0 = m_nComponents + iMax * 2; - std::size_t iVar0 = 3 * m_nComponents + iMax * 2; - - fpar[Acts::eBoundLoc0] = output[iLoc0]; - fpar[Acts::eBoundLoc1] = output[iLoc0 + 1]; - fcov(Acts::eBoundLoc0, Acts::eBoundLoc0) = output[iVar0]; - fcov(Acts::eBoundLoc1, Acts::eBoundLoc1) = output[iVar0 + 1]; - - constexpr std::size_t kSize = - std::remove_reference_t::size(); - Acts::ActsVector cpar = P * fpar; - Acts::ActsSquareMatrix ccov = P * fcov * P.transpose(); - - trackState.allocateCalibrated(kSize); - trackState.template calibrated() = cpar; - trackState.template calibratedCovariance() = ccov; - trackState.setProjector(measurement.projector()); - }, - measurements[idxSourceLink.index()]); + // Most probable value computation of mixture density + std::size_t iMax = 0; + if (m_nComponents > 1) { + iMax = std::distance( + output.begin(), + std::max_element(output.begin(), output.begin() + m_nComponents)); + } + std::size_t iLoc0 = m_nComponents + iMax * 2; + std::size_t iVar0 = 3 * m_nComponents + iMax * 2; + + Measurement measurementCopy = measurement; + measurementCopy.effectiveParameters()[boundLoc0] = output[iLoc0]; + measurementCopy.effectiveParameters()[boundLoc1] = output[iLoc0 + 1]; + measurementCopy.effectiveCovariance()(boundLoc0, boundLoc0) = output[iVar0]; + measurementCopy.effectiveCovariance()(boundLoc1, boundLoc1) = + output[iVar0 + 1]; + + Acts::visit_measurement(measurement.size(), [&](auto N) -> void { + constexpr std::size_t kMeasurementSize = decltype(N)::value; + + trackState.allocateCalibrated(kMeasurementSize); + trackState.calibrated() = + measurementCopy.parameters(); + trackState.calibratedCovariance() = + measurementCopy.covariance(); + trackState.setProjector(measurementCopy.subspace().fullProjector()); + }); } diff --git a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp index 297508140f5..ec423531cfa 100644 --- a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/MeasurementHelpers.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/detail/CalculateResiduals.hpp" #include "Acts/EventData/detail/ParameterTraits.hpp" @@ -23,10 +24,9 @@ namespace ActsExamples { -/// A measurement of a fixed-size subspace of the full parameters. +/// A measurement of a variable-size subspace of the full parameters. /// /// @tparam indices_t Parameter index type, determines the full parameter space -/// @tparam kSize Size of the parameter subset. /// /// The measurement intentionally does not store a pointer/reference to the /// reference object in the geometry hierarchy, i.e. the surface or volume. The @@ -44,23 +44,46 @@ namespace ActsExamples { /// on). Both variants add additional complications. Since the geometry object /// is not required anyway (as discussed above), not storing it removes all /// these complications altogether. -template -class FixedSizeMeasurement { +template +class VariableSizeMeasurement { static constexpr std::size_t kFullSize = Acts::detail::kParametersSize; - using Subspace = Acts::detail::FixedSizeSubspace; + using Subspace = Acts::detail::VariableSizeSubspace; public: using Scalar = Acts::ActsScalar; + /// Vector type containing for measured parameter values. - using ParametersVector = Acts::ActsVector; + template + using ParametersVector = Eigen::Matrix; + template + using ParametersVectorMap = Eigen::Map>; + template + using ConstParametersVectorMap = Eigen::Map>; + using EffectiveParametersVector = Eigen::Matrix; + using EffectiveParametersVectorMap = Eigen::Map; + using ConstEffectiveParametersVectorMap = + Eigen::Map; + /// Matrix type for the measurement covariance. - using CovarianceMatrix = Acts::ActsSquareMatrix; - /// Vector type containing all parameters in the same space. + template + using CovarianceMatrix = Eigen::Matrix; + template + using CovarianceMatrixMap = Eigen::Map>; + template + using ConstCovarianceMatrixMap = Eigen::Map>; + using EffectiveCovarianceMatrix = + Eigen::Matrix; + using EffectiveCovarianceMatrixMap = Eigen::Map; + using ConstEffectiveCovarianceMatrixMap = + Eigen::Map; + using FullParametersVector = Acts::ActsVector; - using ProjectionMatrix = Acts::ActsMatrix; - using ExpansionMatrix = Acts::ActsMatrix; + using FullCovarianceMatrix = Acts::ActsSquareMatrix; + + using ProjectionMatrix = Eigen::Matrix; + using ExpansionMatrix = Eigen::Matrix; /// Construct from source link, subset indices, and measured data. /// @@ -73,100 +96,102 @@ class FixedSizeMeasurement { /// /// @note The indices must be ordered and must describe/match the content /// of parameters and covariance. - template - FixedSizeMeasurement(Acts::SourceLink source, - const std::array& indices, - const Eigen::MatrixBase& params, - const Eigen::MatrixBase& cov) - : m_source(std::move(source)), - m_subspace(indices), - m_params(params), - m_cov(cov) { - // TODO we should be able to support arbitrary ordering, by sorting the - // indices and reordering parameters/covariance. since the parameter order - // can be modified by the user, the user can not always know what the - // right order is. another option is too remove the requirement for index - // ordering from the subspace types, but that will make it harder to - // refactor their implementation later on. + template + VariableSizeMeasurement(Acts::SourceLink source, + const std::array& indices, + const Eigen::MatrixBase& params, + const Eigen::MatrixBase& cov) + : m_source(std::move(source)), m_subspace(indices) { + static_assert(kSize == parameters_t::RowsAtCompileTime, + "Parameter size mismatch"); + static_assert(kSize == covariance_t::RowsAtCompileTime, + "Covariance rows mismatch"); + static_assert(kSize == covariance_t::ColsAtCompileTime, + "Covariance cols mismatch"); + + parameters() = params; + covariance() = cov; } /// A measurement can only be constructed with valid parameters. - FixedSizeMeasurement() = delete; - FixedSizeMeasurement(const FixedSizeMeasurement&) = default; - FixedSizeMeasurement(FixedSizeMeasurement&&) = default; - ~FixedSizeMeasurement() = default; - FixedSizeMeasurement& operator=(const FixedSizeMeasurement&) = default; - FixedSizeMeasurement& operator=(FixedSizeMeasurement&&) = default; + VariableSizeMeasurement() = delete; + VariableSizeMeasurement(const VariableSizeMeasurement&) = default; + VariableSizeMeasurement(VariableSizeMeasurement&&) = default; + ~VariableSizeMeasurement() = default; + VariableSizeMeasurement& operator=(const VariableSizeMeasurement&) = default; + VariableSizeMeasurement& operator=(VariableSizeMeasurement&&) = default; + + constexpr std::size_t size() const { return m_subspace.size(); } /// Source link that connects to the underlying detector readout. const Acts::SourceLink& sourceLink() const { return m_source; } - /// Number of measured parameters. - static constexpr std::size_t size() { return kSize; } + const Subspace& subspace() const { return m_subspace; } /// Check if a specific parameter is part of this measurement. bool contains(indices_t i) const { return m_subspace.contains(i); } - /// The measurement indices - constexpr std::array indices() const { - std::array subInds = m_subspace.indices(); - std::array inds{}; - for (std::size_t i = 0; i < kSize; i++) { - inds[i] = static_cast(subInds[i]); - } - return inds; + template + ConstParametersVectorMap parameters() const { + return ConstParametersVectorMap{m_params.data()}; } - - /// Measured parameters values. - const ParametersVector& parameters() const { return m_params; } - - /// Measured parameters covariance. - const CovarianceMatrix& covariance() const { return m_cov; } - - /// Projection matrix from the full space into the measured subspace. - ProjectionMatrix projector() const { - return m_subspace.template projector(); + template + ParametersVectorMap parameters() { + return ParametersVectorMap{m_params.data()}; + } + ConstEffectiveParametersVectorMap effectiveParameters() const { + return ConstEffectiveParametersVectorMap{m_params.data(), + static_cast(size())}; + } + EffectiveParametersVectorMap effectiveParameters() { + return EffectiveParametersVectorMap{m_params.data(), + static_cast(size())}; } - /// Expansion matrix from the measured subspace into the full space. - /// - /// This is equivalent to the transpose of the projection matrix only in the - /// case of a trivial projection matrix. While this is the case here, it is - /// still recommended to use the expansion matrix directly in cases where it - /// is explicitly used. - ExpansionMatrix expander() const { - return m_subspace.template expander(); + template + ConstCovarianceMatrixMap covariance() const { + return ConstCovarianceMatrixMap{m_cov.data()}; + } + template + CovarianceMatrixMap covariance() { + return CovarianceMatrixMap{m_cov.data()}; + } + ConstEffectiveCovarianceMatrixMap effectiveCovariance() const { + return ConstEffectiveCovarianceMatrixMap{m_cov.data(), + static_cast(size()), + static_cast(size())}; + } + EffectiveCovarianceMatrixMap effectiveCovariance() { + return EffectiveCovarianceMatrixMap{m_cov.data(), + static_cast(size()), + static_cast(size())}; } - /// Compute residuals in the measured subspace. - /// - /// @param reference Reference parameters in the full space. - /// - /// This computes the difference `measured - reference` taking into account - /// the allowed parameter ranges. Only the reference values in the measured - /// subspace are used for the computation. - ParametersVector residuals(const FullParametersVector& reference) const { - ParametersVector res = ParametersVector::Zero(); - Acts::detail::calculateResiduals(static_cast(kSize), - m_subspace.indices(), reference, m_params, - res); - return res; + FullParametersVector fullParameters() const { + FullParametersVector result = FullParametersVector::Zero(); + for (std::size_t i = 0; i < size(); ++i) { + result[m_subspace[i]] = effectiveParameters()[i]; + } + return result; } - std::ostream& operator<<(std::ostream& os) const { - Acts::detail::printMeasurement(os, static_cast(kSize), - m_subspace.indices().data(), m_params.data(), - m_cov.data()); - return os; + FullCovarianceMatrix fullCovariance() const { + FullCovarianceMatrix result = FullCovarianceMatrix::Zero(); + for (std::size_t i = 0; i < size(); ++i) { + for (std::size_t j = 0; j < size(); ++j) { + result(m_subspace[i], m_subspace[j]) = effectiveCovariance()(i, j); + } + } + return result; } private: Acts::SourceLink m_source; Subspace m_subspace; - ParametersVector m_params; - CovarianceMatrix m_cov; + std::array m_params{}; + std::array m_cov{}; }; -/// Construct a fixed-size measurement for the given indices. +/// Construct a variable-size measurement for the given indices. /// /// @tparam parameters_t Input parameters vector type /// @tparam covariance_t Input covariance matrix type @@ -178,77 +203,26 @@ class FixedSizeMeasurement { /// @param cov Measured parameters covariance /// @param index0 Required parameter index, measurement must be at least 1d /// @param tailIndices Additional parameter indices for larger measurements -/// @return Fixed-size measurement w/ the correct type and given inputs -/// -/// This helper function can be used to create a fixed-size measurement using an -/// explicit set of indices, e.g. -/// -/// auto m = makeMeasurement(s, p, c, eBoundLoc0, eBoundTime); -/// -/// for a 2d measurement w/ one position and time. +/// @return Variable-size measurement w/ the correct type and given inputs /// /// @note The indices must be ordered and must be consistent with the content of -/// parameters and covariance. +/// parameters and covariance. template -auto makeFixedSizeMeasurement(Acts::SourceLink source, - const Eigen::MatrixBase& params, - const Eigen::MatrixBase& cov, - indices_t index0, tail_indices_t... tailIndices) - -> FixedSizeMeasurement { +VariableSizeMeasurement makeVariableSizeMeasurement( + Acts::SourceLink source, const Eigen::MatrixBase& params, + const Eigen::MatrixBase& cov, indices_t index0, + tail_indices_t... tailIndices) { using IndexContainer = std::array; return {std::move(source), IndexContainer{index0, tailIndices...}, params, cov}; } -namespace detail { -/// @cond - -// Recursive construction of the measurement variant. `kN` is counted down until -// zero while the sizes are accumulated in the parameter pack. -// -// Example: -// -// VariantMeasurementGenerator<..., 4> -// -> VariantMeasurementGenerator<..., 3, 4> -// -> VariantMeasurementGenerator<..., 2, 3, 4> -// -> VariantMeasurementGenerator<..., 1, 2, 3, 4> -// -> VariantMeasurementGenerator<..., 0, 1, 2, 3, 4> -// -template -struct VariantMeasurementGenerator - : VariantMeasurementGenerator {}; -template -struct VariantMeasurementGenerator { - using Type = std::variant...>; -}; - -/// @endcond -} // namespace detail - -/// Variant that can contain all possible measurements in a parameter space. -/// -/// @tparam indices_t Parameter index type, determines the full parameter space -template -using VariantMeasurement = typename detail::VariantMeasurementGenerator< - indices_t, Acts::detail::kParametersSize>::Type; - -/// Variant that can hold all possible bound measurements. -/// -using BoundVariantMeasurement = VariantMeasurement; - -/// Variant that can hold all possible free measurements. -/// -using FreeVariantMeasurement = VariantMeasurement; - -template -std::ostream& operator<<(std::ostream& os, - const VariantMeasurement& vm) { - return std::visit([&](const auto& m) { return (os << m); }, vm); -} +/// Type that can hold all possible bound measurements. +using BoundVariableMeasurement = VariableSizeMeasurement; /// Variable measurement type that can contain all possible combinations. -using Measurement = BoundVariantMeasurement; +using Measurement = BoundVariableMeasurement; /// Container of measurements. /// diff --git a/Examples/Framework/src/EventData/MeasurementCalibration.cpp b/Examples/Framework/src/EventData/MeasurementCalibration.cpp index dd0f45bee68..444021c91ca 100644 --- a/Examples/Framework/src/EventData/MeasurementCalibration.cpp +++ b/Examples/Framework/src/EventData/MeasurementCalibration.cpp @@ -6,6 +6,7 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/SourceLink.hpp" #include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/Measurement.hpp" @@ -30,19 +31,18 @@ void ActsExamples::PassThroughCalibrator::calibrate( assert((idxSourceLink.index() < measurements.size()) && "Source link index is outside the container bounds"); - std::visit( - [&trackState](const auto& meas) { - using MeasurementType = std::decay_t; - constexpr std::size_t size = MeasurementType::size(); + const auto& measurement = measurements[idxSourceLink.index()]; - trackState.allocateCalibrated(meas.size()); - assert(trackState.hasCalibrated()); + Acts::visit_measurement(measurement.size(), [&](auto N) -> void { + constexpr std::size_t kMeasurementSize = decltype(N)::value; - trackState.calibrated() = meas.parameters(); - trackState.calibratedCovariance() = meas.covariance(); - trackState.setProjector(meas.projector()); - }, - measurements[idxSourceLink.index()]); + trackState.allocateCalibrated(kMeasurementSize); + trackState.calibrated() = + measurement.parameters(); + trackState.calibratedCovariance() = + measurement.covariance(); + trackState.setProjector(measurement.subspace().fullProjector()); + }); } ActsExamples::MeasurementCalibratorAdapter::MeasurementCalibratorAdapter( diff --git a/Examples/Framework/src/EventData/ScalingCalibrator.cpp b/Examples/Framework/src/EventData/ScalingCalibrator.cpp index 5172a9eda47..bac4f1d8cae 100644 --- a/Examples/Framework/src/EventData/ScalingCalibrator.cpp +++ b/Examples/Framework/src/EventData/ScalingCalibrator.cpp @@ -151,34 +151,30 @@ void ActsExamples::ScalingCalibrator::calibrate( const Cluster& cl = clusters->at(idxSourceLink.index()); ConstantTuple ct = m_calib_maps.at(mgid).at(cl.sizeLoc0, cl.sizeLoc1); - std::visit( - [&](const auto& meas) { - auto E = meas.expander(); - auto P = meas.projector(); - - Acts::ActsVector fpar = E * meas.parameters(); - - Acts::ActsSquareMatrix fcov = - E * meas.covariance() * E.transpose(); - - fpar[Acts::eBoundLoc0] += ct.x_offset; - fpar[Acts::eBoundLoc1] += ct.y_offset; - fcov(Acts::eBoundLoc0, Acts::eBoundLoc0) *= ct.x_scale; - fcov(Acts::eBoundLoc1, Acts::eBoundLoc1) *= ct.y_scale; - - constexpr std::size_t kSize = - std::remove_reference_t::size(); - std::array indices = meas.indices(); - Acts::ActsVector cpar = P * fpar; - Acts::ActsSquareMatrix ccov = P * fcov * P.transpose(); - - FixedSizeMeasurement cmeas( - Acts::SourceLink{idxSourceLink}, indices, cpar, ccov); - - trackState.allocateCalibrated(cmeas.size()); - trackState.calibrated() = meas.parameters(); - trackState.calibratedCovariance() = meas.covariance(); - trackState.setProjector(meas.projector()); - }, - (measurements)[idxSourceLink.index()]); + const auto& measurement = measurements[idxSourceLink.index()]; + + assert(measurement.contains(Acts::eBoundLoc0) && + "Measurement does not contain the required bound loc0"); + assert(measurement.contains(Acts::eBoundLoc1) && + "Measurement does not contain the required bound loc1"); + + auto boundLoc0 = measurement.subspace().indexOf(Acts::eBoundLoc0); + auto boundLoc1 = measurement.subspace().indexOf(Acts::eBoundLoc1); + + Measurement measurementCopy = measurement; + measurementCopy.effectiveParameters()[boundLoc0] += ct.x_offset; + measurementCopy.effectiveParameters()[boundLoc1] += ct.y_offset; + measurementCopy.effectiveCovariance()(boundLoc0, boundLoc0) *= ct.x_scale; + measurementCopy.effectiveCovariance()(boundLoc1, boundLoc1) *= ct.y_scale; + + Acts::visit_measurement(measurement.size(), [&](auto N) -> void { + constexpr std::size_t kMeasurementSize = decltype(N)::value; + + trackState.allocateCalibrated(kMeasurementSize); + trackState.calibrated() = + measurement.parameters(); + trackState.calibratedCovariance() = + measurementCopy.covariance(); + trackState.setProjector(measurement.subspace().fullProjector()); + }); } diff --git a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp index a5003650124..0f65fca0d20 100644 --- a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp +++ b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp @@ -99,57 +99,52 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementWriter::writeT( writerMeasurementSimHitMap.append({measIdx, simHitIdx}); } - std::visit( - [&](const auto& m) { - Acts::GeometryIdentifier geoId = - m.sourceLink().template get().geometryId(); - // MEASUREMENT information ------------------------------------ - - // Encoded geometry identifier. same for all hits on the module - meas.geometry_id = geoId.value(); - meas.local_key = 0; - // Create a full set of parameters - auto parameters = (m.expander() * m.parameters()).eval(); - meas.local0 = parameters[Acts::eBoundLoc0] / Acts::UnitConstants::mm; - meas.local1 = parameters[Acts::eBoundLoc1] / Acts::UnitConstants::mm; - meas.phi = parameters[Acts::eBoundPhi] / Acts::UnitConstants::rad; - meas.theta = parameters[Acts::eBoundTheta] / Acts::UnitConstants::rad; - meas.time = parameters[Acts::eBoundTime] / Acts::UnitConstants::mm; - - auto covariance = - (m.expander() * m.covariance() * m.expander().transpose()).eval(); - meas.var_local0 = covariance(Acts::eBoundLoc0, Acts::eBoundLoc0); - meas.var_local1 = covariance(Acts::eBoundLoc1, Acts::eBoundLoc1); - meas.var_phi = covariance(Acts::eBoundPhi, Acts::eBoundPhi); - meas.var_theta = covariance(Acts::eBoundTheta, Acts::eBoundTheta); - meas.var_time = covariance(Acts::eBoundTime, Acts::eBoundTime); - for (unsigned int ipar = 0; - ipar < static_cast(Acts::eBoundSize); ++ipar) { - if (m.contains(static_cast(ipar))) { - meas.local_key = ((1 << (ipar + 1)) | meas.local_key); - } - } - - writerMeasurements.append(meas); - - // CLUSTER / channel information ------------------------------ - if (!clusters.empty() && writerCells) { - auto cluster = clusters[measIdx]; - cell.geometry_id = meas.geometry_id; - cell.measurement_id = meas.measurement_id; - for (auto& c : cluster.channels) { - cell.channel0 = c.bin[0]; - cell.channel1 = c.bin[1]; - // TODO store digital timestamp once added to the cell definition - cell.timestamp = 0; - cell.value = c.activation; - writerCells->append(cell); - } - } - // Increase counter - meas.measurement_id += 1; - }, - measurement); + Acts::GeometryIdentifier geoId = + measurement.sourceLink().template get().geometryId(); + // MEASUREMENT information ------------------------------------ + + // Encoded geometry identifier. same for all hits on the module + meas.geometry_id = geoId.value(); + meas.local_key = 0; + // Create a full set of parameters + auto parameters = measurement.fullParameters(); + meas.local0 = parameters[Acts::eBoundLoc0] / Acts::UnitConstants::mm; + meas.local1 = parameters[Acts::eBoundLoc1] / Acts::UnitConstants::mm; + meas.phi = parameters[Acts::eBoundPhi] / Acts::UnitConstants::rad; + meas.theta = parameters[Acts::eBoundTheta] / Acts::UnitConstants::rad; + meas.time = parameters[Acts::eBoundTime] / Acts::UnitConstants::mm; + + auto covariance = measurement.fullCovariance(); + meas.var_local0 = covariance(Acts::eBoundLoc0, Acts::eBoundLoc0); + meas.var_local1 = covariance(Acts::eBoundLoc1, Acts::eBoundLoc1); + meas.var_phi = covariance(Acts::eBoundPhi, Acts::eBoundPhi); + meas.var_theta = covariance(Acts::eBoundTheta, Acts::eBoundTheta); + meas.var_time = covariance(Acts::eBoundTime, Acts::eBoundTime); + for (unsigned int ipar = 0; + ipar < static_cast(Acts::eBoundSize); ++ipar) { + if (measurement.contains(static_cast(ipar))) { + meas.local_key = ((1 << (ipar + 1)) | meas.local_key); + } + } + + writerMeasurements.append(meas); + + // CLUSTER / channel information ------------------------------ + if (!clusters.empty() && writerCells) { + auto cluster = clusters[measIdx]; + cell.geometry_id = meas.geometry_id; + cell.measurement_id = meas.measurement_id; + for (auto& c : cluster.channels) { + cell.channel0 = c.bin[0]; + cell.channel1 = c.bin[1]; + // TODO store digital timestamp once added to the cell definition + cell.timestamp = 0; + cell.value = c.activation; + writerCells->append(cell); + } + } + // Increase counter + meas.measurement_id += 1; } return ActsExamples::ProcessCode::SUCCESS; } diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp index 056b616d72b..6a7adfe98c5 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp @@ -201,54 +201,48 @@ void EDM4hepUtil::writeMeasurement(const Measurement& from, const Cluster* fromCluster, edm4hep::TrackerHitCollection& toClusters, const MapGeometryIdTo& geometryMapper) { - std::visit( - [&](const auto& m) { - Acts::GeometryIdentifier geoId = - m.sourceLink().template get().geometryId(); - - if (geometryMapper) { - // no need for digitization as we only want to identify the sensor - to.setCellID(geometryMapper(geoId)); - } - - auto parameters = (m.expander() * m.parameters()).eval(); - - to.setTime(parameters[Acts::eBoundTime] / Acts::UnitConstants::ns); - - to.setType(Acts::EDM4hepUtil::EDM4HEP_ACTS_POSITION_TYPE); - // TODO set uv (which are in global spherical coordinates with r=1) - to.setPosition({parameters[Acts::eBoundLoc0], - parameters[Acts::eBoundLoc1], - parameters[Acts::eBoundTime]}); - - auto covariance = - (m.expander() * m.covariance() * m.expander().transpose()).eval(); - to.setCovMatrix({ - static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), - static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc0)), - static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), - 0, - 0, - 0, - }); - - if (fromCluster) { - 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); - } - } - }, - from); + Acts::GeometryIdentifier geoId = + from.sourceLink().template get().geometryId(); + + if (geometryMapper) { + // no need for digitization as we only want to identify the sensor + to.setCellID(geometryMapper(geoId)); + } + + const auto& parameters = from.fullParameters(); + const auto& covariance = from.fullCovariance(); + + to.setTime(parameters[Acts::eBoundTime] / Acts::UnitConstants::ns); + + to.setType(Acts::EDM4hepUtil::EDM4HEP_ACTS_POSITION_TYPE); + // TODO set uv (which are in global spherical coordinates with r=1) + to.setPosition({parameters[Acts::eBoundLoc0], parameters[Acts::eBoundLoc1], + parameters[Acts::eBoundTime]}); + + to.setCovMatrix({ + static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), + static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc0)), + static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), + 0, + 0, + 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); + } + } } void EDM4hepUtil::writeTrajectory( diff --git a/Examples/Io/Root/src/RootMeasurementWriter.cpp b/Examples/Io/Root/src/RootMeasurementWriter.cpp index 9cc695596e4..0a476b031a1 100644 --- a/Examples/Io/Root/src/RootMeasurementWriter.cpp +++ b/Examples/Io/Root/src/RootMeasurementWriter.cpp @@ -151,14 +151,13 @@ struct RootMeasurementWriter::DigitizationTree { /// Convenience function to fill bound parameters /// - /// @tparam measurement_t Type of the parameter set - /// - /// @param m The measurement set - template - void fillBoundMeasurement(const measurement_t& m) { - for (auto [i, ib] : Acts::enumerate(m.indices())) { - recBound[ib] = m.parameters()[i]; - varBound[ib] = m.covariance()(i, i); + /// @param m The measurement + void fillBoundMeasurement(const Measurement& m) { + for (unsigned int i = 0; i < m.size(); ++i) { + auto ib = m.subspace()[i]; + + recBound[ib] = m.effectiveParameters()[i]; + varBound[ib] = m.effectiveCovariance()(i, i); residual[ib] = recBound[ib] - trueBound[ib]; pull[ib] = residual[ib] / std::sqrt(varBound[ib]); @@ -269,44 +268,39 @@ ProcessCode RootMeasurementWriter::writeT( for (Index hitIdx = 0u; hitIdx < measurements.size(); ++hitIdx) { const auto& meas = measurements[hitIdx]; - std::visit( - [&](const auto& m) { - Acts::GeometryIdentifier geoId = - m.sourceLink().template get().geometryId(); - // find the corresponding surface - auto surfaceItr = m_cfg.surfaceByIdentifier.find(geoId); - if (surfaceItr == m_cfg.surfaceByIdentifier.end()) { - return; - } - const Acts::Surface& surface = *(surfaceItr->second); - - // Fill the identification - m_outputTree->fillIdentification(ctx.eventNumber, geoId); - - // Find the contributing simulated hits - auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx)); - // Use average truth in the case of multiple contributing sim hits - auto [local, pos4, dir] = averageSimHits(ctx.geoContext, surface, - simHits, indices, logger()); - Acts::RotationMatrix3 rot = - surface - .referenceFrame(ctx.geoContext, pos4.segment<3>(Acts::ePos0), - dir) - .inverse(); - std::pair angles = - Acts::VectorHelpers::incidentAngles(dir, rot); - - m_outputTree->fillTruthParameters(local, pos4, dir, angles); - m_outputTree->fillBoundMeasurement(m); - if (clusters != nullptr) { - const auto& c = (*clusters)[hitIdx]; - m_outputTree->fillCluster(c); - } - - m_outputTree->fill(); - m_outputTree->clear(); - }, - meas); + Acts::GeometryIdentifier geoId = + meas.sourceLink().template get().geometryId(); + // find the corresponding surface + auto surfaceItr = m_cfg.surfaceByIdentifier.find(geoId); + if (surfaceItr == m_cfg.surfaceByIdentifier.end()) { + continue; + } + const Acts::Surface& surface = *(surfaceItr->second); + + // Fill the identification + m_outputTree->fillIdentification(ctx.eventNumber, geoId); + + // Find the contributing simulated hits + auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx)); + // Use average truth in the case of multiple contributing sim hits + auto [local, pos4, dir] = + averageSimHits(ctx.geoContext, surface, simHits, indices, logger()); + Acts::RotationMatrix3 rot = + surface + .referenceFrame(ctx.geoContext, pos4.segment<3>(Acts::ePos0), dir) + .inverse(); + std::pair angles = + Acts::VectorHelpers::incidentAngles(dir, rot); + + m_outputTree->fillTruthParameters(local, pos4, dir, angles); + m_outputTree->fillBoundMeasurement(meas); + if (clusters != nullptr) { + const auto& c = (*clusters)[hitIdx]; + m_outputTree->fillCluster(c); + } + + m_outputTree->fill(); + m_outputTree->clear(); } return ProcessCode::SUCCESS; diff --git a/Tests/UnitTests/Core/Utilities/SubspaceTests.cpp b/Tests/UnitTests/Core/Utilities/SubspaceTests.cpp index 0d91065089e..2fbf0d2ffc4 100644 --- a/Tests/UnitTests/Core/Utilities/SubspaceTests.cpp +++ b/Tests/UnitTests/Core/Utilities/SubspaceTests.cpp @@ -164,8 +164,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(VariableSizeSubspace, ScalarAndSubspace, BOOST_CHECK_EQUAL(variableSubspace.fullSize(), fixedSubspace.fullSize()); auto fixedProjector = fixedSubspace.template projector(); - std::uint64_t fixedProjectorBits = - matrixToBitset(fixedProjector).to_ullong(); Eigen::Matrix fixedFullProjector; @@ -173,15 +171,11 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(VariableSizeSubspace, ScalarAndSubspace, fixedFullProjector.template topLeftCorner() = fixedProjector; - std::uint64_t fixedFullProjectorBits = - matrixToBitset(fixedFullProjector).to_ullong(); - std::uint64_t variableProjectorBits = variableSubspace.projectorBits(); - std::uint64_t variableFullProjectorBits = - variableSubspace.fullProjectorBits(); + auto variableFullProjector = + variableSubspace.template fullProjector(); - BOOST_CHECK_EQUAL(variableProjectorBits, fixedProjectorBits); - BOOST_CHECK_EQUAL(variableFullProjectorBits, fixedFullProjectorBits); + BOOST_CHECK_EQUAL(variableFullProjector, fixedFullProjector); } while (std::next_permutation(fullIndices.begin(), fullIndices.end())); } diff --git a/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp b/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp index 7f2a6d1e711..cd9e8904a8b 100644 --- a/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp +++ b/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2016-2020 CERN for the benefit of the Acts project +// Copyright (C) 2016-2024 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 @@ -37,10 +37,6 @@ namespace { constexpr BoundIndices boundIndices[] = { eBoundLoc0, eBoundLoc1, eBoundTime, eBoundPhi, eBoundTheta, eBoundQOverP, }; -constexpr FreeIndices freeIndices[] = { - eFreePos0, eFreePos1, eFreePos2, eFreeTime, - eFreeDir0, eFreeDir1, eFreeDir2, eFreeQOverP, -}; const TestSourceLink sourceOrig; const Acts::SourceLink source{sourceOrig}; // fix seed for reproducible tests @@ -53,9 +49,9 @@ std::default_random_engine rng(123); BOOST_AUTO_TEST_SUITE(EventDataMeasurement) -BOOST_DATA_TEST_CASE(FixedBoundOne, bd::make(boundIndices), index) { +BOOST_DATA_TEST_CASE(VariableBoundOne, bd::make(boundIndices), index) { auto [params, cov] = generateParametersCovariance(rng); - auto meas = makeFixedSizeMeasurement(source, params, cov, index); + auto meas = makeVariableSizeMeasurement(source, params, cov, index); BOOST_CHECK_EQUAL(meas.size(), 1); for (auto i : boundIndices) { @@ -65,184 +61,51 @@ BOOST_DATA_TEST_CASE(FixedBoundOne, bd::make(boundIndices), index) { BOOST_CHECK(!meas.contains(i)); } } - BOOST_CHECK_EQUAL(meas.parameters(), params); - BOOST_CHECK_EQUAL(meas.covariance(), cov); + BOOST_CHECK_EQUAL(meas.effectiveParameters(), params); + BOOST_CHECK_EQUAL(meas.effectiveCovariance(), cov); BOOST_CHECK_EQUAL(meas.sourceLink().template get(), sourceOrig); } -BOOST_AUTO_TEST_CASE(FixedBoundAll) { +BOOST_AUTO_TEST_CASE(VariableBoundAll) { auto [params, cov] = generateBoundParametersCovariance(rng); - auto meas = makeFixedSizeMeasurement(source, params, cov, eBoundLoc0, - eBoundLoc1, eBoundPhi, eBoundTheta, - eBoundQOverP, eBoundTime); + auto meas = makeVariableSizeMeasurement(source, params, cov, eBoundLoc0, + eBoundLoc1, eBoundPhi, eBoundTheta, + eBoundQOverP, eBoundTime); BOOST_CHECK_EQUAL(meas.size(), eBoundSize); for (auto i : boundIndices) { BOOST_CHECK(meas.contains(i)); } - BOOST_CHECK_EQUAL(meas.parameters(), params); - BOOST_CHECK_EQUAL(meas.covariance(), cov); - BOOST_CHECK_EQUAL(meas.sourceLink().get(), sourceOrig); -} - -namespace { -// example data for phi residual tests. each entry contains -// -// measured, reference, expected residual -// -const std::vector> kPhiDataset = { - // measurement and reference in bounds and close - {0.5, 0.75, -0.25}, - // measurement and reference in bounds but at different edges - {0.25, 2 * M_PI - 0.25, 0.5}, - {2 * M_PI - 0.125, 0.125, -0.25}, - // measurement in bounds, reference ouf-of-bounds, both near lower edge - {0.25, -0.25, 0.5}, - // measurement in bounds, reference ouf-of-bounds, both near upper edge - {2 * M_PI - 0.25, 2 * M_PI + 0.25, -0.5}, - // measurement out-of-bounds, reference in bounds, both near lower edge - {-0.25, 0.25, -0.5}, - // measurement out-of-bounds, reference in bounds, both near upper edge - {2 * M_PI + 0.25, 2 * M_PI - 0.25, 0.5}, -}; -} // namespace - -BOOST_DATA_TEST_CASE(BoundResidualsPhi, bd::make(kPhiDataset), phiMea, phiRef, - phiRes) { - using MeasurementVector = Acts::ActsVector<1>; - using MeasurementCovariance = Acts::ActsSquareMatrix<1>; - - // prepare measurement - MeasurementVector params = MeasurementVector::Zero(); - MeasurementCovariance cov = MeasurementCovariance::Zero(); - params[0] = phiMea; - auto measurement = makeFixedSizeMeasurement(source, params, cov, eBoundPhi); - // prepare reference parameters - Acts::BoundVector reference = Acts::BoundVector::Zero(); - reference[eBoundPhi] = phiRef; - - // compute and check residual - auto res = measurement.residuals(reference); - CHECK_CLOSE_ABS(res[0], phiRes, std::numeric_limits::epsilon()); -} - -BOOST_DATA_TEST_CASE(FixedFreeOne, bd::make(freeIndices), index) { - auto [params, cov] = generateParametersCovariance(rng); - auto meas = makeFixedSizeMeasurement(source, params, cov, index); - - BOOST_CHECK_EQUAL(meas.size(), 1); - for (auto i : freeIndices) { - if (i == index) { - BOOST_CHECK(meas.contains(i)); - } else { - BOOST_CHECK(!meas.contains(i)); - } - } - BOOST_CHECK_EQUAL(meas.parameters(), params); - BOOST_CHECK_EQUAL(meas.covariance(), cov); - BOOST_CHECK_EQUAL(meas.sourceLink().template get(), - sourceOrig); - - // all free parameters are unrestricted and we know the expected residual. - constexpr auto tol = std::numeric_limits::epsilon(); - auto [ref, refCov] = generateFreeParametersCovariance(rng); - auto res = meas.residuals(ref); - CHECK_CLOSE_ABS(res[0], params[0] - ref[index], tol); -} - -BOOST_AUTO_TEST_CASE(FixedFreeAll) { - auto [params, cov] = generateFreeParametersCovariance(rng); - auto meas = makeFixedSizeMeasurement( - source, params, cov, eFreePos0, eFreePos1, eFreePos2, eFreeTime, - eFreeDir0, eFreeDir1, eFreeDir2, eFreeQOverP); - - BOOST_CHECK_EQUAL(meas.size(), eFreeSize); - for (auto i : freeIndices) { - BOOST_CHECK(meas.contains(i)); - } - BOOST_CHECK_EQUAL(meas.parameters(), params); - BOOST_CHECK_EQUAL(meas.covariance(), cov); + BOOST_CHECK_EQUAL(meas.effectiveParameters(), params); + BOOST_CHECK_EQUAL(meas.effectiveCovariance(), cov); BOOST_CHECK_EQUAL(meas.sourceLink().get(), sourceOrig); - - // all free parameters are unrestricted and we know the expected residual. - constexpr auto tol = std::numeric_limits::epsilon(); - auto [ref, refCov] = generateFreeParametersCovariance(rng); - CHECK_CLOSE_ABS(meas.residuals(ref), params - ref, tol); } -BOOST_AUTO_TEST_CASE(VariantBound) { +BOOST_AUTO_TEST_CASE(VariableBoundReassign) { // generate w/ a single parameter auto [par1, cov1] = generateParametersCovariance(rng); - BoundVariantMeasurement meas = - makeFixedSizeMeasurement(source, par1, cov1, eBoundTheta); - std::visit( - [](const auto& m) { - BOOST_CHECK_EQUAL(m.size(), 1); - BOOST_CHECK(!m.contains(eBoundLoc0)); - BOOST_CHECK(!m.contains(eBoundLoc1)); - BOOST_CHECK(!m.contains(eBoundTime)); - BOOST_CHECK(!m.contains(eBoundPhi)); - BOOST_CHECK(m.contains(eBoundTheta)); - BOOST_CHECK(!m.contains(eBoundQOverP)); - }, - meas); + auto meas = makeVariableSizeMeasurement(source, par1, cov1, eBoundTheta); + BOOST_CHECK_EQUAL(meas.size(), 1); + BOOST_CHECK(!meas.contains(eBoundLoc0)); + BOOST_CHECK(!meas.contains(eBoundLoc1)); + BOOST_CHECK(!meas.contains(eBoundTime)); + BOOST_CHECK(!meas.contains(eBoundPhi)); + BOOST_CHECK(meas.contains(eBoundTheta)); + BOOST_CHECK(!meas.contains(eBoundQOverP)); // reassign w/ all parameters auto [parN, covN] = generateBoundParametersCovariance(rng); - meas = makeFixedSizeMeasurement(source, parN, covN, eBoundLoc0, eBoundLoc1, - eBoundPhi, eBoundTheta, eBoundQOverP, - eBoundTime); - std::visit( - [](const auto& m) { - BOOST_CHECK_EQUAL(m.size(), eBoundSize); - BOOST_CHECK(m.contains(eBoundLoc0)); - BOOST_CHECK(m.contains(eBoundLoc1)); - BOOST_CHECK(m.contains(eBoundTime)); - BOOST_CHECK(m.contains(eBoundPhi)); - BOOST_CHECK(m.contains(eBoundTheta)); - BOOST_CHECK(m.contains(eBoundQOverP)); - }, - meas); -} - -BOOST_AUTO_TEST_CASE(VariantFree) { - // generate w/ two parameters - auto [par2, cov2] = generateParametersCovariance(rng); - FreeVariantMeasurement meas = - makeFixedSizeMeasurement(source, par2, cov2, eFreePos2, eFreeTime); - std::visit( - [](const auto& m) { - BOOST_CHECK_EQUAL(m.size(), 2); - BOOST_CHECK(!m.contains(eFreePos0)); - BOOST_CHECK(!m.contains(eFreePos1)); - BOOST_CHECK(m.contains(eFreePos2)); - BOOST_CHECK(m.contains(eFreeTime)); - BOOST_CHECK(!m.contains(eFreeDir0)); - BOOST_CHECK(!m.contains(eFreeDir1)); - BOOST_CHECK(!m.contains(eFreeDir2)); - BOOST_CHECK(!m.contains(eFreeQOverP)); - }, - meas); - - // reassign w/ all parameters - auto [parN, covN] = generateFreeParametersCovariance(rng); - meas = makeFixedSizeMeasurement(source, parN, covN, eFreePos0, eFreePos1, - eFreePos2, eFreeTime, eFreeDir0, eFreeDir1, - eFreeDir2, eFreeQOverP); - std::visit( - [](const auto& m) { - BOOST_CHECK_EQUAL(m.size(), eFreeSize); - BOOST_CHECK(m.contains(eFreePos0)); - BOOST_CHECK(m.contains(eFreePos1)); - BOOST_CHECK(m.contains(eFreePos2)); - BOOST_CHECK(m.contains(eFreeTime)); - BOOST_CHECK(m.contains(eFreeDir0)); - BOOST_CHECK(m.contains(eFreeDir1)); - BOOST_CHECK(m.contains(eFreeDir2)); - BOOST_CHECK(m.contains(eFreeQOverP)); - }, - meas); + meas = makeVariableSizeMeasurement(source, parN, covN, eBoundLoc0, eBoundLoc1, + eBoundPhi, eBoundTheta, eBoundQOverP, + eBoundTime); + BOOST_CHECK_EQUAL(meas.size(), eBoundSize); + BOOST_CHECK(meas.contains(eBoundLoc0)); + BOOST_CHECK(meas.contains(eBoundLoc1)); + BOOST_CHECK(meas.contains(eBoundTime)); + BOOST_CHECK(meas.contains(eBoundPhi)); + BOOST_CHECK(meas.contains(eBoundTheta)); + BOOST_CHECK(meas.contains(eBoundQOverP)); } BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp index 9c896ee6a93..c9e9fa4126e 100644 --- a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp +++ b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2023 CERN for the benefit of the Acts project +// Copyright (C) 2023-2024 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 @@ -50,11 +50,9 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) { Acts::Vector2 p = Acts::Vector2::Random(); Acts::SquareMatrix2 c = Acts::SquareMatrix2::Random(); - // NOTE this fails: - // auto m = Acts::makeMeasurement(sl, p, c, eBoundLoc0, eBoundTime) - // because we don't support non-consecutive parameters here for now - auto m = makeFixedSizeMeasurement(Acts::SourceLink{sl}, p, c, - Acts::eBoundLoc0, Acts::eBoundLoc1); + BoundVariableMeasurement m(Acts::SourceLink{sl}, + std::array{Acts::eBoundLoc0, Acts::eBoundLoc1}, + p, c); measOriginal.push_back(m); @@ -132,19 +130,13 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) { /////////// // Check // /////////// - auto checkMeasurementClose = [](const auto &m1, const auto &m2) { - constexpr auto SizeA = std::decay_t::size(); - constexpr auto SizeB = std::decay_t::size(); - if constexpr (SizeA == SizeB) { - CHECK_CLOSE_REL(m1.parameters(), m2.parameters(), 1e-4); - } - }; - static_assert( std::is_same_v, decltype(measOriginal)>); BOOST_REQUIRE(measRead.size() == measOriginal.size()); for (const auto &[a, b] : Acts::zip(measRead, measOriginal)) { - std::visit(checkMeasurementClose, a, b); + if (a.size() == b.size()) { + CHECK_CLOSE_REL(a.effectiveParameters(), b.effectiveParameters(), 1e-4); + } } static_assert(std::is_same_v,