diff --git a/Alignment/include/ActsAlignment/Kernel/Alignment.hpp b/Alignment/include/ActsAlignment/Kernel/Alignment.hpp index 1ce2eb3b959..83c7b569315 100644 --- a/Alignment/include/ActsAlignment/Kernel/Alignment.hpp +++ b/Alignment/include/ActsAlignment/Kernel/Alignment.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/Alignment.hpp" +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" @@ -189,6 +190,17 @@ struct Alignment { const fit_options_t& fitOptions, AlignmentResult& alignResult, const AlignmentMask& alignMask = AlignmentMask::All) const; + /// @brief calculate the alignment parameters delta from a set of + /// TrackAlignmentStates + /// + /// @param TrackStateCollection The collection of TrackAlignmentStates + /// as input of fitting + /// @param alignResult [in, out] The aligned result + /// @param alignMask The alignment mask (same for all measurements now) + void calculateAlignmentParameters( + const std::vector& trackAlignmentStates, + AlignmentResult& alignResult) const; + /// @brief update the detector element alignment parameters /// /// @param gctx The geometry context diff --git a/Alignment/include/ActsAlignment/Kernel/Alignment.ipp b/Alignment/include/ActsAlignment/Kernel/Alignment.ipp index 4b00839f747..86622f4ff8e 100644 --- a/Alignment/include/ActsAlignment/Kernel/Alignment.ipp +++ b/Alignment/include/ActsAlignment/Kernel/Alignment.ipp @@ -15,6 +15,7 @@ #include "Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp" #include "Acts/Utilities/detail/EigenCompat.hpp" #include "ActsAlignment/Kernel/AlignmentError.hpp" +#include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp" #include @@ -76,11 +77,6 @@ void ActsAlignment::Alignment::calculateAlignmentParameters( // The total alignment degree of freedom alignResult.alignmentDof = alignResult.idxedAlignSurfaces.size() * Acts::eAlignmentSize; - // Initialize derivative of chi2 w.r.t. alignment parameters for all tracks - Acts::DynamicVector sumChi2Derivative = - Acts::DynamicVector::Zero(alignResult.alignmentDof); - Acts::DynamicMatrix sumChi2SecondDerivative = Acts::DynamicMatrix::Zero( - alignResult.alignmentDof, alignResult.alignmentDof); // Copy the fit options fit_options_t fitOptionsWithRefSurface = fitOptions; // Calculate contribution to chi2 derivatives from all input trajectories @@ -88,7 +84,7 @@ void ActsAlignment::Alignment::calculateAlignmentParameters( alignResult.chi2 = 0; alignResult.measurementDim = 0; alignResult.numTracks = trajectoryCollection.size(); - double sumChi2ONdf = 0; + std::vector alignmentStates; for (unsigned int iTraj = 0; iTraj < trajectoryCollection.size(); iTraj++) { const auto& sourceLinks = trajectoryCollection.at(iTraj); const auto& sParameters = startParametersCollection.at(iTraj); @@ -104,6 +100,28 @@ void ActsAlignment::Alignment::calculateAlignmentParameters( continue; } const auto& alignState = evaluateRes.value(); + alignmentStates.push_back(alignState); + } + return calculateAlignmentParameters(alignmentStates, alignResult); +} + +template +void ActsAlignment::Alignment::calculateAlignmentParameters( + const std::vector& trackAlignmentStates, + AlignmentResult& alignResult) const { + // The total alignment degree of freedom + alignResult.alignmentDof = + alignResult.idxedAlignSurfaces.size() * Acts::eAlignmentSize; + // Initialize derivative of chi2 w.r.t. alignment parameters for all tracks + Acts::DynamicVector sumChi2Derivative = + Acts::DynamicVector::Zero(alignResult.alignmentDof); + Acts::DynamicMatrix sumChi2SecondDerivative = Acts::DynamicMatrix::Zero( + alignResult.alignmentDof, alignResult.alignmentDof); + alignResult.chi2 = 0; + alignResult.measurementDim = 0; + alignResult.numTracks = trackAlignmentStates.size(); + double sumChi2ONdf = 0; + for (const auto& alignState : trackAlignmentStates) { for (const auto& [rowSurface, rows] : alignState.alignedSurfaces) { const auto& [dstRow, srcRow] = rows; // Fill the results into full chi2 derivative matrix diff --git a/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp b/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp index 702167cd326..6cf61768e55 100644 --- a/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp +++ b/Alignment/include/ActsAlignment/Kernel/detail/AlignmentEngine.hpp @@ -74,6 +74,43 @@ struct TrackAlignmentState { void resetAlignmentDerivative(Acts::AlignmentToBoundMatrix& alignToBound, AlignmentMask mask); +/// @brief Helper function to calculate first and second derivative +/// of chi2 w.r.t. alignment parameters for a single track once the +/// relevant information has +/// been retrieved from the respective track. +/// Updates in-place the chi2 and the matrix/vector forming the final +/// equation. +/// @param [in,out] alignState: TrackAlignmentState to modify (in-place). +inline void finaliseTrackAlignState(TrackAlignmentState& alignState) { + // Calculate the chi2 and chi2 derivatives based on the alignment matrixs + alignState.chi2 = alignState.residual.transpose() * + alignState.measurementCovariance.inverse() * + alignState.residual; + alignState.alignmentToChi2Derivative = + Acts::DynamicVector::Zero(alignState.alignmentDof); + alignState.alignmentToChi2SecondDerivative = Acts::DynamicMatrix::Zero( + alignState.alignmentDof, alignState.alignmentDof); + // The covariance of residual + alignState.residualCovariance = Acts::DynamicMatrix::Zero( + alignState.measurementDim, alignState.measurementDim); + alignState.residualCovariance = alignState.measurementCovariance - + alignState.projectionMatrix * + alignState.trackParametersCovariance * + alignState.projectionMatrix.transpose(); + + alignState.alignmentToChi2Derivative = + 2 * alignState.alignmentToResidualDerivative.transpose() * + alignState.measurementCovariance.inverse() * + alignState.residualCovariance * + alignState.measurementCovariance.inverse() * alignState.residual; + alignState.alignmentToChi2SecondDerivative = + 2 * alignState.alignmentToResidualDerivative.transpose() * + alignState.measurementCovariance.inverse() * + alignState.residualCovariance * + alignState.measurementCovariance.inverse() * + alignState.alignmentToResidualDerivative; +} + /// /// Calculate the first and second derivative of chi2 w.r.t. alignment /// parameters for a single track @@ -266,34 +303,7 @@ TrackAlignmentState trackAlignmentState( correlation; } } - - // Calculate the chi2 and chi2 derivatives based on the alignment matrixs - alignState.chi2 = alignState.residual.transpose() * - alignState.measurementCovariance.inverse() * - alignState.residual; - alignState.alignmentToChi2Derivative = - Acts::DynamicVector::Zero(alignState.alignmentDof); - alignState.alignmentToChi2SecondDerivative = Acts::DynamicMatrix::Zero( - alignState.alignmentDof, alignState.alignmentDof); - // The covariance of residual - alignState.residualCovariance = Acts::DynamicMatrix::Zero( - alignState.measurementDim, alignState.measurementDim); - alignState.residualCovariance = alignState.measurementCovariance - - alignState.projectionMatrix * - alignState.trackParametersCovariance * - alignState.projectionMatrix.transpose(); - - alignState.alignmentToChi2Derivative = - 2 * alignState.alignmentToResidualDerivative.transpose() * - alignState.measurementCovariance.inverse() * - alignState.residualCovariance * - alignState.measurementCovariance.inverse() * alignState.residual; - alignState.alignmentToChi2SecondDerivative = - 2 * alignState.alignmentToResidualDerivative.transpose() * - alignState.measurementCovariance.inverse() * - alignState.residualCovariance * - alignState.measurementCovariance.inverse() * - alignState.alignmentToResidualDerivative; + finaliseTrackAlignState(alignState); return alignState; } diff --git a/Examples/Algorithms/AlignmentMillePede/CMakeLists.txt b/Examples/Algorithms/AlignmentMillePede/CMakeLists.txt new file mode 100644 index 00000000000..d6296822f7e --- /dev/null +++ b/Examples/Algorithms/AlignmentMillePede/CMakeLists.txt @@ -0,0 +1,17 @@ +acts_add_library( + ExamplesAlignmentMillePede + src/MillePedeAlignmentSandbox.cpp + src/ActsSolverFromMille.cpp +) +target_include_directories( + ActsExamplesAlignmentMillePede + PUBLIC $ +) +target_link_libraries( + ActsExamplesAlignmentMillePede + PUBLIC + Acts::Alignment + Acts::ExamplesFramework + Acts::ExamplesMagneticField + ActsPluginMille +) diff --git a/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/ActsSolverFromMille.hpp b/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/ActsSolverFromMille.hpp new file mode 100644 index 00000000000..6a36c63e542 --- /dev/null +++ b/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/ActsSolverFromMille.hpp @@ -0,0 +1,85 @@ +// 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/Geometry/GeometryIdentifier.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/detail/SteppingLogger.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" +#include "ActsAlignment/Kernel/Alignment.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +#include + +namespace Acts { +class MagneticFieldProvider; +} + +namespace ActsExamples { + +/// @brief Algorithm for reading Mille binaries into ACTS +/// and solving using them to run an alignment fit +/// with the built-in solver. +class ActsSolverFromMille final : public IAlgorithm { + public: + using SteppingLogger = Acts::detail::SteppingLogger; + using EndOfWorld = Acts::EndOfWorldReached; + using Stepper = Acts::EigenStepper<>; + using Propagator = Acts::Propagator; + using Fitter = Acts::KalmanFitter; + using Alignment = ActsAlignment::Alignment; + + using AlignmentParameters = + std::unordered_map; + + /// configuration + struct Config { + /// name of the mille input binary. You can choose + /// between ".root" / ".csv" / ".dat" extensions + /// to get ROOT tree / plain text / classic Millepede + /// binary outputs. All three can be read by the interface. + std::string milleInput; + // the tracking geometry to use + std::shared_ptr trackingGeometry; + // magnetic field + std::shared_ptr magneticField; + // modules to fix in the alignment to suppress global movements + std::set fixModules; + }; + + /// Constructor of the sandbox algorithm + /// @param cfg is the config struct to configure the algorithm + /// @param level is the logging level + explicit ActsSolverFromMille( + Config cfg, std::unique_ptr logger = nullptr); + + /// Framework execute method of the sandbox algorithm + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ProcessCode execute(const AlgorithmContext& ctx) const override; + ProcessCode finalize() override; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + /// configuration instance + Config m_cfg; + + /// alignment module instance - reuse as much as possible + std::shared_ptr m_align; + /// tracking geometry + std::shared_ptr m_trackingGeometry; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/MillePedeAlignmentSandbox.hpp b/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/MillePedeAlignmentSandbox.hpp new file mode 100644 index 00000000000..a69929d1449 --- /dev/null +++ b/Examples/Algorithms/AlignmentMillePede/include/ActsExamples/AlignmentMillePede/MillePedeAlignmentSandbox.hpp @@ -0,0 +1,116 @@ +// 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/Geometry/GeometryIdentifier.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/detail/SteppingLogger.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" +#include "ActsAlignment/Kernel/Alignment.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +#include + +#include "Mille/MilleFactory.h" + +namespace Acts { +class MagneticFieldProvider; +} + +namespace ActsExamples { + +/// @brief Sandbox algorithm for experimenting with +/// writing ACTS Kalman tracks to Millepede +/// and passing them to the (external) +/// Millepede alignment fit. +/// Will consume an input track collection, +/// pass the tracks through the existing +/// Kalman alignment module and then +/// write this information to a user-configured +/// Mille binary that can be read by the solver. +/// +/// You can either pass standard MC tracks and +/// look for a zero-correction, or pass +/// deliberately misaligned tracks and fit back +/// out the injected misalignment. The module +/// will **ignore** any external geometry context, +/// so injected alignment corrections in the upstream +/// job will show up as distortions here. +class MillePedeAlignmentSandbox final : public IAlgorithm { + public: + using SteppingLogger = Acts::detail::SteppingLogger; + using EndOfWorld = Acts::EndOfWorldReached; + using Stepper = Acts::EigenStepper<>; + using Propagator = Acts::Propagator; + using Fitter = Acts::KalmanFitter; + using Alignment = ActsAlignment::Alignment; + + using AlignmentParameters = + std::unordered_map; + + /// configuration + struct Config { + /// name of the mille output binary. You can choose + /// between ".root" / ".csv" / ".dat" extensions + /// to get ROOT tree / plain text / classic Millepede + /// binary outputs. All three can be read by the solver. + std::string milleOutput; + /// Input measurements collection. + std::string inputMeasurements; + /// Input tracks + std::string inputTracks; + // the tracking geometry to use + std::shared_ptr trackingGeometry; + // magnetic field + std::shared_ptr magneticField; + // modules to fix in the alignment to suppress global movements + std::set fixModules; + }; + + /// Constructor of the sandbox algorithm + /// @param cfg is the config struct to configure the algorithm + /// @param level is the logging level + explicit MillePedeAlignmentSandbox( + Config cfg, std::unique_ptr logger = nullptr); + + /// Framework execute method of the sandbox algorithm + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ProcessCode execute(const AlgorithmContext& ctx) const override; + ProcessCode finalize() override; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + /// configuration instance + Config m_cfg; + + /// measurement container containing the measurements on the input tracks + /// below + ReadDataHandle m_inputMeasurements{this, + "InputMeasurements"}; + /// tracks to use for the alignment + ReadDataHandle m_inputTracks{this, "InputTracks"}; + + /// alignment module instance - reuse as much as possible + std::shared_ptr m_align; + /// tracking geometry + std::shared_ptr m_trackingGeometry; + /// the Mille record instance for writing our alignment info. + std::unique_ptr m_milleOut = nullptr; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/AlignmentMillePede/src/ActsSolverFromMille.cpp b/Examples/Algorithms/AlignmentMillePede/src/ActsSolverFromMille.cpp new file mode 100644 index 00000000000..4312e417138 --- /dev/null +++ b/Examples/Algorithms/AlignmentMillePede/src/ActsSolverFromMille.cpp @@ -0,0 +1,149 @@ +// 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/AlignmentMillePede/ActsSolverFromMille.hpp" + +#include "Acts/Definitions/Alignment.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Material/MaterialInteraction.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsAlignment/Kernel/Alignment.hpp" +#include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" +#include "ActsPlugins/Mille/ActsToMille.hpp" + +#include + +#include + +namespace ActsExamples { + +ActsSolverFromMille::ActsSolverFromMille( + Config cfg, std::unique_ptr logger) + : IAlgorithm("ActsSolverFromMille", std::move(logger)), + m_cfg(std::move(cfg)) { + // retrieve tracking geo + m_trackingGeometry = m_cfg.trackingGeometry; + + // instantiate the alignment tool instance + Acts::Navigator::Config navcfg{m_cfg.trackingGeometry}; + navcfg.resolvePassive = false; + navcfg.resolveMaterial = true; + navcfg.resolveSensitive = true; + Acts::Navigator navigator(navcfg, + this->logger().cloneWithSuffix("Navigator")); + Stepper stepper{m_cfg.magneticField}; + m_align = std::make_shared( + Fitter(Propagator(stepper, Acts::Navigator(navcfg)))); +} + +ProcessCode ActsSolverFromMille::execute( + const AlgorithmContext& /*ClangShutUp*/) const { + // this algorithm will not do anything at event-time + return ProcessCode::SUCCESS; +} + +ProcessCode ActsSolverFromMille::finalize() { + auto milleReader = Mille::spawnMilleReader(m_cfg.milleInput); + auto openStat = milleReader->open(m_cfg.milleInput); + if (!openStat) { + ACTS_FATAL("Failed to read the mille binary " << m_cfg.milleInput); + return ProcessCode::ABORT; + } + + // Assign indices to the alignable surfaces + + // We wish to have a relation between alignment parameter indices and real + // geometry. The unordered_map does not give us this - so perform a manual + // sorting. + std::vector> + sortedGeo; + ActsAlignment::AlignmentResult alignResult; + + sortedGeo.insert(sortedGeo.end(), + m_trackingGeometry->geoIdSurfaceMap().begin(), + m_trackingGeometry->geoIdSurfaceMap().end()); + std::sort( + sortedGeo.begin(), sortedGeo.end(), + [&](const std::pair& lhs, + const std::pair& + rhs) { return (lhs.first.layer() < rhs.first.layer()); }); + + const Acts::Surface* firstSurf = nullptr; + unsigned int iSurface = 0; + for (auto& [geoID, surface] : sortedGeo) { + // only consider sensitive surfaces + if (geoID.sensitive() == 0) { + continue; + } + // use the first sensitive surface as trajectory reference in the kalman + if (firstSurf == nullptr) { + firstSurf = surface; + } + if (!m_cfg.fixModules.contains(geoID)) { + alignResult.idxedAlignSurfaces.emplace(surface, iSurface); + iSurface++; + } + } + + ACTS_INFO("Performing alignment fit on collected Mille records"); + std::vector alignmentStates; + ActsAlignment::detail::TrackAlignmentState state; + std::size_t iRec = 0; + while (ActsPlugins::ActsToMille::unpackMilleRecord( + *milleReader, state, alignResult.idxedAlignSurfaces) == + Mille::MilleDecoder::ReadResult::OK) { + if (++iRec % 10000 == 0) { + ACTS_INFO(" Reading input record " << iRec); + } + alignmentStates.push_back(state); + } + + /// TODO: Should try a local iteration without track state info. + /// Can use the linearised info in the Track Alignment State + /// to calculate approximate track parameter & residual updates + /// and then repeat the solution. As in Millepede, probably + /// safe to keep the "big matrix" and only update the right hand side. + m_align->calculateAlignmentParameters(alignmentStates, alignResult); + + /// in a real experiment, the results would be written out + /// and stored e.g. in a DB file for further use / validation. + /// For this initial demo, we just print them out. + + std::cout << "Performed internal alignment. " << std::endl; + std::cout << std::setw(16) << " Tracks used: " << alignmentStates.size() + << std::endl; + std::cout << std::setw(16) + << " avg Chi2/NDF = " << alignResult.averageChi2ONdf << std::endl; + std::cout << std::setw(16) << " Chi2 = " << alignResult.chi2 << std::endl; + std::cout << std::setw(16) << " delta Chi2 = " << alignResult.deltaChi2 + << std::endl; + std::cout << std::setw(16) << " Alignment parameter updates: " << std::endl; + std::vector parLabels{"dx", "dy", "dz", "rx", "ry", "rz"}; + for (auto [surface, index] : alignResult.idxedAlignSurfaces) { + std::cout << std::setw(20) << " Surface with geo ID " + << surface->geometryId() << ": " << std::endl; + for (std::size_t i = 0; i < Acts::eAlignmentSize; ++i) { + std::size_t row = Acts::eAlignmentSize * index + i; + std::cout << std::setw(20) << parLabels[i] << " = " << std::setw(10) + << alignResult.deltaAlignmentParameters(row) << std::setw(6) + << " +/- " << std::setw(10) + << std::sqrt(alignResult.alignmentCovariance(row, row)) + << std::endl; + ; + } + } + + return ProcessCode::SUCCESS; +} + +} // namespace ActsExamples diff --git a/Examples/Algorithms/AlignmentMillePede/src/MillePedeAlignmentSandbox.cpp b/Examples/Algorithms/AlignmentMillePede/src/MillePedeAlignmentSandbox.cpp new file mode 100644 index 00000000000..10f43718419 --- /dev/null +++ b/Examples/Algorithms/AlignmentMillePede/src/MillePedeAlignmentSandbox.cpp @@ -0,0 +1,178 @@ +// 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/AlignmentMillePede/MillePedeAlignmentSandbox.hpp" + +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Material/MaterialInteraction.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/MeasurementCalibration.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" +#include "ActsPlugins/Mille/ActsToMille.hpp" + +#include + +namespace ActsExamples { + +MillePedeAlignmentSandbox::MillePedeAlignmentSandbox( + Config cfg, std::unique_ptr logger) + : IAlgorithm("MillePedeAlignmentSandbox", std::move(logger)), + m_cfg(std::move(cfg)) { + // initialize our handles + if (m_cfg.inputMeasurements.empty()) { + throw std::invalid_argument("Missing input measurement collection"); + } + if (m_cfg.inputTracks.empty()) { + throw std::invalid_argument("Missing input track collection"); + } + m_inputTracks.initialize(m_cfg.inputTracks); + m_inputMeasurements.initialize(m_cfg.inputMeasurements); + + // retrieve tracking geo + m_trackingGeometry = m_cfg.trackingGeometry; + + // instantiate the alignment tool instance + Acts::Navigator::Config navcfg{m_cfg.trackingGeometry}; + navcfg.resolvePassive = false; + navcfg.resolveMaterial = true; + navcfg.resolveSensitive = true; + Acts::Navigator navigator(navcfg, + this->logger().cloneWithSuffix("Navigator")); + Stepper stepper{m_cfg.magneticField}; + m_align = std::make_shared( + Fitter(Propagator(stepper, Acts::Navigator(navcfg)))); + + /// spawn a Mille binary to record our alignment inputs. + /// You can specify root / csv / dat extensions for + /// ROOT NTuple / plain text (careful: large files) or C-binary + /// storage. + /// The file will be automatically closed upon deletion. + m_milleOut = Mille::spawnMilleRecord(m_cfg.milleOutput); +} + +ProcessCode MillePedeAlignmentSandbox::execute( + const AlgorithmContext& ctx) const { + using TrackFitterOptions = + Acts::KalmanFitterOptions; + + // Read input data + const auto& measurements = m_inputMeasurements(ctx); + const auto& tracks = m_inputTracks(ctx); + + // Assign indices to the alignable surfaces + + // We wish to have a relation between alignment parameter indices and real + // geometry. The unordered_map does not give us this - so perform a manual + // sorting. + std::vector> + sortedGeo; + sortedGeo.insert(sortedGeo.end(), + m_trackingGeometry->geoIdSurfaceMap().begin(), + m_trackingGeometry->geoIdSurfaceMap().end()); + std::sort( + sortedGeo.begin(), sortedGeo.end(), + [&](const std::pair& lhs, + const std::pair& + rhs) { return (lhs.first.layer() < rhs.first.layer()); }); + + std::unordered_map indexedAlignSurfaces; + const Acts::Surface* firstSurf = nullptr; + unsigned int iSurface = 0; + for (auto& [geoID, surface] : sortedGeo) { + // only consider sensitive surfaces + if (geoID.sensitive() == 0) { + continue; + } + // use the first sensitive surface as trajectory reference in the kalman + if (firstSurf == nullptr) { + firstSurf = surface; + } + if (!m_cfg.fixModules.contains(geoID)) { + indexedAlignSurfaces.emplace(surface, iSurface); + iSurface++; + } + } + + // Dirty hack: Overwrite the geometry context to remove knowledge + // of injected alignment shift. + // Detector will look misaligned and the injected correction can be fitted + // out. + // TODO: Replace if an appropriate non-hacky tool becomes available. + Acts::GeometryContext dummyGeoCtx = + Acts::GeometryContext::dangerouslyDefaultConstruct(); + + // Pile of boilerplate code to get the Kalman fitter for the + // alignment module ready to go + IndexSourceLinkSurfaceAccessor slack{*m_trackingGeometry}; + Acts::KalmanFitterExtensions extensions; + ActsExamples::PassThroughCalibrator pcalibrator; + MeasurementCalibratorAdapter calibrator(pcalibrator, measurements); + extensions.calibrator.connect<&MeasurementCalibratorAdapter::calibrate>( + &calibrator); + Acts::GainMatrixUpdater kfUpdater; + Acts::GainMatrixSmoother kfSmoother; + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kfUpdater); + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()>( + &kfSmoother); + extensions.surfaceAccessor + .connect<&ActsExamples::IndexSourceLinkSurfaceAccessor::operator()>( + &slack); + TrackFitterOptions kfOptions( + dummyGeoCtx, ctx.magFieldContext, ctx.calibContext, extensions, + Acts::PropagatorPlainOptions(dummyGeoCtx, ctx.magFieldContext), + firstSurf); + + // loop over tracks in the event + std::vector trackSourceLinks; + for (const auto& track : tracks) { + // for starting parameters for the re-fit, ask the + // existing CKF track + Acts::BoundTrackParameters refPar = track.createParametersAtReference(); + + // Collect source links from this track + trackSourceLinks.clear(); + trackSourceLinks.reserve(track.nTrackStates()); + for (const auto& state : track.trackStates()) { + trackSourceLinks.push_back(state.getUncalibratedSourceLink()); + } + + // get the TrackAlignmentState using the existing + // alignment class. This will compute the needed + // residuals and derivatives. + auto aliStates = m_align->evaluateTrackAlignmentState( + dummyGeoCtx, trackSourceLinks, refPar, kfOptions, indexedAlignSurfaces, + ActsAlignment::AlignmentMask::All // use this to restrict alignment + // degrees of freedom if desired + ); + + // and, if successful, dump the information into our Mille record. + if (aliStates.ok()) { + const ActsAlignment::detail::TrackAlignmentState& state = *aliStates; + ActsPlugins::ActsToMille::dumpToMille(state, m_milleOut.get()); + } + } + + return ProcessCode::SUCCESS; +} +ProcessCode MillePedeAlignmentSandbox::finalize() { + m_milleOut.reset(); // ensure that we do the final write of our output + // before subsequent algos finalise. + return ProcessCode::SUCCESS; +} + +} // namespace ActsExamples diff --git a/Examples/Algorithms/CMakeLists.txt b/Examples/Algorithms/CMakeLists.txt index 8114adecb1b..70a833e65a0 100644 --- a/Examples/Algorithms/CMakeLists.txt +++ b/Examples/Algorithms/CMakeLists.txt @@ -15,6 +15,7 @@ add_subdirectory_if(Jets ACTS_BUILD_EXAMPLES_FASTJET) add_subdirectory(TruthTracking) add_subdirectory(Vertexing) add_subdirectory_if(Alignment ACTS_BUILD_ALIGNMENT) +add_subdirectory_if(AlignmentMillePede ACTS_BUILD_PLUGIN_MILLE) add_subdirectory(Utilities) add_subdirectory(AmbiguityResolution) add_subdirectory(HelloWorld) diff --git a/Examples/Configs/telescope-seeding-config.json b/Examples/Configs/telescope-seeding-config.json new file mode 100644 index 00000000000..afbf27dcd95 --- /dev/null +++ b/Examples/Configs/telescope-seeding-config.json @@ -0,0 +1,38 @@ +[ + { + "layer": 12, + "volume": 1 + }, + { + "layer": 10, + "volume": 1 + }, + { + "layer": 8, + "volume": 1 + }, + { + "layer": 14, + "volume": 1 + }, + { + "layer": 18, + "volume": 1 + }, + { + "layer": 16, + "volume": 1 + }, + { + "layer": 2, + "volume": 1 + }, + { + "layer": 4, + "volume": 1 + }, + { + "layer": 6, + "volume": 1 + } +] diff --git a/Examples/Detectors/TelescopeDetector/CMakeLists.txt b/Examples/Detectors/TelescopeDetector/CMakeLists.txt index 4216ee5a598..12d71e34d8b 100644 --- a/Examples/Detectors/TelescopeDetector/CMakeLists.txt +++ b/Examples/Detectors/TelescopeDetector/CMakeLists.txt @@ -3,6 +3,7 @@ acts_add_library( src/TelescopeDetector.cpp src/TelescopeDetectorElement.cpp src/BuildTelescopeDetector.cpp + src/AlignedTelescopeDetector.cpp ) target_include_directories( diff --git a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/AlignedTelescopeDetector.hpp b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/AlignedTelescopeDetector.hpp new file mode 100644 index 00000000000..f480a1599c5 --- /dev/null +++ b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/AlignedTelescopeDetector.hpp @@ -0,0 +1,25 @@ +// 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 "ActsExamples/DetectorCommons/Aligned.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp" + +namespace ActsExamples { +/// Define the aligned telescope detector element and factory type +using AlignedTelescopeDetectorElement = Aligned; + +class AlignedTelescopeDetector : public TelescopeDetector { + public: + using Config = TelescopeDetector::Config; + explicit AlignedTelescopeDetector(const Config& cfg); +}; + +} // namespace ActsExamples diff --git a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp index bf227dfa812..ef3974cc73e 100644 --- a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp +++ b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp @@ -10,6 +10,7 @@ #include "Acts/Geometry/TrackingGeometry.hpp" #include "Acts/Utilities/AxisDefinitions.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp" #include #include @@ -26,6 +27,7 @@ enum class TelescopeSurfaceType { /// Global method to build the telescope tracking geometry /// /// @param gctx is the detector element dependent geometry context +/// @param factory is the factory responsible for creating the detector elements /// @param detectorStore is the store for the detector element /// @param positions are the positions of different layers in the longitudinal /// direction @@ -41,6 +43,7 @@ enum class TelescopeSurfaceType { /// parallel to std::unique_ptr buildTelescopeDetector( const Acts::GeometryContext& gctx, + const TelescopeDetectorElement::DetectorElementFactory& factory, std::vector>& detectorStore, const std::vector& positions, diff --git a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetector.hpp b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetector.hpp index 679e3bb066f..341a7451e79 100644 --- a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetector.hpp +++ b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetector.hpp @@ -42,6 +42,11 @@ class TelescopeDetector : public Detector { std::unique_ptr buildGeant4DetectorConstruction( const Geant4ConstructionOptions& options) const override; + protected: + struct NoBuildTag {}; + explicit TelescopeDetector( + const Config& cfg, NoBuildTag /*unused*/); // used for aligned version + private: Config m_cfg; }; diff --git a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp index e2dcdf0ac7f..6da4c171292 100644 --- a/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp +++ b/Examples/Detectors/TelescopeDetector/include/ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp @@ -38,6 +38,17 @@ class TelescopeDetectorElement : public Acts::SurfacePlacementBase { /// The current interval of validity unsigned int iov = 0; }; + using identifier_type = unsigned long long; + using identifier_diff = long long; + using Identifier = identifier_type; + + using DetectorElementFactory = + std::function( + const Acts::Transform3&, + std::variant, + std::shared_ptr>, + double, std::shared_ptr, + std::vector>&)>; /// Constructor for single sided detector element /// - bound to a Plane Surface @@ -47,7 +58,7 @@ class TelescopeDetectorElement : public Acts::SurfacePlacementBase { /// @param thickness is the module thickness /// @param material is the (optional) Surface material associated to it TelescopeDetectorElement( - std::shared_ptr transform, + const Identifier ID, std::shared_ptr transform, std::shared_ptr pBounds, double thickness, std::shared_ptr material = nullptr); @@ -59,7 +70,7 @@ class TelescopeDetectorElement : public Acts::SurfacePlacementBase { /// @param thickness is the module thickness /// @param material is the (optional) Surface material associated to it TelescopeDetectorElement( - std::shared_ptr transform, + const Identifier ID, std::shared_ptr transform, std::shared_ptr dBounds, double thickness, std::shared_ptr material = nullptr); @@ -81,13 +92,18 @@ class TelescopeDetectorElement : public Acts::SurfacePlacementBase { /// /// @note this is called from the surface().transform() in the PROXY mode const Acts::Transform3& localToGlobalTransform( - const Acts::GeometryContext& gctx) const final; + const Acts::GeometryContext& gctx) const override; + + /// Return local to global transform associated with this detector element + const Acts::Transform3& nominalTransform() const; /// Return the nominal local to global transform /// /// @note the geometry context will hereby be ignored const Acts::Transform3& nominalTransform( - const Acts::GeometryContext& gctx) const; + const Acts::GeometryContext& /*gctx*/) const { + return nominalTransform(); + } /// Return local to global transform associated with this identifier /// @@ -102,7 +118,12 @@ class TelescopeDetectorElement : public Acts::SurfacePlacementBase { /// Is the detector element a sensitive element bool isSensitive() const final { return true; } + /// The identifier of the detector element + Identifier identifier() const; + private: + // The element identifier + Identifier m_elementIdentifier; /// the transform for positioning in 3D space std::shared_ptr m_elementTransform = nullptr; // the aligned transforms @@ -141,8 +162,8 @@ inline const Acts::Transform3& TelescopeDetectorElement::localToGlobalTransform( return nominalTransform(gctx); } -inline const Acts::Transform3& TelescopeDetectorElement::nominalTransform( - const Acts::GeometryContext& /*gctx*/) const { +inline const Acts::Transform3& TelescopeDetectorElement::nominalTransform() + const { return *m_elementTransform; } @@ -161,4 +182,8 @@ TelescopeDetectorElement::alignedTransforms() const { return m_alignedTransforms; } +inline TelescopeDetectorElement::Identifier +TelescopeDetectorElement::identifier() const { + return m_elementIdentifier; +} } // namespace ActsExamples diff --git a/Examples/Detectors/TelescopeDetector/src/AlignedTelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/AlignedTelescopeDetector.cpp new file mode 100644 index 00000000000..52dd3cf3b0f --- /dev/null +++ b/Examples/Detectors/TelescopeDetector/src/AlignedTelescopeDetector.cpp @@ -0,0 +1,55 @@ +// 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/TelescopeDetector/AlignedTelescopeDetector.hpp" + +#include "ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp" + +namespace ActsExamples { + +AlignedTelescopeDetector::AlignedTelescopeDetector(const Config& cfg) + : TelescopeDetector(cfg, NoBuildTag{}) { + m_nominalGeometryContext = + Acts::GeometryContext::dangerouslyDefaultConstruct(); + + // Set the detector element factory + auto alignedDetectorElementFactory = + [](const Acts::Transform3& transform, + std::variant, + std::shared_ptr> + bounds, + double thickness, + std::shared_ptr material, + std::vector>& + detStore) { + auto ID = + static_cast(detStore.size()); + std::shared_ptr detElem; + if (bounds.index() == 0) { + detElem = std::make_shared( + ID, std::make_shared(transform), + std::get>(bounds), + thickness, std::move(material)); + } else { + detElem = std::make_shared( + ID, std::make_shared(transform), + std::get>(bounds), + thickness, std::move(material)); + } + detStore.push_back(detElem); + return detElem; + }; + m_trackingGeometry = buildTelescopeDetector( + m_nominalGeometryContext, alignedDetectorElementFactory, m_detectorStore, + cfg.positions, cfg.stereos, cfg.offsets, cfg.bounds, cfg.thickness, + static_cast(cfg.surfaceType), + static_cast(cfg.binValue)); +} + +} // namespace ActsExamples diff --git a/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp index 9a5d31bd81c..6b22be10bc4 100644 --- a/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp +++ b/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp @@ -29,7 +29,6 @@ #include "Acts/Surfaces/SurfaceArray.hpp" #include "Acts/Surfaces/SurfacePlacementBase.hpp" #include "Acts/Utilities/Logger.hpp" -#include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp" #include #include @@ -38,6 +37,7 @@ std::unique_ptr ActsExamples::buildTelescopeDetector( const Acts::GeometryContext& gctx, + const TelescopeDetectorElement::DetectorElementFactory& factory, std::vector>& detectorStore, const std::vector& positions, @@ -92,13 +92,11 @@ ActsExamples::buildTelescopeDetector( // Create the detector element std::shared_ptr detElement = nullptr; if (surfaceType == TelescopeSurfaceType::Plane) { - detElement = std::make_shared( - std::make_shared(trafo), pBounds, 1._um, - surfaceMaterial); + detElement = + factory(trafo, pBounds, 1._um, surfaceMaterial, detectorStore); } else { - detElement = std::make_shared( - std::make_shared(trafo), rBounds, 1._um, - surfaceMaterial); + detElement = + factory(trafo, rBounds, 1._um, surfaceMaterial, detectorStore); } // Get the surface auto surface = detElement->surface().getSharedPtr(); diff --git a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp index 58b035a8d90..e6f7a00b8b2 100644 --- a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp +++ b/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp @@ -10,6 +10,7 @@ #include "Acts/Geometry/GeometryContext.hpp" #include "ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetectorElement.hpp" #include @@ -40,12 +41,62 @@ TelescopeDetector::TelescopeDetector(const Config& cfg) m_nominalGeometryContext = Acts::GeometryContext::dangerouslyDefaultConstruct(); + auto detectorElementFactory = + [](const Acts::Transform3& transform, + std::variant, + std::shared_ptr> + bounds, + double thickness, + std::shared_ptr material, + std::vector>& + detStore) { + auto ID = + static_cast(detStore.size()); + std::shared_ptr detElem; + if (bounds.index() == 0) { + detElem = std::make_shared( + ID, std::make_shared(transform), + std::get>(bounds), + thickness, std::move(material)); + } else { + detElem = std::make_shared( + ID, std::make_shared(transform), + std::get>(bounds), + thickness, std::move(material)); + } + detStore.push_back(detElem); + return detElem; + }; m_trackingGeometry = buildTelescopeDetector( - m_nominalGeometryContext, m_detectorStore, m_cfg.positions, m_cfg.stereos, - m_cfg.offsets, m_cfg.bounds, m_cfg.thickness, - static_cast(m_cfg.surfaceType), + m_nominalGeometryContext, detectorElementFactory, m_detectorStore, + m_cfg.positions, m_cfg.stereos, m_cfg.offsets, m_cfg.bounds, + m_cfg.thickness, static_cast(m_cfg.surfaceType), static_cast(m_cfg.binValue)); } +TelescopeDetector::TelescopeDetector(const Config& cfg, NoBuildTag /*unused*/) + : Detector(Acts::getDefaultLogger("TelescopeDetector", cfg.logLevel)), + m_cfg(cfg) { + if (m_cfg.surfaceType > 1) { + throw std::invalid_argument( + "The surface type could either be 0 for plane surface or 1 for disc " + "surface."); + } + if (m_cfg.binValue > 2) { + throw std::invalid_argument("The axis value could only be 0, 1, or 2."); + } + // Check if the bounds values are valid + if (m_cfg.surfaceType == 1 && m_cfg.bounds[0] >= m_cfg.bounds[1]) { + throw std::invalid_argument( + "The minR should be smaller than the maxR for disc surface bounds."); + } + + if (m_cfg.positions.size() != m_cfg.stereos.size()) { + throw std::invalid_argument( + "The number of provided positions must match the number of " + "provided stereo angles."); + } +} + } // namespace ActsExamples diff --git a/Examples/Detectors/TelescopeDetector/src/TelescopeDetectorElement.cpp b/Examples/Detectors/TelescopeDetector/src/TelescopeDetectorElement.cpp index c5efa55b59f..c824730e5a8 100644 --- a/Examples/Detectors/TelescopeDetector/src/TelescopeDetectorElement.cpp +++ b/Examples/Detectors/TelescopeDetector/src/TelescopeDetectorElement.cpp @@ -14,10 +14,11 @@ namespace ActsExamples { TelescopeDetectorElement::TelescopeDetectorElement( - std::shared_ptr transform, + const Identifier ID, std::shared_ptr transform, std::shared_ptr pBounds, double thickness, std::shared_ptr material) - : m_elementTransform(std::move(transform)), + : m_elementIdentifier(ID), + m_elementTransform(std::move(transform)), m_elementSurface( Acts::Surface::makeShared(pBounds, *this)), m_elementThickness(thickness), @@ -28,10 +29,11 @@ TelescopeDetectorElement::TelescopeDetectorElement( } TelescopeDetectorElement::TelescopeDetectorElement( - std::shared_ptr transform, + const Identifier ID, std::shared_ptr transform, std::shared_ptr dBounds, double thickness, std::shared_ptr material) - : m_elementTransform(std::move(transform)), + : m_elementIdentifier(ID), + m_elementTransform(std::move(transform)), m_elementSurface( Acts::Surface::makeShared(dBounds, *this)), m_elementThickness(thickness), diff --git a/Examples/Scripts/Python/millepede_alignment.py b/Examples/Scripts/Python/millepede_alignment.py new file mode 100644 index 00000000000..38c701c5b31 --- /dev/null +++ b/Examples/Scripts/Python/millepede_alignment.py @@ -0,0 +1,319 @@ +#!/usr/bin/env python3 + +import os +import argparse +import pathlib +from pathlib import Path + +import acts +import acts.examples + +from acts.examples import ( + TelescopeDetector, + AlignedTelescopeDetector, + Sequencer, + StructureSelector, + RandomNumbers, + GaussianVertexGenerator, +) + +from acts.examples.alignment import ( + AlignmentDecorator, + GeoIdAlignmentStore, + AlignmentGeneratorGlobalShift, +) +from acts.examples.alignmentmillepede import ( + MillePedeAlignmentSandbox, + ActsSolverFromMille, +) + +from acts.examples.simulation import ( + MomentumConfig, + EtaConfig, + PhiConfig, + ParticleConfig, + ParticleSelectorConfig, + addParticleGun, + addFatras, + addDigitization, + addDigiParticleSelection, +) +from acts.examples.reconstruction import ( + addSeeding, + CkfConfig, + addCKFTracks, + TrackSelectorConfig, + SeedingAlgorithm, + TrackSelectorConfig, + addSeeding, + SeedingAlgorithm, + SeedFinderConfigArg, + SeedFinderOptionsArg, + SeedingAlgorithm, + CkfConfig, + addCKFTracks, + TrackSelectorConfig, +) + + +# Helper to instantiate a telescope detector. +# The square sensors are oriented in the global +# y-direction and cover the x-z plane. +# You can change the number of layers by resizing +# the "bounds", "stereos" and "positions" arrays accordingly. +# +# By default, will be digitised as 25 x 100 pixel grid. +# +# The "layer" field of the geo ID of this detector +# will move in steps of 2 for each module (2,4,6,..,18 for 9 layers). +# Everything is located in volume "1". +# +# In the alignment, at least 4 layers are expected (layer ID 8 will +# be fixed as the alignment reference). +def getTelescopeDetector(misaligned: bool = False): + bounds = [200, 200] + positions = [30, 60, 90, 120, 150, 180, 210, 240, 270] + stereos = [0] * len(positions) + if not misaligned: + detector = TelescopeDetector( + bounds=bounds, positions=positions, stereos=stereos, binValue=1 + ) + else: + detector = AlignedTelescopeDetector( + bounds=bounds, positions=positions, stereos=stereos, binValue=1 + ) + + return detector + + +# Add the alignment sandbox algorithm +def addAlignmentSandbox( + s: Sequencer, + trackingGeometry: acts.TrackingGeometry, + magField: acts.MagneticFieldProvider, + fixModules: set, + inputMeasurements: str = "measurements", + inputTracks: str = "ckf_tracks", + logLevel: acts.logging.Level = acts.logging.WARNING, + milleOutput: str = "MilleBinary.root", +): + sandbox = MillePedeAlignmentSandbox( + level=logLevel, + milleOutput=milleOutput, + inputMeasurements=inputMeasurements, + inputTracks=inputTracks, + trackingGeometry=trackingGeometry, + magneticField=magField, + fixModules=fixModules, + ) + s.addAlgorithm(sandbox) + return s + + +def addSolverFromMille( + s: Sequencer, + trackingGeometry: acts.TrackingGeometry, + magField: acts.MagneticFieldProvider, + fixModules: set, + logLevel: acts.logging.Level = acts.logging.WARNING, + milleInput: str = "MilleBinary.root", +): + + solver = ActsSolverFromMille( + level=logLevel, + milleInput=milleInput, + trackingGeometry=trackingGeometry, + magneticField=magField, + fixModules=fixModules, + ) + s.addAlgorithm(solver) + return s + + +u = acts.UnitConstants + +# Can also use zero-field - but not healthy for track covariance. +# field = acts.NullBField() +field = acts.ConstantBField(acts.Vector3(0, 0, 2 * acts.UnitConstants.T)) + + +parser = argparse.ArgumentParser( + description="MillePede alignment demo with the Telescope Detector" +) +parser.add_argument( + "--output", + "-o", + help="Output directory", + type=pathlib.Path, + default=pathlib.Path.cwd() / "mpali_output", +) +parser.add_argument("--events", "-n", help="Number of events", type=int, default=2000) +parser.add_argument("--skip", "-s", help="Number of events", type=int, default=0) + + +args = parser.parse_args() + +outputDir = args.output +# ensure out output dir exists +os.makedirs(outputDir, exist_ok=True) + +# Instantiate the telescope detector - with alignment enabled +detector = getTelescopeDetector(misaligned=True) +trackingGeometry = detector.trackingGeometry() +decorators = detector.contextDecorators() + +# inject a known misalignment. + +# Misalign the second tracking layer (ID = 4). +layerToBump = acts.GeometryIdentifier(layer=4, volume=1, sensitive=1) + +# shift this layer by 200 microns in the global Z direction +leShift = AlignmentGeneratorGlobalShift() +leShift.shift = acts.Vector3(0, 0, 200.0e-3) + +# now add some boilerplate code to make this happen +alignDecoConfig = AlignmentDecorator.Config() +alignDecoConfig.nominalStore = GeoIdAlignmentStore( + StructureSelector(trackingGeometry).selectedTransforms( + acts.GeometryContext.dangerouslyDefaultConstruct(), layerToBump + ) +) +alignDecoConfig.iovGenerators = [((0, 10000000), leShift)] +alignDeco = AlignmentDecorator(alignDecoConfig, acts.logging.WARNING) +contextDecorators = [alignDeco] + +# decide on at least on detector module to fix in place +# as a reference for the alignment. +# By default, fix the innermost layer. +fixModules = {acts.GeometryIdentifier(layer=2, volume=1, sensitive=1)} + + +# More Boilerplate code - for setting up the sequence + +rnd = RandomNumbers(seed=42) +s = Sequencer( + events=args.events, + skip=args.skip, + numThreads=1, + outputDir=str(outputDir), +) + +# Add a context with the alignment shift - sim, digi and +# initial reco will "see" the distorted detector +s.addContextDecorator(alignDeco) + +# Run particle gun and fire some muons at our telescope +addParticleGun( + s, + MomentumConfig( + 10 * u.GeV, + 100 * u.GeV, + transverse=True, + ), + EtaConfig(-0.3, 0.3), + # aim roughly along +Y... + PhiConfig(60 * u.degree, 120 * u.degree), + ParticleConfig(1, acts.PdgParticle.eMuon, randomizeCharge=True), + vtxGen=GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4(5.0 * u.mm, 0.0 * u.mm, 5.0 * u.mm, 0.0 * u.ns), + ), + multiplicity=1, + rnd=rnd, +) +# fast sim +addFatras( + s, + trackingGeometry, + field, + enableInteractions=True, + outputDirRoot=None, + outputDirCsv=None, + outputDirObj=None, + rnd=rnd, +) + +# digitise with the default 25x100 pixel config +srcdir = Path(__file__).resolve().parent.parent.parent.parent +addDigitization( + s, + trackingGeometry, + field, + digiConfigFile=srcdir / "Examples/Configs/telescope-digi-smearing-config.json", + outputDirRoot=None, + outputDirCsv=None, + rnd=rnd, +) +addDigiParticleSelection( + s, + ParticleSelectorConfig( + measurements=(3, None), + removeNeutral=True, + ), +) + +# Run grid seeder +addSeeding( + s, + trackingGeometry, + field, + # Settings copied from existing snippet, slightly adapted + seedFinderConfigArg=SeedFinderConfigArg( + r=(20 * u.mm, 200 * u.mm), + deltaR=(1 * u.mm, 300 * u.mm), + collisionRegion=(-250 * u.mm, 250 * u.mm), + z=(-100 * u.mm, 100 * u.mm), + maxSeedsPerSpM=1, + sigmaScattering=5, + radLengthPerSeed=0.1, + minPt=0.5 * u.GeV, + impactMax=3 * u.mm, + ), + # why do we need to specify this here again? Not taken from event context? + seedFinderOptionsArg=SeedFinderOptionsArg(bFieldInZ=2 * u.T), + seedingAlgorithm=SeedingAlgorithm.GridTriplet, + initialSigmas=[ + 3 * u.mm, + 3 * u.mm, + 1 * u.degree, + 1 * u.degree, + 0 * u.e / u.GeV, + 1 * u.ns, + ], + initialSigmaQoverPt=0.1 * u.e / u.GeV, + initialSigmaPtRel=0.1, + initialVarInflation=[1.0] * 6, + # This file should be adapted if you add layers and want to include + # them in the seeding. + geoSelectionConfigFile=srcdir / "Examples/Configs/telescope-seeding-config.json", + outputDirRoot=None, +) + +# Add CKF track finding +addCKFTracks( + s, + trackingGeometry, + field, + TrackSelectorConfig(), + CkfConfig( + chi2CutOffMeasurement=150.0, + chi2CutOffOutlier=250.0, + numMeasurementsCutOff=50, + seedDeduplication=(True), + stayOnSeed=True, + ), + outputDirRoot=None, + writePerformance=False, + writeTrackSummary=False, +) + +# And add our alignment sandbox +addAlignmentSandbox( + s, trackingGeometry, field, fixModules, milleOutput=outputDir / "MyBinary.root" +) +# And finally read back and solve in ACTS +addSolverFromMille( + s, trackingGeometry, field, fixModules, milleInput=outputDir / "MyBinary.root" +) + +s.run() diff --git a/Plugins/Mille/CMakeLists.txt b/Plugins/Mille/CMakeLists.txt index 64b88487489..e628bce57f2 100644 --- a/Plugins/Mille/CMakeLists.txt +++ b/Plugins/Mille/CMakeLists.txt @@ -1,5 +1,6 @@ acts_add_library( PluginMille + src/Helpers.cpp src/ActsToMille.cpp ACTS_INCLUDE_FOLDER include/ActsPlugins ) diff --git a/Plugins/Mille/include/ActsPlugins/Mille/ActsToMille.hpp b/Plugins/Mille/include/ActsPlugins/Mille/ActsToMille.hpp index f5a93a2ec95..16cd49b2ef0 100644 --- a/Plugins/Mille/include/ActsPlugins/Mille/ActsToMille.hpp +++ b/Plugins/Mille/include/ActsPlugins/Mille/ActsToMille.hpp @@ -10,19 +10,43 @@ #include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp" +#include "Mille/IMilleReader.h" +#include "Mille/MilleDecoder.h" #include "Mille/MilleRecord.h" namespace ActsPlugins::ActsToMille { -/// The MilleRecord is the primary interface for -/// writing out alignment fit inputs +/// The MilleRecord is Millepede's interface for +/// writing out alignment fit inputs. It can be instantiated +/// using Mille::spawnMilleRecord(desired_file_name), +/// provided by Mille/MilleFactory.h using Mille::MilleRecord; -/// Placeholder method to test header -/// discovery and linkage. -/// Using TrackAlignmentState as a potential -/// candidate for a future (internal) interface. -void dumpToMille(const ActsAlignment::detail::TrackAlignmentState&, +/// @brief Dump a Kalman track encoded as a TrackAlignmentState into +/// a Mille record. +/// @param state: Alignment state to dump. +/// @param record: Mille record to write to - should be valid pointer +/// Note: Not very efficient - we have to "un-fit" the kalman track. +/// Used for R&D, recommending the GBL track model (under development) +// for production use. +void dumpToMille(const ActsAlignment::detail::TrackAlignmentState& state, MilleRecord* record); +/// @brief read one record (= track or (constrained) track pair) from +/// a Mille binary into the equivalent matrices of a TrackAlignmentState. +/// Allows to use Mille to collect tracks across multiple events and +/// align them with the ACTS solver, and to validate the outputs of dumpToMille. +/// @param reader: A Mille Reader, connected to a valid input file. +/// @param targetState: The TrackAlignmentState to populate. +/// @param idxedAlignSurfaces: [optional]: Indexed alignment surfaces from the geometry. If passed, +/// the internal `alignedSurfaces` member of the state will be configured to +/// link back to the correct surfaces. +/// @return a ReadResult enum with 3 possible states to indicate the outcome- ok / end-of-file / read-error. +/// The targetState will only be modified if the result is 'ok'. +Mille::MilleDecoder::ReadResult unpackMilleRecord( + Mille::IMilleReader& reader, + ActsAlignment::detail::TrackAlignmentState& targetState, + const std::unordered_map& + idxedAlignSurfaces); + } // namespace ActsPlugins::ActsToMille diff --git a/Plugins/Mille/include/ActsPlugins/Mille/Helpers.hpp b/Plugins/Mille/include/ActsPlugins/Mille/Helpers.hpp new file mode 100644 index 00000000000..30a2fb13be5 --- /dev/null +++ b/Plugins/Mille/include/ActsPlugins/Mille/Helpers.hpp @@ -0,0 +1,44 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +namespace ActsPlugins::ActsToMille { + +/// @brief: Regularise a covariance matrix for decomposition into Mille. +/// Required especially when running fits with non-timing detectors. +/// @param inputCov: Input covariance matrix, to be regularised +/// @param conditionCutOff: Lowest value (relative to leading EV) to +// clamp small / negative eigenvalues to +/// @param removeLargeLeading: If true, will also regularise a single huge +/// leading eigenvalue, symptomatic of time-fits +/// on non-timing-detectors +/// @param return A new matrix, which has been regularised by clamping +/// eigenvalues to the iterval [conditionCutOff x max_EV, max_EV], +/// where max_EV is either the leading eigenvalue or, if +/// removeLargeLeading is set, either the first or the second +/// leading EV (second if largest is more than 100 times the +/// second) +Acts::DynamicMatrix regulariseCovariance(const Acts::DynamicMatrix& inputCov, + double conditionCutOff = 1e-10, + bool removeLargeLeading = true); + +/// Calculates the solution X to the matrix equation C = (A + X)^-1, +/// for a poorly conditioned C and known A. Required to decompose the ACTS +/// Kalman covariance matrix into a series of Mille pseudo-measurements. +/// Uses cholesky factorisation of C to solve (1 - CA) = CX +/// to avoid directly inverting C. +/// @param target: The target matrix C to be decomposed - expected to be symmetric +/// and positive (semi)definite +/// @param existing_sol: The partial existing solution A - expected to be symmetric +/// @return the missing piece X solving the above equation. +Acts::DynamicMatrix getInverseComplement( + const Acts::DynamicMatrix& target, const Acts::DynamicMatrix& existing_sol); + +} // namespace ActsPlugins::ActsToMille diff --git a/Plugins/Mille/src/ActsToMille.cpp b/Plugins/Mille/src/ActsToMille.cpp index 0b4cb3d85d1..53b5cc4dc24 100644 --- a/Plugins/Mille/src/ActsToMille.cpp +++ b/Plugins/Mille/src/ActsToMille.cpp @@ -8,21 +8,329 @@ #include "ActsPlugins/Mille/ActsToMille.hpp" +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Alignment.hpp" +#include "ActsPlugins/Mille/Helpers.hpp" + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "Mille/MilleDecoder.h" + namespace ActsPlugins::ActsToMille { -/// this is a placeholder, which will be replaced by actual functionality -void dumpToMille(const ActsAlignment::detail::TrackAlignmentState&, +namespace { + +/// for the ACTS-internal indices, start counting at 0 +unsigned long internalIndexSurfToParam(unsigned long surfaceIndex, + unsigned long dofIndex) { + return surfaceIndex * Acts::eAlignmentSize + dofIndex; +} +/// for global alignment parameters, start counting at 1 (fotran convention) for +/// Mille +unsigned long globalIndexSurfToParam(unsigned long surfaceIndex, + unsigned long dofIndex) { + return surfaceIndex * Acts::eAlignmentSize + dofIndex + 1; +} + +} // namespace + +void dumpToMille(const ActsAlignment::detail::TrackAlignmentState& state, MilleRecord* record) { - // call into Mille to test for linkage - if (record) { - // write a single dummy measurement to test the I/O functionality + // check for a valid record + if (record == nullptr) { + std::cerr << " Missing Mille record " << std::endl; + return; + } + + // prepare the vectors to interface to Mille + std::vector localIndices(state.trackParametersDim, 0); + std::vector globalIndices(state.alignmentDof, 0.); + std::vector localDeriv(state.trackParametersDim, 0.); + std::vector globalDeriv(state.alignmentDof, 0.); + + // prepare the track parameter index array (always the same) + std::iota(localIndices.begin(), localIndices.end(), 1); + + // map the alignment parameter labels. + // Important: Millepede expects indices to start with 1 + std::vector> aliParLocalToGlobal; + for (auto& [surf, indices] : state.alignedSurfaces) { + auto& [globalSurfIndex, localStartIndex] = indices; + for (std::size_t iPar = 0; iPar < Acts::eAlignmentSize; ++iPar) { + aliParLocalToGlobal.emplace_back( + internalIndexSurfToParam(localStartIndex, iPar), + globalIndexSurfToParam(globalSurfIndex, iPar)); + } + } + + /// 1) write out the local measurements on the surfaces and their direct + /// derivatives. This will populate the upper / left three quadrants of the + /// alignment matrix, including direct correlations between alignment and + /// track parameters. + /// TODO: Add explicit diagonalisation for correlated (stereo) measurements. + for (std::size_t iMeas = 0; iMeas < state.measurementDim; ++iMeas) { + // arrange the global parameters correctly + for (auto& [srcGlobal, destGlobal] : aliParLocalToGlobal) { + // index for each global derivative + globalIndices[srcGlobal] = destGlobal; + // value for each global derivative + globalDeriv[srcGlobal] = + state.alignmentToResidualDerivative(iMeas, srcGlobal); + } + // local derivatives due to measurement uncertainties + for (std::size_t iTrkPar = 0; iTrkPar < state.trackParametersDim; + ++iTrkPar) { + localDeriv[iTrkPar] = state.projectionMatrix(iMeas, iTrkPar); + } + // write a measurement to the ongoing Mille record. record->addData( - -0.005f, // measurement - 0.008f, // uncertainty - std::vector{1}, // indices of local parameters - std::vector{1.}, // derivates by local parameters - std::vector{42}, // indices of alignment parameters - std::vector{-1.}); // derivatives by alignment parameters + // residual + state.residual(iMeas), + // sigma + std::sqrt(state.measurementCovariance(iMeas, iMeas)), + // local parameter indices + localIndices, + // local derivatives + localDeriv, globalIndices, globalDeriv); + } + + /// 2) Write out additional pseudo-measurements representing the (local) track + /// parameter correlations from the Kalman fit (linearisation point). + /// These enter the bottom right quadrant of the alignment matrix and + /// represent the additional contributions to the track fit chi2 + /// arising from the correlations between the parameters on different + /// surfaces encoded in the Kalman track model. + /// Step 1) already added the terms arising from measurement uncertainties, + /// now we need to extend this to the full Kalman covariance. + /// This requires "unfitting" the tracks, making this a rather slow + /// and numerically tricky step. + + /// compute the part of the weight matrix arising from Step 1). + /// This is already present in the Mille record and should not + /// be duplicated + const Acts::DynamicMatrix weightMatMeasurements = + state.projectionMatrix.transpose() * + state.measurementCovariance.inverse() * state.projectionMatrix; + + // regularise the (full) Kalman covariance. This is needed to stabilise + // poorly constrained directions (usually: time) + const Acts::DynamicMatrix regularisedCov = + regulariseCovariance(state.trackParametersCovariance); + + // now we can get the piece of the weight matrix not already covered by + // the measurement uncertainties + const Acts::DynamicMatrix correlationTerm = + getInverseComplement(regularisedCov, weightMatMeasurements); + + // Decompose the matrix we need to add into a sum of rank-1 matrices, + // C_add = sum (lambda_i v_i v_i^T), which can be interpreted + // as pseudo-measurements with sigma_i 1/sqrt(lambda_i) and local derivatives + // v_i. This relies on C_add being symmetric positive (semi)definite. + Eigen::SelfAdjointEigenSolver eigenSolver( + correlationTerm); + if (eigenSolver.info() != Eigen::Success) { + std::cout << " FAILED to find decompose correlation term" << std::endl; + return; + } + const Acts::DynamicVector eigenVals = eigenSolver.eigenvalues(); + const Acts::DynamicMatrix eigenVecs = eigenSolver.eigenvectors(); + + // no dependence on global parameters - these terms only enter the + // track covariance sub-matrix of the alignment problem (bottom right + // quadrant) + globalDeriv.clear(); + globalIndices.clear(); + + /// convert each EV to a pseudo-measurement + for (std::size_t iMeas = 0; iMeas < state.trackParametersDim; ++iMeas) { + // fill the local derivatives from the current eigenvector + for (std::size_t iTrkPar = 0; iTrkPar < localDeriv.size(); ++iTrkPar) { + localDeriv[iTrkPar] = eigenVecs(iTrkPar, iMeas); + } + // and write a pseudo-measurement to Mille. + record->addData( + // residual == 0 for pseudo-measurements + 0, + // EV == weight = 1/sigma^2 + 1. / std::sqrt(eigenVals(iMeas)), + // local parameter indices + localIndices, + // local derivatives + localDeriv, globalIndices, globalDeriv); + } + // track is fully written - end the record in Mille + record->writeRecord(); +} + +Mille::MilleDecoder::ReadResult unpackMilleRecord( + Mille::IMilleReader& reader, + ActsAlignment::detail::TrackAlignmentState& targetState, + const std::unordered_map& + idxedAlignSurfaces) { + /// book a decoder + Mille::MilleDecoder decoder; + // vector to hold the extracted measurements + std::vector measurements; + // attempt to decode the next record from the binary + auto res = decoder.decode(reader, measurements); + + // if we are EoF or encountered an error, return the result. + if (res != Mille::MilleDecoder::ReadResult::OK) { + return res; + } + + // initialise the components of the target state which we will touch + targetState.measurementDim = 0; + + // Still here - we got a valid record! Let's write it into our target state. + // In the following, we emulate what MillePede-II is doing internally. + // This is somewhat approximate, as we do not run any of the cleaning / + // conditioning performed by MP-II. + + // Step 1: Parameter discovery + // Goal: Identify all existing parameters and assign internal indices. + int firstLocal = 99999; + int lastLocal = 0; + std::set seenGlobalLabels; + std::set seenSurfaceLabels; + // a measurement with a residual of identical zero is typically a + // correlation constraint encoded as pseudo-measurement + targetState.measurementDim = + std::count_if(measurements.begin(), measurements.end(), + [](const Mille::MilleMeasurement& measurement) { + return measurement.measurement != 0; + }); + + // discover labels in use + for (const Mille::MilleMeasurement& measurement : measurements) { + auto [minLabel, maxLabel] = std::minmax_element( + measurement.localLabels.begin(), measurement.localLabels.end()); + firstLocal = std::min(firstLocal, *minLabel); + lastLocal = std::max(lastLocal, *maxLabel); + seenGlobalLabels.insert(measurement.globalLabels.begin(), + measurement.globalLabels.end()); + } + targetState.trackParametersDim = lastLocal - firstLocal + 1; + targetState.alignmentDof = seenGlobalLabels.size(); + + for (int label : seenGlobalLabels) { + if ((label - 1) % Acts::eAlignmentSize == 0) { + seenSurfaceLabels.insert((label - 1) / Acts::eAlignmentSize); + } + }; + + /// the trackAlignmentState uses an internal indexing for alignment + /// parameters - so remap the indices to replicate this internal logic. + std::map globalToInternal; + + /// try to map indices from the global indexed surface list, if we have it. + unsigned long iExtra = 0; + for (auto [surface, index] : idxedAlignSurfaces) { + if (seenSurfaceLabels.contains(index)) { + targetState.alignedSurfaces.emplace(surface, + std::make_pair(index, iExtra++)); + } + } + + // now use this to fill our internal "global to local" mapping function + for (auto [surface, indices] : targetState.alignedSurfaces) { + auto [globIx, intIx] = indices; + for (std::size_t iAli = 0; iAli < Acts::eAlignmentSize; ++iAli) { + globalToInternal.emplace(globalIndexSurfToParam(globIx, iAli), + internalIndexSurfToParam(intIx, iAli)); + } } + + // Now we have the needed information to initialise our matrices + targetState.measurementCovariance = Acts::DynamicMatrix::Zero( + targetState.measurementDim, targetState.measurementDim); + + targetState.projectionMatrix = Acts::DynamicMatrix::Zero( + targetState.measurementDim, targetState.trackParametersDim); + + targetState.alignmentToResidualDerivative = Acts::DynamicMatrix::Zero( + targetState.measurementDim, targetState.alignmentDof); + + targetState.trackParametersCovariance = Acts::DynamicMatrix::Zero( + targetState.trackParametersDim, targetState.trackParametersDim); + + targetState.residual = Acts::DynamicVector::Zero(targetState.measurementDim); + + /// Second loop - fill the matrices + + std::size_t iMeas = 0; + for (const auto& measurement : measurements) { + // need distinction: Measurement on surface vs. correlation term. + // The reason is that ACTS only counts surface measurements, and stores + // the correlation information directly in the track parameter covariance. + // MillePede considers constraints to be additional measurements. + // A MillePede pseudomeasurement has residual 0 and no global derivatives. + bool isMeasurementOnSurface = (measurement.measurement != 0 || + !measurement.globalDerivatives.empty()); + // surface measurements populate the residual vector and measurement + // covariance matrix + if (isMeasurementOnSurface) { + targetState.residual(iMeas) = measurement.measurement; + targetState.measurementCovariance(iMeas, iMeas) = + measurement.uncertainty * measurement.uncertainty; + } + // loop over all track parameters affecting this measurement + for (std::size_t iLoc = 0; iLoc < measurement.localLabels.size(); ++iLoc) { + // find out where to book it in the ACTS matrix + unsigned int localIndex = measurement.localLabels[iLoc] - firstLocal; + // if we are a surface measurement, fill the projection matrix + if (isMeasurementOnSurface) + targetState.projectionMatrix(iMeas, localIndex) = + measurement.localDerivatives[iLoc]; + // now fill the covariance matrix by looping over all products of (local) + // derivatives + for (std::size_t jLoc = 0; jLoc < measurement.localLabels.size(); + ++jLoc) { + // again determine where to book the column index + unsigned int localIndex2 = measurement.localLabels[jLoc] - firstLocal; + // and update the covariance. + targetState.trackParametersCovariance(localIndex, localIndex2) += + measurement.localDerivatives[iLoc] * + measurement.localDerivatives[jLoc] / measurement.uncertainty / + measurement.uncertainty; + } + } + // loop over global (= alignment) parameters affecting the measurement + for (std::size_t iGlob = 0; iGlob < measurement.globalLabels.size(); + ++iGlob) { + // find out where to book - here we need to map to the ACTS track-level + // indexing scheme + int internalAliIndex = globalToInternal[measurement.globalLabels[iGlob]]; + // and update the alignment-to-residual derivative matrix. + targetState.alignmentToResidualDerivative(iMeas, internalAliIndex) = + measurement.globalDerivatives[iGlob]; + } + // increment the measurement-on-surface index every time we finish + // processing one. + if (isMeasurementOnSurface) + ++iMeas; + } + + /// (carefully) invert the covariance - upstairs, we filled it as a weight + /// matrix + auto solver = targetState.trackParametersCovariance.ldlt(); + targetState.trackParametersCovariance = + solver.solve(Acts::DynamicMatrix::Identity( + targetState.trackParametersDim, targetState.trackParametersDim)); + + /// and calculate the dependent members + /// (first and second derivatives, chi2) of the state. + ActsAlignment::detail::finaliseTrackAlignState(targetState); + + return res; } + } // namespace ActsPlugins::ActsToMille diff --git a/Plugins/Mille/src/Helpers.cpp b/Plugins/Mille/src/Helpers.cpp new file mode 100644 index 00000000000..fd9ce920e3f --- /dev/null +++ b/Plugins/Mille/src/Helpers.cpp @@ -0,0 +1,69 @@ +// 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 "ActsPlugins/Mille/Helpers.hpp" + +#include + +#include +#include + +namespace ActsPlugins::ActsToMille { + +Acts::DynamicMatrix regulariseCovariance(const Acts::DynamicMatrix& inputCov, + double conditionCutOff, + bool removeLargeLeading) { + Acts::DynamicMatrix out = + Acts::DynamicMatrix::Zero(inputCov.rows(), inputCov.cols()); + + /// add a tiny diagonal matrix for additional stabilisation + auto eigensolver = Eigen::SelfAdjointEigenSolver( + inputCov + + 1.e-10 * Acts::DynamicMatrix::Identity(inputCov.rows(), inputCov.cols())); + + if (eigensolver.info() != Eigen::Success) { + std::cout << " FAILED to find eigenvec" << std::endl; + return out; + } + auto eigenVals = eigensolver.eigenvalues(); + auto eigenVecs = eigensolver.eigenvectors(); + + std::size_t maxIndex = eigenVals.size() - 1; + double lambdaMax = eigenVals(eigenVals.size() - 1); + // check for a huge leading eigenvalue - this happens when the time coordinate + // is unconstrained + if (removeLargeLeading && lambdaMax > 100 * eigenVals(maxIndex - 1)) { + --maxIndex; + lambdaMax = eigenVals(maxIndex); + } + double lambdaMin = conditionCutOff * lambdaMax; + // clamp the EV to the permitted interval + for (Eigen::Index i = 0; i < eigenVals.size(); ++i) { + out(i, i) = std::clamp(eigenVals(i), lambdaMin, lambdaMax); + } + // then return the regularised matrix as V D V^T + out = out * eigenVecs.transpose(); + out = eigenVecs * out; + return out; +} + +Acts::DynamicMatrix getInverseComplement( + const Acts::DynamicMatrix& target, + const Acts::DynamicMatrix& existing_sol) { + Acts::DynamicMatrix Rhs = + Acts::DynamicMatrix::Identity(target.rows(), target.cols()) - + target * existing_sol; + // call decomposition from Eigen (prefer over llt for semi-def. matrices) + auto LDLT = target.ldlt(); + Acts::DynamicMatrix solution = LDLT.solve(Rhs); + // finally, symmetrise to correct for floating point effects + solution = 0.5 * (solution + solution.transpose()); + return solution; +} + +} // end namespace ActsPlugins::ActsToMille diff --git a/Python/Examples/CMakeLists.txt b/Python/Examples/CMakeLists.txt index fe7f18d4f28..97f9642ade9 100644 --- a/Python/Examples/CMakeLists.txt +++ b/Python/Examples/CMakeLists.txt @@ -148,6 +148,7 @@ target_link_libraries( # Add all plugin related examples bindings conditionally to their build flag add_examples_binding_if(Alignment "Acts::ExamplesAlignment;Acts::ExamplesDetectorsCommon" ACTS_BUILD_ALIGNMENT FALSE) +add_examples_binding_if(AlignmentMillePede Acts::ExamplesAlignmentMillePede ACTS_BUILD_PLUGIN_MILLE FALSE) add_examples_binding_if(DD4hep Acts::ExamplesDetectorDD4hep ACTS_BUILD_EXAMPLES_DD4HEP TRUE) add_examples_binding_if(EDM4hep Acts::ExamplesIoEDM4hep ACTS_BUILD_EXAMPLES_EDM4HEP FALSE) add_examples_binding_if(Geant4 diff --git a/Python/Examples/src/Detector.cpp b/Python/Examples/src/Detector.cpp index 69a9d462527..61049aa40dc 100644 --- a/Python/Examples/src/Detector.cpp +++ b/Python/Examples/src/Detector.cpp @@ -16,6 +16,7 @@ #include "ActsExamples/Framework/IContextDecorator.hpp" #include "ActsExamples/GenericDetector/AlignedGenericDetector.hpp" #include "ActsExamples/GenericDetector/GenericDetector.hpp" +#include "ActsExamples/TelescopeDetector/AlignedTelescopeDetector.hpp" #include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp" #include "ActsExamples/Utilities/Options.hpp" #include "ActsPython/Utilities/Helpers.hpp" @@ -91,6 +92,12 @@ void addDetector(py::module& mex) { ACTS_PYTHON_STRUCT(c, positions, stereos, offsets, bounds, thickness, surfaceType, binValue, materialDecorator, logLevel); } + { + auto ad = py::class_>( + mex, "AlignedTelescopeDetector") + .def(py::init()); + } } } // namespace ActsPython diff --git a/Python/Examples/src/plugins/AlignmentMillePede.cpp b/Python/Examples/src/plugins/AlignmentMillePede.cpp new file mode 100644 index 00000000000..3992094d874 --- /dev/null +++ b/Python/Examples/src/plugins/AlignmentMillePede.cpp @@ -0,0 +1,30 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#include "ActsExamples/AlignmentMillePede/ActsSolverFromMille.hpp" +#include "ActsExamples/AlignmentMillePede/MillePedeAlignmentSandbox.hpp" +#include "ActsPython/Utilities/Macros.hpp" + +#include +#include + +namespace py = pybind11; + +using namespace Acts; +using namespace ActsExamples; +using namespace ActsPython; + +PYBIND11_MODULE(ActsExamplesPythonBindingsAlignmentMillePede, m) { + ACTS_PYTHON_DECLARE_ALGORITHM(MillePedeAlignmentSandbox, m, + "MillePedeAlignmentSandbox", milleOutput, + inputMeasurements, inputTracks, + trackingGeometry, magneticField, fixModules); + ACTS_PYTHON_DECLARE_ALGORITHM(ActsSolverFromMille, m, "ActsSolverFromMille", + milleInput, trackingGeometry, magneticField, + fixModules); +} diff --git a/Tests/CommonHelpers/include/ActsTests/CommonHelpers/AlignmentHelpers.hpp b/Tests/CommonHelpers/include/ActsTests/CommonHelpers/AlignmentHelpers.hpp new file mode 100644 index 00000000000..32830c978a0 --- /dev/null +++ b/Tests/CommonHelpers/include/ActsTests/CommonHelpers/AlignmentHelpers.hpp @@ -0,0 +1,289 @@ +// 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/. + +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/VectorMultiTrajectory.hpp" +#include "Acts/EventData/detail/TestSourceLink.hpp" +#include "Acts/Geometry/CuboidVolumeBounds.hpp" +#include "Acts/Geometry/CuboidVolumeBuilder.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/LayerArrayCreator.hpp" +#include "Acts/Geometry/LayerCreator.hpp" +#include "Acts/Geometry/PlaneLayer.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Geometry/TrackingGeometryBuilder.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/MagneticField/ConstantBField.hpp" +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" +#include "Acts/Material/ISurfaceMaterial.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Navigator.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/StraightLineStepper.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Surfaces/SurfaceArray.hpp" +#include "Acts/TrackFitting/GainMatrixSmoother.hpp" +#include "Acts/TrackFitting/GainMatrixUpdater.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" +#include "Acts/Utilities/CalibrationContext.hpp" +#include "ActsTests/CommonHelpers/DetectorElementStub.hpp" +#include "ActsTests/CommonHelpers/MeasurementsCreator.hpp" +#include "ActsTests/CommonHelpers/PredefinedMaterials.hpp" + +#include +#include + +namespace ActsTests::AlignmentUtils { +using namespace Acts::UnitLiterals; + +using StraightPropagator = + Acts::Propagator; +using ConstantFieldStepper = Acts::EigenStepper<>; +using ConstantFieldPropagator = + Acts::Propagator; + +using KalmanUpdater = Acts::GainMatrixUpdater; +using KalmanSmoother = Acts::GainMatrixSmoother; +using KalmanFitterType = + Acts::KalmanFitter; + +/// helper struct packaging commonly used data members for alignment +/// unit tests +struct aliTestUtils { + KalmanUpdater kfUpdater; + KalmanSmoother kfSmoother; + const Acts::GeometryContext geoCtx = + Acts::GeometryContext::dangerouslyDefaultConstruct(); + const Acts::MagneticFieldContext magCtx; + const Acts::CalibrationContext calCtx; + std::normal_distribution normalDist{0., 1.}; + std::default_random_engine rng{42}; + // detector resolutions + const MeasurementResolution resPixel = {MeasurementType::eLoc01, + {30_um, 50_um}}; + const MeasurementResolutionMap resolutions = { + {Acts::GeometryIdentifier(), resPixel}, + }; +}; + +// Create a test context + +Acts::KalmanFitterExtensions getExtensions( + aliTestUtils& utils) { + Acts::KalmanFitterExtensions extensions; + extensions.calibrator.connect<&Acts::detail::Test::testSourceLinkCalibrator< + Acts::VectorMultiTrajectory>>(); + extensions.updater + .connect<&KalmanUpdater::operator()>( + &utils.kfUpdater); + extensions.smoother + .connect<&KalmanSmoother::operator()>( + &utils.kfSmoother); + return extensions; +} + +/// +/// @brief Construct a telescope-like detector +/// +struct TelescopeDetector { + /// Default constructor for the Cubic tracking geometry + /// + /// @param gctx the geometry context for this geometry at building time + explicit TelescopeDetector( + std::reference_wrapper gctx) + : geoContext(gctx) { + // Construct the rotation + rotation.col(0) = Acts::Vector3(0, 0, -1); + rotation.col(1) = Acts::Vector3(0, 1, 0); + rotation.col(2) = Acts::Vector3(1, 0, 0); + + // Boundaries of the surfaces + rBounds = std::make_shared(0.1_m, 0.1_m); + + // Material of the surfaces + Acts::MaterialSlab matProp(makeSilicon(), 80_um); + + surfaceMaterial = + std::make_shared(matProp); + } + + /// + /// Call operator to build the standard cubic tracking geometry + /// + std::shared_ptr operator()() { + using namespace Acts::UnitLiterals; + + unsigned int nLayers = 6; + std::vector positions = {-500_mm, -300_mm, -100_mm, + 100_mm, 300_mm, 500_mm}; + auto length = positions.back() - positions.front(); + + std::vector layers(nLayers); + for (unsigned int i = 0; i < nLayers; ++i) { + // The transform + Acts::Translation3 trans(0., 0., positions[i]); + Acts::Transform3 trafo(rotation * trans); + auto detElement = std::make_shared( + trafo, rBounds, 1._um, surfaceMaterial); + // The surface is not right!!! + auto surface = detElement->surface().getSharedPtr(); + // Add it to the event store + detectorStore.push_back(std::move(detElement)); + std::unique_ptr surArray( + new Acts::SurfaceArray(surface)); + // The layer thickness should not be too large + layers[i] = Acts::PlaneLayer::create( + trafo, rBounds, std::move(surArray), + 1._mm); // Associate the layer to the surface + auto mutableSurface = const_cast(surface.get()); + mutableSurface->associateLayer(*layers[i]); + } + + // The volume transform + Acts::Translation3 transVol(0, 0, 0); + Acts::Transform3 trafoVol(rotation * transVol); + auto boundsVol = std::make_shared( + rBounds->halfLengthX() + 10._mm, rBounds->halfLengthY() + 10._mm, + length + 10._mm); + + Acts::LayerArrayCreator::Config lacConfig; + Acts::LayerArrayCreator layArrCreator( + lacConfig, + Acts::getDefaultLogger("LayerArrayCreator", Acts::Logging::INFO)); + Acts::LayerVector layVec; + for (unsigned int i = 0; i < nLayers; i++) { + layVec.push_back(layers[i]); + } + + // Create the layer array + std::unique_ptr layArr(layArrCreator.layerArray( + geoContext, layVec, positions.front() - 2._mm, positions.back() + 2._mm, + Acts::BinningType::arbitrary, Acts::AxisDirection::AxisX)); + + // Build the tracking volume + auto trackVolume = std::make_shared( + trafoVol, boundsVol, nullptr, std::move(layArr), nullptr, + Acts::MutableTrackingVolumeVector{}, "Telescope"); + + return std::make_shared(trackVolume); + } + + Acts::RotationMatrix3 rotation = Acts::RotationMatrix3::Identity(); + std::shared_ptr rBounds = nullptr; + std::shared_ptr surfaceMaterial = nullptr; + + std::vector> detectorStore; + + std::reference_wrapper geoContext; +}; + +// Construct a straight-line propagator. +StraightPropagator makeStraightPropagator( + std::shared_ptr geo) { + Acts::Navigator::Config cfg{std::move(geo)}; + cfg.resolvePassive = false; + cfg.resolveMaterial = true; + cfg.resolveSensitive = true; + Acts::Navigator navigator(cfg); + Acts::StraightLineStepper stepper; + return StraightPropagator(stepper, std::move(navigator)); +} + +// Construct a propagator using a constant magnetic field along z. +ConstantFieldPropagator makeConstantFieldPropagator( + std::shared_ptr geo, double bz, + std::unique_ptr logger) { + Acts::Navigator::Config cfg{std::move(geo)}; + cfg.resolvePassive = false; + cfg.resolveMaterial = true; + cfg.resolveSensitive = true; + Acts::Navigator navigator(cfg, logger->cloneWithSuffix("Nav")); + auto field = + std::make_shared(Acts::Vector3(0.0, 0.0, bz)); + ConstantFieldStepper stepper(std::move(field)); + return ConstantFieldPropagator(std::move(stepper), std::move(navigator), + logger->cloneWithSuffix("Prop")); +} + +// Construct initial track parameters. +Acts::BoundTrackParameters makeParameters(aliTestUtils& utils) { + using namespace Acts::UnitLiterals; + // create covariance matrix from reasonable standard deviations + Acts::BoundVector stddev; + stddev[Acts::eBoundLoc0] = 100_um; + stddev[Acts::eBoundLoc1] = 100_um; + stddev[Acts::eBoundTime] = 25_ns; + stddev[Acts::eBoundPhi] = 0.5_degree; + stddev[Acts::eBoundTheta] = 0.5_degree; + stddev[Acts::eBoundQOverP] = 1 / 100_GeV; + Acts::BoundMatrix cov = stddev.cwiseProduct(stddev).asDiagonal(); + + auto loc0 = 0. + stddev[Acts::eBoundLoc0] * utils.normalDist(utils.rng); + auto loc1 = 0. + stddev[Acts::eBoundLoc1] * utils.normalDist(utils.rng); + auto t = 42_ns + stddev[Acts::eBoundTime] * utils.normalDist(utils.rng); + auto phi = 0_degree + stddev[Acts::eBoundPhi] * utils.normalDist(utils.rng); + auto theta = + 90_degree + stddev[Acts::eBoundTheta] * utils.normalDist(utils.rng); + auto qOverP = + 1_e / 1_GeV + stddev[Acts::eBoundQOverP] * utils.normalDist(utils.rng); + + // define a track in the transverse plane along x + Acts::Vector4 mPos4(-1_m, loc0, loc1, t); + + return Acts::BoundTrackParameters::createCurvilinear( + mPos4, phi, theta, qOverP, cov, Acts::ParticleHypothesis::pion()); +} + +struct KalmanFitterInputTrajectory { + // The source links + std::vector sourceLinks; + // The start parameters + std::optional startParameters; +}; + +/// +/// Function to create trajectories for kalman fitter +/// +std::vector createTrajectories( + std::shared_ptr geo, + std::size_t nTrajectories, aliTestUtils& util) { + // simulation propagator + const auto simPropagator = makeStraightPropagator(std::move(geo)); + + std::vector trajectories; + trajectories.reserve(nTrajectories); + + for (unsigned int iTrack = 0; iTrack < nTrajectories; iTrack++) { + auto start = makeParameters(util); + // Launch and collect - the measurements + auto measurements = + createMeasurements(simPropagator, util.geoCtx, util.magCtx, start, + util.resolutions, util.rng); + + // Extract measurements from result of propagation. + KalmanFitterInputTrajectory traj; + traj.startParameters = start; + traj.sourceLinks = measurements.sourceLinks; + + trajectories.push_back(std::move(traj)); + } + return trajectories; +} +} // namespace ActsTests::AlignmentUtils + +/// diff --git a/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp b/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp index df943fae8c0..ea2e83330f0 100644 --- a/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp +++ b/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp @@ -8,272 +8,27 @@ #include -#include "Acts/Definitions/Units.hpp" #include "Acts/EventData/detail/TestSourceLink.hpp" -#include "Acts/Geometry/CuboidVolumeBounds.hpp" #include "Acts/Geometry/CuboidVolumeBuilder.hpp" #include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Geometry/LayerArrayCreator.hpp" -#include "Acts/Geometry/LayerCreator.hpp" -#include "Acts/Geometry/PlaneLayer.hpp" -#include "Acts/Geometry/TrackingGeometry.hpp" -#include "Acts/Geometry/TrackingGeometryBuilder.hpp" -#include "Acts/Geometry/TrackingVolume.hpp" -#include "Acts/MagneticField/ConstantBField.hpp" -#include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" -#include "Acts/Material/ISurfaceMaterial.hpp" -#include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/Navigator.hpp" -#include "Acts/Propagator/Propagator.hpp" -#include "Acts/Propagator/StraightLineStepper.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" -#include "Acts/Surfaces/RectangleBounds.hpp" -#include "Acts/Surfaces/SurfaceArray.hpp" -#include "Acts/TrackFitting/GainMatrixSmoother.hpp" -#include "Acts/TrackFitting/GainMatrixUpdater.hpp" -#include "Acts/TrackFitting/KalmanFitter.hpp" -#include "Acts/Utilities/CalibrationContext.hpp" #include "ActsAlignment/Kernel/Alignment.hpp" -#include "ActsTests/CommonHelpers/DetectorElementStub.hpp" +#include "ActsTests/CommonHelpers/AlignmentHelpers.hpp" #include "ActsTests/CommonHelpers/FloatComparisons.hpp" -#include "ActsTests/CommonHelpers/MeasurementsCreator.hpp" -#include "ActsTests/CommonHelpers/PredefinedMaterials.hpp" -#include #include -namespace { - using namespace Acts; -using namespace ActsAlignment; using namespace ActsTests; +using namespace ActsTests::AlignmentUtils; using namespace Acts::detail::Test; -using namespace Acts::UnitLiterals; - -using StraightPropagator = Propagator; -using ConstantFieldStepper = EigenStepper<>; -using ConstantFieldPropagator = Propagator; - -using KalmanUpdater = GainMatrixUpdater; -using KalmanSmoother = GainMatrixSmoother; -using KalmanFitterType = - KalmanFitter; - -KalmanUpdater kfUpdater; -KalmanSmoother kfSmoother; - -// Create a test context -const auto geoCtx = GeometryContext::dangerouslyDefaultConstruct(); -const MagneticFieldContext magCtx; -const CalibrationContext calCtx; - -std::normal_distribution normalDist(0., 1.); -std::default_random_engine rng(42); - -KalmanFitterExtensions getExtensions() { - KalmanFitterExtensions extensions; - extensions.calibrator - .connect<&testSourceLinkCalibrator>(); - extensions.updater.connect<&KalmanUpdater::operator()>( - &kfUpdater); - extensions.smoother - .connect<&KalmanSmoother::operator()>(&kfSmoother); - return extensions; -} - -/// -/// @brief Construct a telescope-like detector -/// -struct TelescopeDetector { - /// Default constructor for the Cubic tracking geometry - /// - /// @param gctx the geometry context for this geometry at building time - explicit TelescopeDetector(std::reference_wrapper gctx) - : geoContext(gctx) { - // Construct the rotation - rotation.col(0) = Vector3(0, 0, -1); - rotation.col(1) = Vector3(0, 1, 0); - rotation.col(2) = Vector3(1, 0, 0); - - // Boundaries of the surfaces - rBounds = std::make_shared(0.1_m, 0.1_m); - - // Material of the surfaces - MaterialSlab matProp(makeSilicon(), 80_um); - - surfaceMaterial = std::make_shared(matProp); - } - - /// - /// Call operator to build the standard cubic tracking geometry - /// - std::shared_ptr operator()() { - using namespace UnitLiterals; - - unsigned int nLayers = 6; - std::vector positions = {-500_mm, -300_mm, -100_mm, - 100_mm, 300_mm, 500_mm}; - auto length = positions.back() - positions.front(); - - std::vector layers(nLayers); - for (unsigned int i = 0; i < nLayers; ++i) { - // The transform - Translation3 trans(0., 0., positions[i]); - Transform3 trafo(rotation * trans); - auto detElement = std::make_shared( - trafo, rBounds, 1._um, surfaceMaterial); - // The surface is not right!!! - auto surface = detElement->surface().getSharedPtr(); - // Add it to the event store - detectorStore.push_back(std::move(detElement)); - auto surArray = std::make_unique(surface); - // The layer thickness should not be too large - layers[i] = - PlaneLayer::create(trafo, rBounds, std::move(surArray), - 1._mm); // Associate the layer to the surface - auto mutableSurface = const_cast(surface.get()); - mutableSurface->associateLayer(*layers[i]); - } - - // The volume transform - Translation3 transVol(0, 0, 0); - Transform3 trafoVol(rotation * transVol); - auto boundsVol = std::make_shared( - rBounds->halfLengthX() + 10._mm, rBounds->halfLengthY() + 10._mm, - length + 10._mm); - - LayerArrayCreator::Config lacConfig; - LayerArrayCreator layArrCreator( - lacConfig, getDefaultLogger("LayerArrayCreator", Logging::INFO)); - LayerVector layVec; - for (unsigned int i = 0; i < nLayers; i++) { - layVec.push_back(layers[i]); - } - - // Create the layer array - std::unique_ptr layArr(layArrCreator.layerArray( - geoContext, layVec, positions.front() - 2._mm, positions.back() + 2._mm, - BinningType::arbitrary, AxisDirection::AxisX)); - - // Build the tracking volume - auto trackVolume = std::make_shared( - trafoVol, boundsVol, nullptr, std::move(layArr), nullptr, - MutableTrackingVolumeVector{}, "Telescope"); - - return std::make_shared(trackVolume); - } - - RotationMatrix3 rotation = RotationMatrix3::Identity(); - std::shared_ptr rBounds = nullptr; - std::shared_ptr surfaceMaterial = nullptr; - - std::vector> detectorStore; - - std::reference_wrapper geoContext; -}; - -// Construct a straight-line propagator. -StraightPropagator makeStraightPropagator( - std::shared_ptr geo) { - Navigator::Config cfg{std::move(geo)}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Navigator navigator(cfg); - StraightLineStepper stepper; - return StraightPropagator(stepper, std::move(navigator)); -} - -// Construct a propagator using a constant magnetic field along z. -ConstantFieldPropagator makeConstantFieldPropagator( - std::shared_ptr geo, double bz, - std::unique_ptr logger) { - Navigator::Config cfg{std::move(geo)}; - cfg.resolvePassive = false; - cfg.resolveMaterial = true; - cfg.resolveSensitive = true; - Navigator navigator(cfg, logger->cloneWithSuffix("Nav")); - auto field = std::make_shared(Vector3(0.0, 0.0, bz)); - ConstantFieldStepper stepper(std::move(field)); - return ConstantFieldPropagator(std::move(stepper), std::move(navigator), - logger->cloneWithSuffix("Prop")); -} - -// Construct initial track parameters. -BoundTrackParameters makeParameters() { - // create covariance matrix from reasonable standard deviations - BoundVector stddev; - stddev[eBoundLoc0] = 100_um; - stddev[eBoundLoc1] = 100_um; - stddev[eBoundTime] = 25_ns; - stddev[eBoundPhi] = 0.5_degree; - stddev[eBoundTheta] = 0.5_degree; - stddev[eBoundQOverP] = 1 / 100_GeV; - BoundMatrix cov = stddev.cwiseProduct(stddev).asDiagonal(); - - auto loc0 = 0. + stddev[eBoundLoc0] * normalDist(rng); - auto loc1 = 0. + stddev[eBoundLoc1] * normalDist(rng); - auto t = 42_ns + stddev[eBoundTime] * normalDist(rng); - auto phi = 0_degree + stddev[eBoundPhi] * normalDist(rng); - auto theta = 90_degree + stddev[eBoundTheta] * normalDist(rng); - auto qOverP = 1_e / 1_GeV + stddev[eBoundQOverP] * normalDist(rng); - - // define a track in the transverse plane along x - Vector4 mPos4(-1_m, loc0, loc1, t); - - return BoundTrackParameters::createCurvilinear(mPos4, phi, theta, qOverP, cov, - ParticleHypothesis::pion()); -} - -// detector resolutions -const MeasurementResolution resPixel = {MeasurementType::eLoc01, - {30_um, 50_um}}; -const MeasurementResolutionMap resolutions = { - {GeometryIdentifier(), resPixel}, -}; - -struct KalmanFitterInputTrajectory { - // The source links - std::vector sourceLinks; - // The start parameters - std::optional startParameters; -}; - -/// -/// Function to create trajectories for kalman fitter -/// -std::vector createTrajectories( - std::shared_ptr geo, std::size_t nTrajectories) { - // simulation propagator - const auto simPropagator = makeStraightPropagator(std::move(geo)); - - std::vector trajectories; - trajectories.reserve(nTrajectories); - - for (unsigned int iTrack = 0; iTrack < nTrajectories; iTrack++) { - auto start = makeParameters(); - // Launch and collect - the measurements - auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start, - resolutions, rng); - - // Extract measurements from result of propagation. - KalmanFitterInputTrajectory traj; - traj.startParameters = start; - traj.sourceLinks = measurements.sourceLinks; - - trajectories.push_back(std::move(traj)); - } - return trajectories; -} -} // namespace - +using namespace Acts::UnitConstants; /// /// @brief Unit test for KF-based alignment algorithm /// BOOST_AUTO_TEST_CASE(ZeroFieldKalmanAlignment) { + aliTestUtils utils; // Build detector - TelescopeDetector detector(geoCtx); + TelescopeDetector detector(utils.geoCtx); const auto geometry = detector(); // reconstruction propagator and fitter @@ -284,28 +39,30 @@ BOOST_AUTO_TEST_CASE(ZeroFieldKalmanAlignment) { // alignment auto alignLogger = getDefaultLogger("Alignment", Logging::INFO); - const auto alignZero = Alignment(std::move(kfZero), std::move(alignLogger)); + const auto alignZero = + ActsAlignment::Alignment(std::move(kfZero), std::move(alignLogger)); // Create 10 trajectories - const auto& trajectories = createTrajectories(geometry, 10); + const auto& trajectories = createTrajectories(geometry, 10, utils); // Construct the KalmanFitter options - auto extensions = getExtensions(); + auto extensions = getExtensions(utils); TestSourceLink::SurfaceAccessor surfaceAccessor{*geometry}; extensions.surfaceAccessor .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); - KalmanFitterOptions kfOptions(geoCtx, magCtx, calCtx, extensions, - PropagatorPlainOptions(geoCtx, magCtx)); + KalmanFitterOptions kfOptions( + utils.geoCtx, utils.magCtx, utils.calCtx, extensions, + PropagatorPlainOptions(utils.geoCtx, utils.magCtx)); // Construct a non-updating alignment updater - AlignedTransformUpdater voidAlignUpdater = + ActsAlignment::AlignedTransformUpdater voidAlignUpdater = [](SurfacePlacementBase* /*element*/, const GeometryContext& /*gctx*/, const Transform3& /*transform*/) { return true; }; // Construct the alignment options - AlignmentOptions> alignOptions( - kfOptions, voidAlignUpdater); + ActsAlignment::AlignmentOptions> + alignOptions(kfOptions, voidAlignUpdater); alignOptions.maxIterations = 1; // Set the surfaces to be aligned (fix the layer 8) @@ -327,7 +84,7 @@ BOOST_AUTO_TEST_CASE(ZeroFieldKalmanAlignment) { auto evaluateRes = alignZero.evaluateTrackAlignmentState( kfOptions.geoContext, inputTraj.sourceLinks, *inputTraj.startParameters, - kfOptions, idxedAlignSurfaces, AlignmentMask::All); + kfOptions, idxedAlignSurfaces, ActsAlignment::AlignmentMask::All); BOOST_CHECK(evaluateRes.ok()); const auto& alignState = evaluateRes.value(); diff --git a/Tests/UnitTests/Plugins/Mille/ActsKalmanToMilleConsistency.cpp b/Tests/UnitTests/Plugins/Mille/ActsKalmanToMilleConsistency.cpp new file mode 100644 index 00000000000..e1745e3347c --- /dev/null +++ b/Tests/UnitTests/Plugins/Mille/ActsKalmanToMilleConsistency.cpp @@ -0,0 +1,160 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/EventData/detail/TestSourceLink.hpp" +#include "Acts/Geometry/CuboidVolumeBuilder.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "ActsAlignment/Kernel/Alignment.hpp" +#include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp" +#include "ActsPlugins/Mille/ActsToMille.hpp" +#include "ActsPlugins/Mille/Helpers.hpp" +#include "ActsTests/CommonHelpers/AlignmentHelpers.hpp" + +#include + +#include +#include + +using namespace Acts; +using namespace ActsTests; +using namespace ActsTests::AlignmentUtils; +using namespace Acts::detail::Test; +using namespace Acts::UnitConstants; + +BOOST_AUTO_TEST_CASE(ZeroFieldKalmanToMille) { + /// Part 1): Build some test data. + /// This is shamelessly "borrowed" from + /// the existing alignment unit test. + + aliTestUtils utils; + // Build detector + TelescopeDetector detector(utils.geoCtx); + const auto geometry = detector(); + + // reconstruction propagator and fitter + auto kfLogger = getDefaultLogger("KalmanFilter", Logging::INFO); + const auto kfZeroPropagator = + makeConstantFieldPropagator(geometry, 0_T, std::move(kfLogger)); + auto kfZero = KalmanFitterType(kfZeroPropagator); + + // alignment + auto alignLogger = getDefaultLogger("Alignment", Logging::INFO); + const auto alignZero = + ActsAlignment::Alignment(std::move(kfZero), std::move(alignLogger)); + + // Create 10 trajectories + const auto& trajectories = createTrajectories(geometry, 10, utils); + + // Construct the KalmanFitter options + + auto extensions = getExtensions(utils); + TestSourceLink::SurfaceAccessor surfaceAccessor{*geometry}; + extensions.surfaceAccessor + .connect<&TestSourceLink::SurfaceAccessor::operator()>(&surfaceAccessor); + KalmanFitterOptions kfOptions( + utils.geoCtx, utils.magCtx, utils.calCtx, extensions, + PropagatorPlainOptions(utils.geoCtx, utils.magCtx)); + + // Construct a non-updating alignment updater + ActsAlignment::AlignedTransformUpdater voidAlignUpdater = + [](SurfacePlacementBase* /*element*/, const GeometryContext& /*gctx*/, + const Transform3& /*transform*/) { return true; }; + + // Construct the alignment options + ActsAlignment::AlignmentOptions> + alignOptions(kfOptions, voidAlignUpdater); + alignOptions.maxIterations = 1; + + // Set the surfaces to be aligned (fix the layer 8) + unsigned int iSurface = 0; + std::unordered_map idxedAlignSurfaces; + // Loop over the detector elements + for (auto& det : detector.detectorStore) { + const auto& surface = det->surface(); + if (surface.geometryId().layer() != 8) { + alignOptions.alignedDetElements.push_back(det.get()); + idxedAlignSurfaces.emplace(&surface, iSurface); + iSurface++; + } + } + + // Test the method to evaluate alignment state for a single track + const auto& inputTraj = trajectories.front(); + kfOptions.referenceSurface = &(*inputTraj.startParameters).referenceSurface(); + + auto evaluateRes = alignZero.evaluateTrackAlignmentState( + kfOptions.geoContext, inputTraj.sourceLinks, *inputTraj.startParameters, + kfOptions, idxedAlignSurfaces, ActsAlignment::AlignmentMask::All); + BOOST_CHECK(evaluateRes.ok()); + + /// Part 2: Now dump the alignment state to Mille. + + auto milleRecord = Mille::spawnMilleRecord("myRecord.root", true); + + const auto& alignState = evaluateRes.value(); + ActsPlugins::ActsToMille::dumpToMille(alignState, milleRecord.get()); + + // trigger file close by destroying the Mille record + milleRecord->flushOutputFile(); + milleRecord.reset(); + + /// Part 3: We read back the record from the binary file, + /// convert it to an alignment state, and compare to + /// what we started from + + // read back the record + auto milleReader = Mille::spawnMilleReader("myRecord.root"); + BOOST_CHECK(milleReader->open("myRecord.root")); + + // Decode it back into a TrackAlignmentState + + ActsAlignment::detail::TrackAlignmentState millePedeState; + // we need to externally supply the alignment parameter indexing logic + millePedeState.alignedSurfaces = alignState.alignedSurfaces; + BOOST_CHECK(ActsPlugins::ActsToMille::unpackMilleRecord( + *milleReader, millePedeState, idxedAlignSurfaces) == + Mille::MilleDecoder::ReadResult::OK); + + // now compare the results! + + BOOST_CHECK_EQUAL(millePedeState.alignmentDof, alignState.alignmentDof); + BOOST_CHECK_EQUAL(millePedeState.measurementDim, alignState.measurementDim); + BOOST_CHECK_EQUAL(millePedeState.trackParametersDim, + alignState.trackParametersDim); + + BOOST_CHECK_EQUAL(millePedeState.residual, alignState.residual); + BOOST_CHECK_EQUAL(millePedeState.chi2, alignState.chi2); + BOOST_CHECK_EQUAL(millePedeState.measurementCovariance, + alignState.measurementCovariance); + BOOST_CHECK_EQUAL(millePedeState.projectionMatrix, + alignState.projectionMatrix); + + BOOST_CHECK_EQUAL(millePedeState.alignmentToResidualDerivative, + alignState.alignmentToResidualDerivative); + + // for the track parameter covariance, we only compare the regularised version + // (differences expected in poorly constrained parameters. Also increase the + // tolerance a bit. + BOOST_CHECK(millePedeState.trackParametersCovariance.isApprox( + ActsPlugins::ActsToMille::regulariseCovariance( + alignState.trackParametersCovariance), + 1.e-6)); + + // for anything derived from the TP covariance, we can also only ask for + // approximate equivalence as a result. + BOOST_CHECK(millePedeState.residualCovariance.isApprox( + alignState.residualCovariance, 1.e-6)); + BOOST_CHECK(millePedeState.alignmentToChi2Derivative.isApprox( + alignState.alignmentToChi2Derivative, 1.e-6)); + BOOST_CHECK(millePedeState.alignmentToChi2SecondDerivative.isApprox( + alignState.alignmentToChi2SecondDerivative, 1.e-6)); +} diff --git a/Tests/UnitTests/Plugins/Mille/ActsMilleBasicCalls.cpp b/Tests/UnitTests/Plugins/Mille/ActsMilleBasicCalls.cpp index 1a123505482..b2dc6a8b1be 100644 --- a/Tests/UnitTests/Plugins/Mille/ActsMilleBasicCalls.cpp +++ b/Tests/UnitTests/Plugins/Mille/ActsMilleBasicCalls.cpp @@ -23,9 +23,9 @@ BOOST_AUTO_TEST_CASE(OpenMilleRecord) { std::unique_ptr theRecord = Mille::spawnMilleRecord(fname); BOOST_CHECK_NE(theRecord.get(), nullptr); - // for now, will write one entry of dummy information. - ActsAlignment::detail::TrackAlignmentState dummyState; - dumpToMille(dummyState, theRecord.get()); + theRecord->addData(0.02, 0.005, std::vector{0u}, + std::vector{1.0}, std::vector{42}, + std::vector{-1.0}); theRecord->writeRecord(); theRecord->flushOutputFile(); theRecord.reset(nullptr); diff --git a/Tests/UnitTests/Plugins/Mille/CMakeLists.txt b/Tests/UnitTests/Plugins/Mille/CMakeLists.txt index c55ea456a35..eaebe537293 100644 --- a/Tests/UnitTests/Plugins/Mille/CMakeLists.txt +++ b/Tests/UnitTests/Plugins/Mille/CMakeLists.txt @@ -1,3 +1,4 @@ set(unittest_extra_libraries ActsPluginMille) +add_unittest(ActsKalmanToMilleConsistency ActsKalmanToMilleConsistency.cpp) add_unittest(ActsMilleBasicCalls ActsMilleBasicCalls.cpp)