feat: Demonstrator for MillePede alignment of ACTS Kalman tracks

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<detail::TrackAlignmentState>& trackAlignmentStates,
+      AlignmentResult& alignResult) const;
+
   /// @brief update the detector element alignment parameters
   ///
   /// @param gctx The geometry context

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 <queue>
 
@@ -76,19 +77,14 @@ void ActsAlignment::Alignment<fitter_t>::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
   // @Todo: How to update the source link error iteratively?
   alignResult.chi2 = 0;
   alignResult.measurementDim = 0;
   alignResult.numTracks = trajectoryCollection.size();
-  double sumChi2ONdf = 0;
+  std::vector<detail::TrackAlignmentState> 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<fitter_t>::calculateAlignmentParameters(
       continue;
     }
     const auto& alignState = evaluateRes.value();
+    alignmentStates.push_back(alignState);
+  }
+  return calculateAlignmentParameters(alignmentStates, alignResult);
+}
+
+template <typename fitter_t>
+void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
+    const std::vector<detail::TrackAlignmentState>& 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

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;
 }

Examples/Algorithms/AlignmentMillePede/CMakeLists.txt

@@ -0,0 +1,17 @@
+acts_add_library(
+    ExamplesAlignmentMillePede
+    src/MillePedeAlignmentSandbox.cpp
+    src/ActsSolverFromMille.cpp
+)
+target_include_directories(
+    ActsExamplesAlignmentMillePede
+    PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+)
+target_link_libraries(
+    ActsExamplesAlignmentMillePede
+    PUBLIC
+        Acts::Alignment
+        Acts::ExamplesFramework
+        Acts::ExamplesMagneticField
+        ActsPluginMille
+)

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 <memory>
+
+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<Stepper, Acts::Navigator>;
+  using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>;
+  using Alignment = ActsAlignment::Alignment<Fitter>;
+
+  using AlignmentParameters =
+      std::unordered_map<Acts::SurfacePlacementBase*, Acts::Transform3>;
+
+  /// 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<const Acts::TrackingGeometry> trackingGeometry;
+    // magnetic field
+    std::shared_ptr<const Acts::MagneticFieldProvider> magneticField;
+    // modules to fix in the alignment to suppress global movements
+    std::set<Acts::GeometryIdentifier> 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<const Acts::Logger> 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<Alignment> m_align;
+  /// tracking geometry
+  std::shared_ptr<const Acts::TrackingGeometry> m_trackingGeometry;
+};
+
+}  // namespace ActsExamples

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 <memory>
+
+#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<Stepper, Acts::Navigator>;
+  using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>;
+  using Alignment = ActsAlignment::Alignment<Fitter>;
+
+  using AlignmentParameters =
+      std::unordered_map<Acts::SurfacePlacementBase*, Acts::Transform3>;
+
+  /// 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<const Acts::TrackingGeometry> trackingGeometry;
+    // magnetic field
+    std::shared_ptr<const Acts::MagneticFieldProvider> magneticField;
+    // modules to fix in the alignment to suppress global movements
+    std::set<Acts::GeometryIdentifier> 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<const Acts::Logger> 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<MeasurementContainer> m_inputMeasurements{this,
+                                                           "InputMeasurements"};
+  /// tracks to use for the alignment
+  ReadDataHandle<ConstTrackContainer> m_inputTracks{this, "InputTracks"};
+
+  /// alignment module instance - reuse as much as possible
+  std::shared_ptr<Alignment> m_align;
+  /// tracking geometry
+  std::shared_ptr<const Acts::TrackingGeometry> m_trackingGeometry;
+  /// the Mille record instance for writing our alignment info.
+  std::unique_ptr<Mille::MilleRecord> m_milleOut = nullptr;
+};
+
+}  // namespace ActsExamples