Add alignment workflow for the OpenDataDetector

.zenodo.json

@@ -1,8 +1,8 @@
 {
   "description": "Experiment-independent toolkit for (charged) particle track reconstruction in (high energy) physics experiments implemented in modern C++ ",
   "license": "MPL-2.0",
-  "title": "acts-project/acts: v11.0.0",
-  "version": "v11.0.0",
+  "title": "acts-project/acts: v44.0.1",
+  "version": "v44.0.1",
   "upload_type": "software",
   "creators": [
     {
@@ -130,17 +130,6 @@
       "scheme": "url",
       "identifier": "https://github.com/acts-project/acts/tree/v10.0.0",
       "relation": "isSupplementTo"
-    },
-    {
-      "scheme": "doi",
-      "identifier": "10.5281/zenodo.4818137",
-      "relation": "obsoletes",
-      "resource_type": "software"
-    },
-    {
-      "scheme": "doi",
-      "identifier": "10.5281/zenodo.5141418",
-      "relation": "isVersionOf"
     }
   ]
 }

Alignment/include/ActsAlignment/Kernel/Alignment.ipp

@@ -13,6 +13,13 @@
 #include "Acts/EventData/VectorMultiTrajectory.hpp"
 #include "Acts/EventData/VectorTrackContainer.hpp"
 
+#include <algorithm>
+#include <iomanip>
+#include <stdexcept>
+#include <type_traits>
+#include <typeinfo>
+#include <unordered_map>
+
 template <typename fitter_t>
 template <typename source_link_t, typename start_parameters_t,
           typename fit_options_t>
@@ -31,23 +38,48 @@ ActsAlignment::Alignment<fitter_t>::evaluateTrackAlignmentState(
   Acts::SourceLinkAdapterIterator begin{sourceLinks.begin()};
   Acts::SourceLinkAdapterIterator end{sourceLinks.end()};
 
+  // Validate input data before fitting
+  if (sourceLinks.empty()) {
+    ACTS_WARNING("Empty source links for track");
+    return ActsAlignment::AlignmentError::NoAlignmentDofOnTrack;
+  }
+
   // Perform the fit
+  ACTS_DEBUG("Fitting track with " << sourceLinks.size() << " measurements");
   auto fitRes = m_fitter.fit(begin, end, sParameters, fitOptions, tracks);
 
   if (!fitRes.ok()) {
-    ACTS_WARNING("Fit failure");
+    ACTS_WARNING("Fit failure for track: " << fitRes.error());
     return fitRes.error();
   }
+
   // The fit results
   const auto& track = fitRes.value();
+
+  // Check if track has valid tip index
+  if (track.tipIndex() == Acts::MultiTrajectoryTraits::kInvalid) {
+    ACTS_WARNING("Invalid tip index for fitted track");
+    return ActsAlignment::AlignmentError::NoAlignmentDofOnTrack;
+  }
+
   // Calculate the global track parameters covariance with the fitted track
-  const auto& globalTrackParamsCov =
+  const auto globalTrackParamsCov =
       Acts::detail::globalTrackParametersCovariance(
           tracks.trackStateContainer(), track.tipIndex());
+
+  // Check if the result is valid
+  if (globalTrackParamsCov.first.rows() == 0 ||
+      globalTrackParamsCov.first.cols() == 0) {
+    ACTS_WARNING("Empty global track parameters covariance for track "
+                 << track.index());
+    return ActsAlignment::AlignmentError::NoAlignmentDofOnTrack;
+  }
+
   // Calculate the alignment state
   const auto alignState = detail::trackAlignmentState(
       gctx, tracks.trackStateContainer(), track.tipIndex(),
       globalTrackParamsCov, idxedAlignSurfaces, alignMask);
+
   if (alignState.alignmentDof == 0) {
     ACTS_VERBOSE("No alignment dof on track!");
     return AlignmentError::NoAlignmentDofOnTrack;
@@ -77,6 +109,7 @@ void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
   Acts::ActsDynamicMatrix sumChi2SecondDerivative =
       Acts::ActsDynamicMatrix::Zero(alignResult.alignmentDof,
                                     alignResult.alignmentDof);
+
   // Copy the fit options
   fit_options_t fitOptionsWithRefSurface = fitOptions;
   // Calculate contribution to chi2 derivatives from all input trajectories
@@ -85,6 +118,7 @@ void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
   alignResult.measurementDim = 0;
   alignResult.numTracks = trajectoryCollection.size();
   double sumChi2ONdf = 0;
+
   for (unsigned int iTraj = 0; iTraj < trajectoryCollection.size(); iTraj++) {
     const auto& sourceLinks = trajectoryCollection.at(iTraj);
     const auto& sParameters = startParametersCollection.at(iTraj);
@@ -146,7 +180,8 @@ void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
       -sumChi2SecondDerivative.fullPivLu().solve(sumChi2Derivative);
   ACTS_VERBOSE("sumChi2SecondDerivative = \n" << sumChi2SecondDerivative);
   ACTS_VERBOSE("sumChi2Derivative = \n" << sumChi2Derivative);
-  ACTS_VERBOSE("alignResult.deltaAlignmentParameters \n");
+  ACTS_VERBOSE("deltaAlignmentParameters = \n"
+               << alignResult.deltaAlignmentParameters);
 
   // Alignment parameters covariance
   alignResult.alignmentCovariance = 2 * sumChi2SecondDerivativeInverse;
@@ -180,7 +215,11 @@ ActsAlignment::Alignment<fitter_t>::updateAlignmentParameters(
         deltaAlignmentParam.segment<3>(Acts::eAlignmentRotation0);
 
     // 3. The new transform
-    const Acts::Vector3 newCenter = oldCenter + deltaCenter;
+    const auto& rotation = oldTransform.rotation();
+    const Acts::Vector3 newCenter =
+        oldCenter + deltaCenter(0) * rotation.col(0) +
+        deltaCenter(1) * rotation.col(1) + deltaCenter(2) * rotation.col(2);
+
     Acts::Transform3 newTransform = oldTransform;
     newTransform.translation() = newCenter;
     // Rotation first around fixed local x, then around fixed local y, and last
@@ -192,7 +231,6 @@ ActsAlignment::Alignment<fitter_t>::updateAlignmentParameters(
         Acts::AngleAxis3(deltaEulerAngles(1), Acts::Vector3::UnitY());
     newTransform *=
         Acts::AngleAxis3(deltaEulerAngles(0), Acts::Vector3::UnitX());
-
     // 4. Update the aligned transform
     //@Todo: use a better way to handle this (need dynamic cast to inherited
     // detector element type)
@@ -245,10 +283,18 @@ ActsAlignment::Alignment<fitter_t>::align(
     if (iter_it != alignOptions.iterationState.end()) {
       alignMask = iter_it->second;
     }
+
     // Calculate the alignment parameters delta etc.
     calculateAlignmentParameters(
         trajectoryCollection, startParametersCollection,
         alignOptions.fitOptions, alignResult, alignMask);
+
+    // Check if calculation failed (e.g., no successful tracks)
+    if (!alignResult.result.ok()) {
+      ACTS_ERROR("Alignment calculation failed in iteration " << iIter);
+      return alignResult.result.error();
+    }
+
     // Screen out the information
     ACTS_INFO("iIter = " << iIter << ", total chi2 = " << alignResult.chi2
                          << ", total measurementDim = "
@@ -286,7 +332,10 @@ ActsAlignment::Alignment<fitter_t>::align(
 
     ACTS_INFO("The solved delta of alignmentParameters = \n "
               << alignResult.deltaAlignmentParameters);
-    // Not coveraged yet, update the detector element alignment parameters
+
+    // ========================================================================
+    // UPDATE: Apply the alignment corrections to the detector elements
+    // ========================================================================
     auto updateRes = updateAlignmentParameters(
         alignOptions.fitOptions.geoContext, alignOptions.alignedDetElements,
         alignOptions.alignedTransformUpdater, alignResult);
@@ -295,6 +344,9 @@ ActsAlignment::Alignment<fitter_t>::align(
       return updateRes.error();
     }
     alignmentParametersUpdated = true;
+    ACTS_VERBOSE("Alignment parameters updated for iteration " << iIter);
+    // ========================================================================
+
   }  // end of all iterations
 
   // Alignment failure if not converged

CITATION.cff

@@ -92,7 +92,7 @@ authors:
 - given-names: Carlo
   family-names: Varni
   affiliation: UC Berkeley
-version: 10.0.0
+version: v44.0.1
 date-released: 2021-07-28
 repository-code: https://github.com/acts-project/acts
 license: MPL-2.0

Core/src/Surfaces/Surface.cpp

@@ -92,18 +92,21 @@ AlignmentToBoundMatrix Surface::alignmentToBoundDerivativeWithoutCorrection(
   // The local frame rotation
   const auto& rotation = transform(gctx).rotation();
   // The axes of local frame
-  const auto& localXAxis = rotation.col(0);
-  const auto& localYAxis = rotation.col(1);
-  const auto& localZAxis = rotation.col(2);
+  // const auto& localXAxis = rotation.col(0);
+  // const auto& localYAxis = rotation.col(1);
+  // const auto& localZAxis = rotation.col(2);
   // Calculate the derivative of local frame axes w.r.t its rotation
   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
       detail::rotationToLocalAxesDerivative(rotation);
   // Calculate the derivative of local 3D Cartesian coordinates w.r.t.
   // alignment parameters (without path correction)
   AlignmentToPositionMatrix alignToLoc3D = AlignmentToPositionMatrix::Zero();
-  alignToLoc3D.block<1, 3>(eX, eAlignmentCenter0) = -localXAxis.transpose();
-  alignToLoc3D.block<1, 3>(eY, eAlignmentCenter0) = -localYAxis.transpose();
-  alignToLoc3D.block<1, 3>(eZ, eAlignmentCenter0) = -localZAxis.transpose();
+  alignToLoc3D.block<1, 3>(eX, eAlignmentCenter0) =
+      -Vector3::UnitX().transpose();  // [-1, 0, 0]
+  alignToLoc3D.block<1, 3>(eY, eAlignmentCenter0) =
+      -Vector3::UnitY().transpose();  // [0, -1, 0]
+  alignToLoc3D.block<1, 3>(eZ, eAlignmentCenter0) =
+      -Vector3::UnitZ().transpose();  // [0, 0, 1]
   alignToLoc3D.block<1, 3>(eX, eAlignmentRotation0) =
       pcRowVec * rotToLocalXAxis;
   alignToLoc3D.block<1, 3>(eY, eAlignmentRotation0) =
@@ -141,7 +144,7 @@ AlignmentToPathMatrix Surface::alignmentToPathDerivative(
   // Initialize the derivative of propagation path w.r.t. local frame
   // translation (origin) and rotation
   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
-  alignToPath.segment<3>(eAlignmentCenter0) = localZAxis.transpose() / dz;
+  alignToPath.segment<3>(eAlignmentCenter0) = Vector3::UnitZ().transpose() / dz;
   alignToPath.segment<3>(eAlignmentRotation0) =
       -pcRowVec * rotToLocalZAxis / dz;
 

Examples/Algorithms/Alignment/include/ActsExamples/Alignment/AlignmentAlgorithm.hpp

@@ -82,7 +82,8 @@ class AlignmentAlgorithm final : public IAlgorithm {
   /// contains shared_ptr anyway.
   static std::shared_ptr<AlignmentFunction> makeAlignmentFunction(
       std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
-      std::shared_ptr<const Acts::MagneticFieldProvider> magneticField);
+      std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
+      Acts::Logging::Level logLevel = Acts::Logging::INFO);
 
   struct Config {
     /// Input measurements collection.
@@ -93,6 +94,8 @@ class AlignmentAlgorithm final : public IAlgorithm {
     std::string inputInitialTrackParameters;
     /// Output aligned parameters collection.
     std::string outputAlignmentParameters;
+    /// Tracking geometry for surface access
+    std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
     /// Type erased fitter function.
     std::shared_ptr<AlignmentFunction> align;
     /// The aligned transform updater
@@ -119,12 +122,19 @@ class AlignmentAlgorithm final : public IAlgorithm {
   AlignmentAlgorithm(Config cfg, Acts::Logging::Level lvl);
 
   /// Framework execute method of the alignment algorithm
+  /// Collects track data from each event
   ///
   /// @param ctx is the algorithm context that holds event-wise information
   /// @return a process code to steer the algorithm flow
   ActsExamples::ProcessCode execute(
       const ActsExamples::AlgorithmContext& ctx) const override;
 
+  /// Framework finalize method
+  /// Performs the actual alignment using all collected tracks
+  ///
+  /// @return a process code to steer the algorithm flow
+  ProcessCode finalize() override;
+
  private:
   Config m_cfg;
 
@@ -136,6 +146,14 @@ class AlignmentAlgorithm final : public IAlgorithm {
                                                          "InputProtoTracks"};
   WriteDataHandle<AlignmentParameters> m_outputAlignmentParameters{
       this, "OutputAlignmentParameters"};
+
+  // Collected track data from all events
+  mutable std::vector<std::vector<IndexSourceLink>> m_collectedSourceLinks;
+  mutable TrackParametersContainer m_collectedInitialParameters;
+  mutable std::shared_ptr<MeasurementContainer> m_collectedMeasurements;
+  mutable Acts::GeometryContext m_savedGeoContext;
+  mutable Acts::MagneticFieldContext m_savedMagFieldContext;
+  mutable Acts::CalibrationContext m_savedCalibContext;
 };
 
 }  // namespace ActsExamples