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