ActsAlignment/Kernel/Alignment.ipp

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 "ActsAlignment/Kernel/Alignment.hpp"
0012
0013 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0014 #include "Acts/EventData/VectorTrackContainer.hpp"
0015 #include "Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp"
0016 #include "Acts/Utilities/Logger.hpp"
0017 #include "Acts/Utilities/detail/EigenCompat.hpp"
0018 #include "ActsAlignment/Kernel/AlignmentError.hpp"
0019 #include "ActsAlignment/Kernel/detail/AlignmentEngine.hpp"
0020
0021 #include <queue>
0022
0023 template <typename fitter_t>
0024 template <typename source_link_t, typename fit_options_t>
0025 Acts::Result<ActsAlignment::detail::TrackAlignmentState>
0026 ActsAlignment::Alignment<fitter_t>::evaluateTrackAlignmentState(
0027     const Acts::GeometryContext& gctx,
0028     const std::vector<source_link_t>& sourceLinks,
0029     const Acts::BoundTrackParameters& sParameters,
0030     const fit_options_t& fitOptions,
0031     const std::unordered_map<const Acts::Surface*, std::size_t>&
0032         idxedAlignSurfaces,
0033     const ActsAlignment::AlignmentMask& alignMask,
0034     const ActsAlignment::AlignmentHierarchy* hierarchy) const {
0035   Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0036                               Acts::VectorMultiTrajectory{}};
0037
0038   // Convert to Acts::SourceLink during iteration
0039   Acts::SourceLinkAdapterIterator begin{sourceLinks.begin()};
0040   Acts::SourceLinkAdapterIterator end{sourceLinks.end()};
0041
0042   // Perform the fit
0043   auto fitRes = m_fitter.fit(begin, end, sParameters, fitOptions, tracks);
0044
0045   if (!fitRes.ok()) {
0046     ACTS_WARNING("Fit failure: " << fitRes.error().message());
0047     return fitRes.error();
0048   }
0049   // The fit results
0050   const auto& track = fitRes.value();
0051   // Calculate the global track parameters covariance with the fitted track
0052   const auto& globalTrackParamsCov =
0053       Acts::detail::globalTrackParametersCovariance(
0054           tracks.trackStateContainer(), track.tipIndex());
0055   // Calculate the alignment state
0056   const auto alignState = detail::trackAlignmentState(
0057       gctx, tracks.trackStateContainer(), track.tipIndex(),
0058       globalTrackParamsCov, idxedAlignSurfaces, alignMask, hierarchy);
0059   if (alignState.alignmentDof == 0) {
0060     ACTS_VERBOSE("No alignment dof on track!");
0061     return AlignmentError::NoAlignmentDofOnTrack;
0062   }
0063   return alignState;
0064 }
0065
0066 template <typename fitter_t>
0067 template <typename trajectory_container_t,
0068           typename start_parameters_container_t, typename fit_options_t>
0069 void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
0070     const trajectory_container_t& trajectoryCollection,
0071     const start_parameters_container_t& startParametersCollection,
0072     const fit_options_t& fitOptions,
0073     ActsAlignment::AlignmentResult& alignResult,
0074     const ActsAlignment::AlignmentMask& alignMask,
0075     const ActsAlignment::AlignmentHierarchy* hierarchy) const {
0076   // The number of trajectories must be equal to the number of starting
0077   // parameters
0078   assert(trajectoryCollection.size() == startParametersCollection.size());
0079
0080   // The total alignment degree of freedom
0081   alignResult.alignmentDof =
0082       alignResult.idxedAlignSurfaces.size() * Acts::eAlignmentSize;
0083   // Copy the fit options
0084   fit_options_t fitOptionsWithRefSurface = fitOptions;
0085   // Calculate contribution to chi2 derivatives from all input trajectories
0086   // @Todo: How to update the source link error iteratively?
0087   alignResult.chi2 = 0;
0088   alignResult.measurementDim = 0;
0089   alignResult.numTracks = trajectoryCollection.size();
0090   std::vector<detail::TrackAlignmentState> alignmentStates;
0091   for (unsigned int iTraj = 0; iTraj < trajectoryCollection.size(); iTraj++) {
0092     const auto& sourceLinks = trajectoryCollection.at(iTraj);
0093     const auto& sParameters = startParametersCollection.at(iTraj);
0094     // Set the target surface
0095     fitOptionsWithRefSurface.referenceSurface = &sParameters.referenceSurface();
0096     // The result for one single track
0097     auto evaluateRes = evaluateTrackAlignmentState(
0098         fitOptions.geoContext, sourceLinks, sParameters,
0099         fitOptionsWithRefSurface, alignResult.idxedAlignSurfaces, alignMask,
0100         hierarchy);
0101     if (!evaluateRes.ok()) {
0102       ACTS_DEBUG("Evaluation of alignment state for track " << iTraj
0103                                                             << " failed");
0104       continue;
0105     }
0106     const auto& alignState = evaluateRes.value();
0107     alignmentStates.push_back(alignState);
0108   }
0109   return calculateAlignmentParameters(alignmentStates, alignResult);
0110 }
0111
0112 template <typename fitter_t>
0113 void ActsAlignment::Alignment<fitter_t>::calculateAlignmentParameters(
0114     const std::vector<detail::TrackAlignmentState>& trackAlignmentStates,
0115     AlignmentResult& alignResult) const {
0116   // Delegate to the out-of-line, fitter-independent implementation so its Eigen
0117   // algebra is not re-instantiated for every Alignment<fitter_t>.
0118   detail::solveAlignmentParameters(trackAlignmentStates, alignResult, logger());
0119 }
0120
0121 template <typename fitter_t>
0122 double ActsAlignment::Alignment<fitter_t>::decompositionAnalysis(
0123     const AlignmentResult& res, std::ostream& out) {
0124   if (res.sumChi2SecondDerivative.cols() == 0) {
0125     ACTS_ERROR(
0126         "Please run Alignment::calculateAlignmentParameters before calling "
0127         "Alignment::decompositionAnalysis.");
0128     return -1;
0129   }
0130   Eigen::SelfAdjointEigenSolver<Acts::DynamicMatrix> eigenSolver(
0131       res.sumChi2SecondDerivative);
0132   if (eigenSolver.info() != Eigen::Success) {
0133     std::cout << " FAILED to find decompose correlation term" << std::endl;
0134     return -1;
0135   }
0136   const Acts::DynamicVector eigenVals = eigenSolver.eigenvalues();
0137   const Acts::DynamicMatrix eigenVecs = eigenSolver.eigenvectors();
0138
0139   std::map<double, int> sortedEV;
0140   for (int k = 0; k < eigenVals.size(); ++k) {
0141     sortedEV.emplace(eigenVals(k), k);
0142   }
0143   double firstEV = -1, lastEV = -1;
0144   for (auto& [EV, index] : sortedEV) {
0145     if (EV > 0 && firstEV < 0) {
0146       firstEV = EV;
0147     }
0148     lastEV = EV;
0149     out << " Eigenvector " << index << " has eigenvalue " << EV << std::endl;
0150     for (Eigen::Index row = 0; row < eigenVecs.rows(); ++row) {
0151       out << "        " << std::setw(12) << "  " << std::setw(3) << row + 1
0152           << "  " << std::setw(12) << eigenVecs(row, index) << std::endl;
0153     }
0154     out << std::endl;
0155   }
0156   return lastEV / firstEV;
0157 }
0158
0159 template <typename fitter_t>
0160 Acts::Result<void>
0161 ActsAlignment::Alignment<fitter_t>::updateAlignmentParameters(
0162     const Acts::GeometryContext& gctx,
0163     const std::vector<Acts::SurfacePlacementBase*>& alignedDetElements,
0164     const ActsAlignment::AlignedTransformUpdaterConcept auto&
0165         alignedTransformUpdater,
0166     ActsAlignment::AlignmentResult& alignResult) const {
0167   // Update the aligned transform
0168   Acts::AlignmentVector deltaAlignmentParam = Acts::AlignmentVector::Zero();
0169   for (const auto& [surface, index] : alignResult.idxedAlignSurfaces) {
0170     // 1. The original transform
0171     const Acts::Transform3& oldTransform =
0172         surface->localToGlobalTransform(gctx);
0173
0174     // 2. The delta transform
0175     deltaAlignmentParam = alignResult.deltaAlignmentParameters.segment(
0176         Acts::eAlignmentSize * index, Acts::eAlignmentSize);
0177     // Translation along the local x/y/z axes
0178     const Acts::Vector3 deltaCenterLocal =
0179         deltaAlignmentParam.segment<3>(Acts::eAlignmentCenter0);
0180     // The delta Euler angles
0181     Acts::Vector3 deltaEulerAngles =
0182         deltaAlignmentParam.segment<3>(Acts::eAlignmentRotation0);
0183
0184     // 3. The new transform. translate() shifts in the local frame.
0185     Acts::Transform3 newTransform = oldTransform;
0186     newTransform.translate(deltaCenterLocal);
0187     // Rotation first around fixed local x, then around fixed local y, and last
0188     // around fixed local z, this is the same as first around local z, then
0189     // around new local y, and last around new local x below
0190     newTransform *=
0191         Acts::AngleAxis3(deltaEulerAngles(2), Acts::Vector3::UnitZ());
0192     newTransform *=
0193         Acts::AngleAxis3(deltaEulerAngles(1), Acts::Vector3::UnitY());
0194     newTransform *=
0195         Acts::AngleAxis3(deltaEulerAngles(0), Acts::Vector3::UnitX());
0196
0197     // 4. Update the aligned transform
0198     //@Todo: use a better way to handle this (need dynamic cast to inherited
0199     // detector element type)
0200     ACTS_VERBOSE("Delta of alignment parameters at element "
0201                  << index << "= \n"
0202                  << deltaAlignmentParam);
0203     bool updated = alignedTransformUpdater(alignedDetElements.at(index), gctx,
0204                                            newTransform);
0205     if (!updated) {
0206       ACTS_ERROR("Update alignment parameters for detector element failed");
0207       return AlignmentError::AlignmentParametersUpdateFailure;
0208     }
0209   }
0210
0211   return Acts::Result<void>::success();
0212 }
0213
0214 template <typename fitter_t>
0215 template <typename trajectory_container_t,
0216           typename start_parameters_container_t, typename fit_options_t>
0217 Acts::Result<ActsAlignment::AlignmentResult>
0218 ActsAlignment::Alignment<fitter_t>::align(
0219     const trajectory_container_t& trajectoryCollection,
0220     const start_parameters_container_t& startParametersCollection,
0221     const ActsAlignment::AlignmentOptions<fit_options_t>& alignOptions) const {
0222   // Construct an AlignmentResult object
0223   AlignmentResult alignResult;
0224
0225   // Build the hierarchy registry and validate it. Mixed mode is supported:
0226   // a detector element may float as a standalone module (i.e. not appear in
0227   // any structure) and is then aligned at module level.
0228   AlignmentHierarchy hierarchy(alignOptions.alignedStructures);
0229   const auto validation = hierarchy.validate();
0230   if (!validation.ok()) {
0231     for (const auto* surface : validation.overlapping) {
0232       ACTS_ERROR("Surface " << surface->geometryId()
0233                             << " is assigned to multiple alignable structures");
0234     }
0235     return AlignmentError::HierarchyValidationFailure;
0236   }
0237
0238   // No-overlap-across-levels: for every iteration, a DoF that is floating at
0239   // structure level must not also float at module level on a child element.
0240   // Otherwise the Hessian becomes singular along the correlated direction.
0241   bool hadMaskConflicts = false;
0242   for (unsigned int iIter = 0; iIter < alignOptions.maxIterations; ++iIter) {
0243     AlignmentMask moduleMask = AlignmentMask::All;
0244     auto iter_it = alignOptions.iterationState.find(iIter);
0245     if (iter_it != alignOptions.iterationState.end()) {
0246       moduleMask = iter_it->second;
0247     }
0248     const auto conflicts = hierarchy.detectMaskConflicts(
0249         moduleMask, alignOptions.alignedDetElements);
0250     for (const auto& c : conflicts) {
0251       ACTS_ERROR("Detector element "
0252                  << c.detElement
0253                  << " (Surface ID: " << c.detElement->surface().geometryId()
0254                  << ") has DoFs floating at both structure level (ID "
0255                  << c.structure->geometryId()
0256                  << ") and module level in iteration " << iIter);
0257       hadMaskConflicts = true;
0258     }
0259   }
0260   if (hadMaskConflicts) {
0261     return AlignmentError::HierarchyValidationFailure;
0262   }
0263
0264   // Assign index to the alignable surface
0265   for (unsigned int iDetElement = 0;
0266        iDetElement < alignOptions.alignedDetElements.size(); iDetElement++) {
0267     alignResult.idxedAlignSurfaces.emplace(
0268         &alignOptions.alignedDetElements.at(iDetElement)->surface(),
0269         iDetElement);
0270   }
0271   ACTS_VERBOSE("There are " << alignResult.idxedAlignSurfaces.size()
0272                             << " detector elements to be aligned");
0273
0274   // Start the iteration to minimize the chi2
0275   bool converged = false;
0276   bool alignmentParametersUpdated = false;
0277   std::queue<double> recentChi2ONdf;
0278   ACTS_INFO("Max number of iterations: " << alignOptions.maxIterations);
0279   for (unsigned int iIter = 0; iIter < alignOptions.maxIterations; iIter++) {
0280     // Perform the fit to the trajectories and update alignment parameters
0281     // Initialize the alignment mask (all dof in default)
0282     AlignmentMask alignMask = AlignmentMask::All;
0283     // Set the alignment mask
0284     auto iter_it = alignOptions.iterationState.find(iIter);
0285     if (iter_it != alignOptions.iterationState.end()) {
0286       alignMask = iter_it->second;
0287     }
0288     // Calculate the alignment parameters delta etc.
0289     calculateAlignmentParameters(
0290         trajectoryCollection, startParametersCollection,
0291         alignOptions.fitOptions, alignResult, alignMask, &hierarchy);
0292     // Screen out the information
0293     ACTS_INFO("iIter = " << iIter << ", total chi2 = " << alignResult.chi2
0294                          << ", total measurementDim = "
0295                          << alignResult.measurementDim
0296                          << " and average chi2/ndf = "
0297                          << alignResult.averageChi2ONdf);
0298     // Check if it has converged against the provided precision
0299     // 1. either the delta average chi2/ndf in the last few
0300     // iterations is within tolerance
0301     if (recentChi2ONdf.size() >=
0302         alignOptions.deltaAverageChi2ONdfCutOff.first) {
0303       if (std::abs(recentChi2ONdf.front() - alignResult.averageChi2ONdf) <=
0304           alignOptions.deltaAverageChi2ONdfCutOff.second) {
0305         ACTS_INFO(
0306             "Alignment has converged with change of chi2/ndf < "
0307             << alignOptions.deltaAverageChi2ONdfCutOff.second << " in the last "
0308             << alignOptions.deltaAverageChi2ONdfCutOff.first << " iterations"
0309             << " after " << iIter << " iteration(s)");
0310         converged = true;
0311         break;
0312       }
0313       recentChi2ONdf.pop();
0314     }
0315     // 2. or the average chi2/ndf (is this correct?)
0316     if (alignResult.averageChi2ONdf <= alignOptions.averageChi2ONdfCutOff) {
0317       ACTS_INFO("Alignment has converged with average chi2/ndf < "
0318                 << alignOptions.averageChi2ONdfCutOff << " after " << iIter
0319                 << " iteration(s)");
0320       converged = true;
0321       break;
0322     }
0323     // Remove the first element
0324     // Store the result in the queue
0325     recentChi2ONdf.push(alignResult.averageChi2ONdf);
0326
0327     ACTS_INFO("The solved delta of alignmentParameters = \n "
0328               << alignResult.deltaAlignmentParameters);
0329     // Not coveraged yet, update the detector element alignment parameters
0330     auto updateRes = updateAlignmentParameters(
0331         alignOptions.fitOptions.geoContext, alignOptions.alignedDetElements,
0332         alignOptions.alignedTransformUpdater, alignResult);
0333     if (!updateRes.ok()) {
0334       ACTS_ERROR("Update alignment parameters failed: " << updateRes.error());
0335       return updateRes.error();
0336     }
0337     alignmentParametersUpdated = true;
0338   }  // end of all iterations
0339
0340   // Alignment failure if not converged
0341   if (!converged) {
0342     ACTS_ERROR("Alignment is not converged.");
0343     alignResult.result = AlignmentError::ConvergeFailure;
0344   }
0345
0346   // Screen out the final aligned parameters
0347   // @todo
0348   if (alignmentParametersUpdated) {
0349     for (const auto& det : alignOptions.alignedDetElements) {
0350       const auto& surface = &det->surface();
0351       const auto& transform =
0352           det->localToGlobalTransform(alignOptions.fitOptions.geoContext);
0353       // write it to the result
0354       alignResult.alignedParameters.emplace(det, transform);
0355       const auto& translation = transform.translation();
0356       const auto& rotation = transform.rotation();
0357       const Acts::Vector3 rotAngles =
0358           Acts::detail::EigenCompat::canonicalEulerAngles(rotation, 2, 1, 0);
0359       ACTS_VERBOSE("Detector element with surface "
0360                    << surface->geometryId()
0361                    << " has aligned geometry position as below:");
0362       ACTS_VERBOSE("Center (cenX, cenY, cenZ) = " << translation.transpose());
0363       ACTS_VERBOSE(
0364           "Euler angles (rotZ, rotY, rotX) = " << rotAngles.transpose());
0365       ACTS_VERBOSE("Rotation matrix = \n" << rotation);
0366     }
0367   } else {
0368     ACTS_DEBUG("Alignment parameters is not updated.");
0369   }
0370
0371   return alignResult;
0372 }