Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-07 08:09:57

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