Kernel/detail/AlignmentEngine.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Alignment.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/MultiTrajectory.hpp"
0014 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/Logger.hpp"
0018 #include "ActsAlignment/Kernel/AlignmentMask.hpp"
0019
0020 #include <unordered_map>
0021 #include <vector>
0022
0023 namespace ActsAlignment {
0024 struct AlignmentResult;
0025 }  // namespace ActsAlignment
0026
0027 namespace ActsAlignment::detail {
0028
0029 ///
0030 ///@brief struct to store info needed for track-based alignment
0031 ///
0032 struct TrackAlignmentState {
0033   // The dimension of measurements
0034   std::size_t measurementDim = 0;
0035
0036   // The dimension of track parameters
0037   std::size_t trackParametersDim = 0;
0038
0039   // The contributed alignment degree of freedom
0040   std::size_t alignmentDof = 0;
0041
0042   // The measurements covariance
0043   Acts::DynamicMatrix measurementCovariance;
0044
0045   // The track parameters covariance
0046   Acts::DynamicMatrix trackParametersCovariance;
0047
0048   // The projection matrix
0049   Acts::DynamicMatrix projectionMatrix;
0050
0051   // The residual
0052   Acts::DynamicVector residual;
0053
0054   // The covariance of residual
0055   Acts::DynamicMatrix residualCovariance;
0056
0057   // The chi2
0058   double chi2 = 0;
0059
0060   // The derivative of residual w.r.t. alignment parameters
0061   Acts::DynamicMatrix alignmentToResidualDerivative;
0062
0063   // The derivative of chi2 w.r.t. alignment parameters
0064   Acts::DynamicVector alignmentToChi2Derivative;
0065
0066   // The second derivative of chi2 w.r.t. alignment parameters
0067   Acts::DynamicMatrix alignmentToChi2SecondDerivative;
0068
0069   // The alignable surfaces on the track and their indices in both the global
0070   // alignable surfaces pool and those relevant with this track
0071   std::unordered_map<const Acts::Surface*, std::pair<std::size_t, std::size_t>>
0072       alignedSurfaces;
0073 };
0074
0075 /// Reset some columns of the alignment to bound derivative to zero if the
0076 /// relevant degree of freedom is fixed
0077 ///
0078 /// @param alignToBound The alignment to bound parameters derivative
0079 /// @param mask The alignment mask
0080 void resetAlignmentDerivative(Acts::AlignmentToBoundMatrix& alignToBound,
0081                               AlignmentMask mask);
0082
0083 /// @brief Helper function to calculate first and second derivative
0084 /// of chi2 w.r.t. alignment parameters for a single track once the
0085 /// relevant information has
0086 /// been retrieved from the respective track.
0087 /// Updates in-place the chi2 and the matrix/vector forming the final
0088 /// equation.
0089 /// @param [in,out] alignState: TrackAlignmentState to modify (in-place).
0090 ///
0091 /// @note This operates purely on the (dynamically sized) matrices already
0092 /// stored in @p alignState and does not depend on the track/trajectory type.
0093 /// It is deliberately defined out-of-line (in AlignmentEngine.cpp) so its
0094 /// Eigen expression templates are instantiated once in the Acts alignment
0095 /// library instead of in every translation unit that instantiates
0096 /// @c trackAlignmentState.
0097 void finaliseTrackAlignState(TrackAlignmentState& alignState);
0098
0099 /// @brief Assemble the summed chi2 derivatives from the per-track alignment
0100 /// states and solve for the alignment parameter delta.
0101 ///
0102 /// @param trackAlignmentStates the per-track alignment states
0103 /// @param [in,out] alignResult the alignment result to fill
0104 /// @param logger the logger to use
0105 ///
0106 /// @note This operates purely on the (dynamically sized) matrices in the
0107 /// inputs and does not depend on the fitter/track type. It is defined
0108 /// out-of-line (in AlignmentEngine.cpp) so its Eigen algebra (matrix inverse,
0109 /// full-pivot LU solve) is instantiated once in the Acts alignment library
0110 /// rather than for every @c Alignment<fitter_t> specialization.
0111 void solveAlignmentParameters(
0112     const std::vector<TrackAlignmentState>& trackAlignmentStates,
0113     AlignmentResult& alignResult, const Acts::Logger& logger);
0114
0115 ///
0116 /// Calculate the first and second derivative of chi2 w.r.t. alignment
0117 /// parameters for a single track
0118 ///
0119 /// Suppose there are n measurements on the track, and m (m<=n) of them are on
0120 /// alignable surface, then (eAlignmentSize*m) alignment parameters
0121 /// will be involved for this particular track, i.e. this track will contribute
0122 /// to at most (eAlignmentSize*m*2) elements of the full chi2
0123 /// second derivative matrix
0124 ///
0125 /// @tparam source_link_t The source link type of the trajectory
0126 ///
0127 /// @param gctx The current geometry context object
0128 /// @param multiTraj The MultiTrajectory containing the trajectory to be
0129 /// investigated
0130 /// @param entryIndex The trajectory entry index
0131 /// @param globalTrackParamsCov The global track parameters covariance for a
0132 /// single track and the starting row/column for smoothed states. This contains
0133 /// all smoothed track states including those non-measurement states. Selection
0134 /// of certain rows/columns for measurement states is needed.
0135 /// @param idxedAlignSurfaces The indexed surfaces to be aligned
0136 ///
0137 /// @return The track alignment state containing fundamental alignment
0138 /// ingredients
0139 template <typename traj_t>
0140 TrackAlignmentState trackAlignmentState(
0141     const Acts::GeometryContext& gctx, const traj_t& multiTraj,
0142     Acts::TrackIndexType entryIndex,
0143     const std::pair<Acts::DynamicMatrix,
0144                     std::unordered_map<std::size_t, std::size_t>>&
0145         globalTrackParamsCov,
0146     const std::unordered_map<const Acts::Surface*, std::size_t>&
0147         idxedAlignSurfaces,
0148     const AlignmentMask& alignMask) {
0149   // Construct an alignment state
0150   TrackAlignmentState alignState;
0151
0152   // Remember the index within the trajectory and whether it's alignable
0153   std::vector<std::pair<Acts::TrackIndexType, bool>> measurementStates;
0154   measurementStates.reserve(15);
0155   // Number of smoothed states on the track
0156   // std::size_t nSmoothedStates = 0; // commented because clang-tidy complains
0157   // about unused Number of alignable surfaces on the track
0158   std::size_t nAlignSurfaces = 0;
0159
0160   // Visit the track states on the track
0161   multiTraj.visitBackwards(entryIndex, [&](const auto& ts) {
0162     // Remember the number of smoothed states
0163     if (ts.hasSmoothed()) {
0164       // nSmoothedStates++; // commented because clang-tidy complains about
0165       // unused
0166     } else {
0167       // @note: this should in principle never happen now. But still keep it as a note
0168       return true;
0169     }
0170
0171     // Only measurement states matter (we can't align non-measurement states,
0172     // no?)
0173     if (!ts.typeFlags().isMeasurement()) {
0174       return true;
0175     }
0176     // Check if the reference surface is to be aligned
0177     bool isAlignable = false;
0178     const auto surface = &ts.referenceSurface();
0179     if (auto it = idxedAlignSurfaces.find(surface);
0180         it != idxedAlignSurfaces.end()) {
0181       isAlignable = true;
0182       // Remember the surface and its index
0183       alignState.alignedSurfaces[surface].first = it->second;
0184       nAlignSurfaces++;
0185     }
0186     // Remember the index of the state within the trajectory and whether it's
0187     // alignable
0188     measurementStates.push_back({ts.index(), isAlignable});
0189     // Add up measurement dimension
0190     alignState.measurementDim += ts.calibratedSize();
0191     return true;
0192   });
0193
0194   // Return now if the track contains no alignable surfaces
0195   if (nAlignSurfaces == 0) {
0196     return alignState;
0197   }
0198
0199   // The alignment degree of freedom
0200   alignState.alignmentDof = Acts::eAlignmentSize * nAlignSurfaces;
0201   // Dimension of global track parameters (from only measurement states)
0202   alignState.trackParametersDim = Acts::eBoundSize * measurementStates.size();
0203
0204   // Initialize the alignment matrices with components from the measurement
0205   // states
0206   // The measurement covariance
0207   alignState.measurementCovariance = Acts::DynamicMatrix::Zero(
0208       alignState.measurementDim, alignState.measurementDim);
0209   // The bound parameters to measurement projection matrix
0210   alignState.projectionMatrix = Acts::DynamicMatrix::Zero(
0211       alignState.measurementDim, alignState.trackParametersDim);
0212   // The derivative of residual w.r.t. alignment parameters
0213   alignState.alignmentToResidualDerivative = Acts::DynamicMatrix::Zero(
0214       alignState.measurementDim, alignState.alignmentDof);
0215   // The track parameters covariance
0216   alignState.trackParametersCovariance = Acts::DynamicMatrix::Zero(
0217       alignState.trackParametersDim, alignState.trackParametersDim);
0218   // The residual
0219   alignState.residual = Acts::DynamicVector::Zero(alignState.measurementDim);
0220
0221   // Unpack global track parameters covariance and the starting row/column for
0222   // all smoothed states.
0223   // Note that the dimension of provided global track parameters covariance
0224   // should be same as eBoundSize * nSmoothedStates
0225   const auto& [sourceTrackParamsCov, stateRowIndices] = globalTrackParamsCov;
0226
0227   // Loop over the measurement states to fill those alignment matrices
0228   // This is done in reverse order
0229   std::size_t iMeasurement = alignState.measurementDim;
0230   std::size_t iParams = alignState.trackParametersDim;
0231   std::size_t iSurface = nAlignSurfaces;
0232   for (const auto& [rowStateIndex, isAlignable] : measurementStates) {
0233     const auto& state = multiTraj.getTrackState(rowStateIndex);
0234     const std::size_t measdim = state.calibratedSize();
0235     // Update index of current measurement and parameter
0236     iMeasurement -= measdim;
0237     iParams -= Acts::eBoundSize;
0238     // (a) Get and fill the measurement covariance matrix
0239     const Acts::DynamicMatrix measCovariance =
0240         state.effectiveCalibratedCovariance();
0241     alignState.measurementCovariance.block(iMeasurement, iMeasurement, measdim,
0242                                            measdim) = measCovariance;
0243
0244     // (b) Get and fill the bound parameters to measurement projection matrix
0245     const Acts::DynamicMatrix H =
0246         state.projectorSubspaceHelper().fullProjector().topLeftCorner(
0247             measdim, Acts::eBoundSize);
0248     alignState.projectionMatrix.block(iMeasurement, iParams, measdim,
0249                                       Acts::eBoundSize) = H;
0250     // (c) Get and fill the residual
0251     alignState.residual.segment(iMeasurement, measdim) =
0252         state.effectiveCalibrated() - H * state.smoothed();
0253
0254     // (d) Get the derivative of alignment parameters w.r.t. measurement
0255     // or residual
0256     if (isAlignable) {
0257       iSurface -= 1;
0258       const auto surface = &state.referenceSurface();
0259       alignState.alignedSurfaces.at(surface).second = iSurface;
0260       // The free parameters transformed from the smoothed parameters
0261       const Acts::FreeVector freeParams =
0262           Acts::MultiTrajectoryHelpers::freeSmoothed(gctx, state);
0263       // The position
0264       const Acts::Vector3 position = freeParams.segment<3>(Acts::eFreePos0);
0265       // The direction
0266       const Acts::Vector3 direction = freeParams.segment<3>(Acts::eFreeDir0);
0267       // The derivative of free parameters w.r.t. path length. @note Here, we
0268       // assume a linear track model, i.e. neglecting the change of track
0269       // direction. Otherwise, we need to know the magnetic field at the free
0270       // parameters
0271       Acts::FreeVector pathDerivative = Acts::FreeVector::Zero();
0272       pathDerivative.head<3>() = direction;
0273       // Get the derivative of bound parameters w.r.t. alignment parameters
0274       Acts::AlignmentToBoundMatrix alignToBound =
0275           surface->alignmentToBoundDerivative(gctx, position, direction,
0276                                               pathDerivative);
0277       // Set the degree of freedom per surface.
0278       // @Todo: don't allocate memory for fixed degree of freedom and consider surface/layer/volume wise align mask (instead of using global mask as now)
0279       resetAlignmentDerivative(alignToBound, alignMask);
0280
0281       // Residual is calculated as the (measurement - parameters), thus we need
0282       // a minus sign below
0283       alignState.alignmentToResidualDerivative.block(
0284           iMeasurement, iSurface * Acts::eAlignmentSize, measdim,
0285           Acts::eAlignmentSize) = -H * alignToBound;
0286     }
0287
0288     // (e) Extract and fill the track parameters covariance matrix for only
0289     // measurement states
0290     // @Todo: add helper function to select rows/columns of a matrix
0291     for (unsigned int iColState = 0; iColState < measurementStates.size();
0292          iColState++) {
0293       std::size_t colStateIndex = measurementStates.at(iColState).first;
0294       // Retrieve the block from the source covariance matrix
0295       Acts::BoundMatrix correlation =
0296           sourceTrackParamsCov.block<Acts::eBoundSize, Acts::eBoundSize>(
0297               stateRowIndices.at(rowStateIndex),
0298               stateRowIndices.at(colStateIndex));
0299       // Fill the block of the target covariance matrix
0300       std::size_t iCol =
0301           alignState.trackParametersDim - (iColState + 1) * Acts::eBoundSize;
0302       alignState.trackParametersCovariance
0303           .block<Acts::eBoundSize, Acts::eBoundSize>(iParams, iCol) =
0304           correlation;
0305     }
0306   }
0307   finaliseTrackAlignState(alignState);
0308
0309   return alignState;
0310 }
0311
0312 }  // namespace ActsAlignment::detail