Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-14 09:39:38 

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/Seeding/CompositeSpacePointLineFitter.hpp"
0012 
0013 #include "Acts/Utilities/AlgebraHelpers.hpp"
0014 #include "Acts/Utilities/StringHelpers.hpp"
0015 
0016 #include <format>
0017 
0018 namespace Acts::Experimental {
0019 
0020 template <CompositeSpacePointContainer Cont_t>
0021 CompositeSpacePointLineFitter::DoFcounts
0022 CompositeSpacePointLineFitter::countDoF(const Cont_t& measurements) const {
0023   return countDoF(measurements, Selector_t<SpacePoint_t<Cont_t>>{});
0024 }
0025 
0026 template <CompositeSpacePointContainer Cont_t>
0027 CompositeSpacePointLineFitter::DoFcounts
0028 CompositeSpacePointLineFitter::countDoF(
0029     const Cont_t& measurements,
0030     const Selector_t<SpacePoint_t<Cont_t>>& selector) const {
0031   DoFcounts counts{};
0032   std::size_t nValid{0};
0033   for (const auto& sp : measurements) {
0034     if (selector.connected() && !selector(*sp)) {
0035       ACTS_VERBOSE(__func__ << "() " << __LINE__
0036                             << " -  Skip invalid measurement @"
0037                             << toString(sp->localPosition()));
0038       continue;
0039     }
0040     ++nValid;
0041     if (sp->measuresLoc0()) {
0042       ++counts.nonBending;
0043     }
0044     if (sp->measuresLoc1()) {
0045       ++counts.bending;
0046     }
0047     if (sp->isStraw()) {
0048       ++counts.straw;
0049     } else if (sp->hasTime()) {
0050       ++counts.time;
0051     }
0052   }
0053   ACTS_DEBUG(__func__ << "() " << __LINE__ << " - " << nValid << "/"
0054                       << measurements.size()
0055                       << " valid measurements passed. Found "
0056                       << counts.nonBending << "/" << counts.bending << "/"
0057                       << counts.time << " measuring loc0/loc1/time with "
0058                       << counts.straw << " straw measurements.");
0059   return counts;
0060 }
0061 
0062 template <CompositeSpacePointContainer Cont_t>
0063 CompositeSpacePointLineFitter::FastFitResult
0064 CompositeSpacePointLineFitter::fastPrecFit(
0065     const Cont_t& measurements, const Line_t& initialGuess,
0066     const std::size_t nStraws,
0067     const std::vector<FitParIndex>& parsToUse) const {
0068   if (std::ranges::none_of(parsToUse, [](const FitParIndex idx) {
0069         using enum FitParIndex;
0070         return idx == theta || idx == y0;
0071       })) {
0072     return std::nullopt;
0073   }
0074   ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fetched " << nStraws
0075                       << " measurements for fast fit without time.");
0076 
0077   if (nStraws > 2) {
0078     std::vector<int> signs = detail::CompSpacePointAuxiliaries::strawSigns(
0079         initialGuess, measurements);
0080     FastFitResult precResult{m_fastFitter.fit(measurements, signs)};
0081 
0082     // Fast fit gave a bad chi2 -> Maybe L/R solution is swapped?
0083     if (precResult && precResult->chi2 / static_cast<double>(precResult->nDoF) >
0084                           m_cfg.badFastChi2SignSwap) {
0085       // Swap signs
0086       ACTS_DEBUG(__func__ << "() " << __LINE__
0087                           << " - The fit result is of poor quality "
0088                           << (*precResult) << " attempt with L<->R swapping");
0089       for (auto& s : signs) {
0090         s = -s;
0091       }
0092       // Retry & check whether the chi2 is better
0093       FastFitResult swappedPrecResult = {m_fastFitter.fit(measurements, signs)};
0094       if (swappedPrecResult && swappedPrecResult->chi2 < precResult->chi2) {
0095         ACTS_DEBUG(__func__ << "() " << __LINE__
0096                             << " - Swapped fit is of better quality "
0097                             << (*swappedPrecResult));
0098         swappedPrecResult->nIter += precResult->nIter;
0099         return swappedPrecResult;
0100       } else {
0101         precResult->nIter += swappedPrecResult->nIter;
0102         ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fit did not improve "
0103                             << (*swappedPrecResult));
0104       }
0105     }
0106     return precResult;
0107   }
0108 
0109   // Not enough straw measurements, proceeding with strip fast fit
0110   using ResidualIdx = detail::CompSpacePointAuxiliaries::ResidualIdx;
0111   return m_fastFitter.fit(measurements, ResidualIdx::bending);
0112 }
0113 
0114 template <CompositeSpacePointContainer Cont_t,
0115           CompositeSpacePointFastCalibrator<
0116               CompositeSpacePointLineFitter::SpacePoint_t<Cont_t>>
0117               Calibrator_t>
0118 CompositeSpacePointLineFitter::FastFitResultT0
0119 CompositeSpacePointLineFitter::fastPrecFit(
0120     const Acts::CalibrationContext& ctx, const Calibrator_t& calibrator,
0121     const Cont_t& measurements, const Line_t& initialGuess,
0122     const double initialT0, const std::vector<FitParIndex>& parsToUse) const {
0123   if (std::ranges::none_of(parsToUse, [](const FitParIndex idx) {
0124         using enum FitParIndex;
0125         return idx == theta || idx == y0;
0126       })) {
0127     return std::nullopt;
0128   }
0129   ACTS_DEBUG(__func__ << "() " << __LINE__
0130                       << " Starting precision fast fit with time.");
0131 
0132   std::vector<int> signs =
0133       detail::CompSpacePointAuxiliaries::strawSigns(initialGuess, measurements);
0134   FastFitResultT0 precResult{
0135       m_fastFitter.fit(ctx, calibrator, measurements, signs, initialT0)};
0136 
0137   // Fast fit gave a bad chi2 -> Maybe L/R solution is swapped?
0138   if (precResult && precResult->chi2 / static_cast<double>(precResult->nDoF) >
0139                         m_cfg.badFastChi2SignSwap) {
0140     // Swap signs
0141     ACTS_DEBUG(__func__ << "() " << __LINE__
0142                         << " - The fit result is of poor quality "
0143                         << (*precResult) << " attempt with L<->R swapping");
0144     for (auto& s : signs) {
0145       s = -s;
0146     }
0147     // Retry & check whether the chi2 is better
0148     FastFitResultT0 swappedPrecResult{
0149         m_fastFitter.fit(ctx, calibrator, measurements, signs, initialT0)};
0150     if (swappedPrecResult && swappedPrecResult->chi2 < precResult->chi2) {
0151       ACTS_DEBUG(__func__ << "() " << __LINE__
0152                           << " - Swapped fit is of better quality "
0153                           << (*swappedPrecResult));
0154       swappedPrecResult->nIter += precResult->nIter;
0155       return swappedPrecResult;
0156     } else {
0157       precResult->nIter += swappedPrecResult->nIter;
0158       ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fit did not improve "
0159                           << (*swappedPrecResult));
0160     }
0161   }
0162   return precResult;
0163 }
0164 
0165 template <CompositeSpacePointContainer Cont_t>
0166 CompositeSpacePointLineFitter::FitParameters
0167 CompositeSpacePointLineFitter::fastFit(
0168     const Cont_t& measurements, const Line_t& initialGuess,
0169     const std::size_t nStraws,
0170     const std::vector<FitParIndex>& parsToUse) const {
0171   using enum FitParIndex;
0172 
0173   FitParameters result{};
0174   const FastFitResult precResult{
0175       fastPrecFit(measurements, initialGuess, nStraws, parsToUse)};
0176   if (!precResult) {
0177     ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fast fit failed.");
0178     return result;
0179   }
0180 
0181   // Copy spatial parameters & covariance from precision fit and
0182   // perform non-bending fit (if required)
0183   mergePrecAndNonPrec(result, precResult, measurements, parsToUse);
0184 
0185   return result;
0186 }
0187 
0188 template <CompositeSpacePointContainer Cont_t,
0189           CompositeSpacePointFastCalibrator<
0190               CompositeSpacePointLineFitter::SpacePoint_t<Cont_t>>
0191               Calibrator_t>
0192 CompositeSpacePointLineFitter::FitParameters
0193 CompositeSpacePointLineFitter::fastFit(
0194     const Acts::CalibrationContext& ctx, const Calibrator_t& calibrator,
0195     const Cont_t& measurements, const Line_t& initialGuess,
0196     const double startT0, const std::vector<FitParIndex>& parsToUse) const {
0197   using enum FitParIndex;
0198 
0199   FitParameters result{};
0200   const FastFitResultT0 precResult{fastPrecFit(
0201       ctx, calibrator, measurements, initialGuess, startT0, parsToUse)};
0202   if (!precResult) {
0203     ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fast fit failed.");
0204     return result;
0205   }
0206   // Copy t0 parameter & covariance from precision fit
0207   result.parameters[toUnderlying(t0)] = precResult->t0;
0208   result.covariance(toUnderlying(t0), toUnderlying(t0)) =
0209       Acts::square(precResult->dT0);
0210 
0211   // Copy spatial parameters & covariance from precision fit and
0212   // perform non-bending fit (if required)
0213   const bool fitT0{true};
0214   mergePrecAndNonPrec(result, precResult, measurements, parsToUse, fitT0);
0215 
0216   return result;
0217 }
0218 
0219 template <CompositeSpacePointContainer Cont_t>
0220 void CompositeSpacePointLineFitter::mergePrecAndNonPrec(
0221     FitParameters& result, const FastFitResult& precResult,
0222     const Cont_t& measurements, const std::vector<FitParIndex>& parsToUse,
0223     const bool fitT0) const {
0224   using enum FitParIndex;
0225   using namespace Acts::UnitLiterals;
0226 
0227   // Copy the parameters & covariance from precision fit result
0228   result.parameters[toUnderlying(y0)] = precResult->y0;
0229   result.parameters[toUnderlying(theta)] = precResult->theta;
0230   result.covariance(toUnderlying(y0), toUnderlying(y0)) =
0231       Acts::square(precResult->dY0);
0232   result.covariance(toUnderlying(theta), toUnderlying(theta)) =
0233       Acts::square(precResult->dTheta);
0234   result.nDoF = static_cast<unsigned>(precResult->nDoF);
0235   result.nIter = static_cast<unsigned>(precResult->nIter);
0236   result.chi2 = precResult->chi2;
0237   result.converged = true;
0238 
0239   // Check whether a non-bending fit is required
0240   if (std::ranges::none_of(parsToUse, [](const FitParIndex idx) {
0241         return idx == phi || idx == x0;
0242       })) {
0243     result.parameters[toUnderlying(phi)] = 90._degree;
0244     ACTS_DEBUG(
0245         __func__ << "() " << __LINE__
0246                  << " - No measurements in non precision direction parsed.");
0247     return;
0248   }
0249 
0250   ACTS_DEBUG(__func__ << "() " << __LINE__
0251                       << " - Start fast non-precision fit.");
0252   using ResidualIdx = detail::CompSpacePointAuxiliaries::ResidualIdx;
0253   const FastFitResult nonPrecResult =
0254       m_fastFitter.fit(measurements, ResidualIdx::nonBending);
0255   if (!nonPrecResult) {
0256     result.parameters[toUnderlying(phi)] = 90._degree;
0257     ACTS_DEBUG(__func__ << "() " << __LINE__
0258                         << " - Fast non-precision fit failed.");
0259     return;
0260   }
0261   result.parameters[toUnderlying(x0)] = nonPrecResult->y0;
0262   result.covariance(toUnderlying(x0), toUnderlying(x0)) =
0263       Acts::square(nonPrecResult->dY0);
0264   result.nDoF += nonPrecResult->nDoF;
0265   result.nIter += nonPrecResult->nIter;
0266 
0267   // Combine the two results into a single direction
0268   double tanTheta = std::tan(result.parameters[toUnderlying(theta)]);
0269   double tanPhi = std::tan(nonPrecResult->theta);
0270   auto dir = makeDirectionFromAxisTangents(tanPhi, tanTheta);
0271   result.parameters[toUnderlying(phi)] = VectorHelpers::phi(dir);
0272   result.parameters[toUnderlying(theta)] = VectorHelpers::theta(dir);
0273 
0274   // Recompute the chi2
0275   Line_t recoLine{result.parameters};
0276   result.chi2 = 0.;
0277   std::ranges::for_each(
0278       measurements, [&recoLine, &result, &fitT0](const auto& m) {
0279         result.chi2 +=
0280             fitT0 ? detail::CompSpacePointAuxiliaries::chi2Term(
0281                         recoLine, result.parameters[toUnderlying(t0)], *m)
0282                   : detail::CompSpacePointAuxiliaries::chi2Term(recoLine, *m);
0283       });
0284   ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Fast fit done. Obtained result"
0285                       << result);
0286 }
0287 
0288 template <CompositeSpacePointContainer Cont_t,
0289           CompositeSpacePointCalibrator<Cont_t, Cont_t> Calibrator_t>
0290 CompositeSpacePointLineFitter::FitResult<Cont_t>
0291 CompositeSpacePointLineFitter::fit(
0292     FitOptions<Cont_t, Calibrator_t>&& fitOpts) const {
0293   using namespace Acts::UnitLiterals;
0294 
0295   if (!fitOpts.calibrator) {
0296     throw std::invalid_argument(
0297         "CompositeSpacePointLineFitter::fit() - Please provide a valid pointer "
0298         "to a calibrator.");
0299   }
0300 
0301   FitResult<Cont_t> result{};
0302   result.measurements = std::move(fitOpts.measurements);
0303   result.parameters = std::move(fitOpts.startParameters);
0304 
0305   // Declare the auxiliaries object to calculate the residuals
0306   detail::CompSpacePointAuxiliaries::Config resCfg{};
0307   resCfg.localToGlobal = std::move(fitOpts.localToGlobal);
0308   resCfg.useHessian = m_cfg.useHessian;
0309   resCfg.calcAlongStraw = m_cfg.calcAlongStraw;
0310   resCfg.calcAlongStrip = m_cfg.calcAlongStrip;
0311   resCfg.includeToF = m_cfg.includeToF;
0312 
0313   DoFcounts hitCounts{countDoF(result.measurements, fitOpts.selector)};
0314   const std::size_t nStraws{hitCounts.straw};
0315   resCfg.parsToUse = extractFitablePars(hitCounts);
0316   if (resCfg.parsToUse.empty()) {
0317     ACTS_WARNING(__func__ << "() " << __LINE__
0318                           << ": No valid degrees of freedom parsed. Please "
0319                           << "check your measurements");
0320     return result;
0321   }
0322 
0323   Line_t line{};
0324 
0325   // Fast fitter shall be used
0326   if (m_cfg.useFastFitter) {
0327     line.updateParameters(result.parameters);
0328     ACTS_DEBUG(__func__ << "() " << __LINE__
0329                         << " - Attempt a fast fit, first.");
0330     const bool fitTime{resCfg.parsToUse.back() == FitParIndex::t0};
0331     FitParameters fastResult{};
0332     if (fitTime && nStraws > 2) {
0333       if constexpr (CompositeSpacePointFastCalibrator<Calibrator_t,
0334                                                       SpacePoint_t<Cont_t>>) {
0335         fastResult = fastFit(fitOpts.calibContext, *fitOpts.calibrator,
0336                              result.measurements, line,
0337                              result.parameters[toUnderlying(FitParIndex::t0)],
0338                              resCfg.parsToUse);
0339         if (fastResult.converged) {
0340           fastResult.parameters[toUnderlying(FitParIndex::t0)] -=
0341               (fitOpts.localToGlobal * line.position()).norm() /
0342               PhysicalConstants::c;
0343         }
0344       } else {
0345         ACTS_WARNING(__func__
0346                      << "() " << __LINE__
0347                      << ": Skipping fast fit - calibrator does not satisfy "
0348                      << "CompositeSpacePointFastCalibrator concept");
0349       }
0350     } else {
0351       fastResult =
0352           fastFit(result.measurements, line, nStraws, resCfg.parsToUse);
0353     }
0354     if (fastResult.converged) {
0355       static_cast<FitParameters&>(result) = std::move(fastResult);
0356       ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fit converged.");
0357       // Use the result from the fast fitter as final answer
0358       if (!m_cfg.fastPreFitter) {
0359         line.updateParameters(result.parameters);
0360         // Recalibrate the measurements before returning
0361         result.measurements = fitOpts.calibrator->calibrate(
0362             fitOpts.calibContext, line.position(), line.direction(),
0363             fitTime ? result.parameters[toUnderlying(FitParIndex::t0)] : 0.,
0364             result.measurements);
0365         return result;
0366       }
0367       // Set convergence flag to false because the full fit comes later.
0368       result.converged = false;
0369     } else {
0370       ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Fit failed.");
0371       return result;
0372     }
0373   }
0374   // Update the drift signs if no re-calibration per iteration is scheduled
0375   if (!m_cfg.recalibrate) {
0376     line.updateParameters(result.parameters);
0377     fitOpts.calibrator->updateSigns(line.position(), line.direction(),
0378                                     result.measurements);
0379   }
0380   /// Proceed with the usual fit
0381   ChiSqCache cache{};
0382   detail::CompSpacePointAuxiliaries pullCalculator{resCfg, logger().clone()};
0383   for (; !result.converged && result.nIter < m_cfg.maxIter; ++result.nIter) {
0384     cache.reset();
0385     // Update the parameters from the last iteration
0386     line.updateParameters(result.parameters);
0387     const double t0 = result.parameters[toUnderlying(FitParIndex::t0)];
0388 
0389     ACTS_VERBOSE(__func__ << "() " << __LINE__ << ": Start next iteration "
0390                           << result << " --> " << toString(line.position())
0391                           << " + " << toString(line.direction()) << ".");
0392     // Update the measurements if calibration loop is switched on
0393     if (m_cfg.recalibrate) {
0394       result.measurements = fitOpts.calibrator->calibrate(
0395           fitOpts.calibContext, line.position(), line.direction(), t0,
0396           result.measurements);
0397       // Check whether the measurements are still valid
0398       auto parsBkp = std::move(resCfg.parsToUse);
0399       resCfg.parsToUse =
0400           extractFitablePars(countDoF(result.measurements, fitOpts.selector));
0401 
0402       // No valid measurement is left
0403       if (resCfg.parsToUse.empty()) {
0404         ACTS_WARNING(__func__ << "() " << __LINE__ << ":  Line parameters "
0405                               << toString(line.position()) << " + "
0406                               << toString(line.direction()) << ", t0: "
0407                               << t0 / 1._ns << " invalidated all measurements");
0408         return result;
0409       }
0410 
0411       if (parsBkp != resCfg.parsToUse) {
0412         // Instantiate an updated pull calculator
0413         pullCalculator =
0414             detail::CompSpacePointAuxiliaries{resCfg, logger().clone()};
0415       }
0416       // update the drift signs
0417       fitOpts.calibrator->updateSigns(line.position(), line.direction(),
0418                                       result.measurements);
0419     }
0420     // Calculate the new chi2
0421     for (const auto& spacePoint : result.measurements) {
0422       // Skip bad measurements
0423       if (fitOpts.selector.connected() && !fitOpts.selector(*spacePoint)) {
0424         continue;
0425       }
0426       // Calculate the residual & derivatives
0427       if (pullCalculator.config().parsToUse.back() == FitParIndex::t0) {
0428         double driftV{fitOpts.calibrator->driftVelocity(fitOpts.calibContext,
0429                                                         *spacePoint)};
0430         double driftA{fitOpts.calibrator->driftAcceleration(
0431             fitOpts.calibContext, *spacePoint)};
0432         pullCalculator.updateFullResidual(line, t0, *spacePoint, driftV,
0433                                           driftA);
0434       } else {
0435         pullCalculator.updateSpatialResidual(line, *spacePoint);
0436       }
0437       // Propagte to the chi2
0438       pullCalculator.updateChiSq(cache, spacePoint->covariance());
0439     }
0440     pullCalculator.symmetrizeHessian(cache);
0441     result.chi2 = cache.chi2;
0442     // Now update the parameters
0443     UpdateStep update{UpdateStep::goodStep};
0444     switch (pullCalculator.config().parsToUse.size()) {
0445       // 2D fit (intercept + inclination angle)
0446       case 2: {
0447         update = updateParameters<2>(pullCalculator.config().parsToUse.front(),
0448                                      cache, result.parameters);
0449         break;
0450       }
0451       // 2D fit + time
0452       case 3: {
0453         update = updateParameters<3>(pullCalculator.config().parsToUse.front(),
0454                                      cache, result.parameters);
0455         break;
0456       }
0457       // 3D spatial fit (x0, y0, theta, phi)
0458       case 4: {
0459         update = updateParameters<4>(pullCalculator.config().parsToUse.front(),
0460                                      cache, result.parameters);
0461         break;
0462       }
0463       // full fit
0464       case 5: {
0465         update = updateParameters<5>(pullCalculator.config().parsToUse.front(),
0466                                      cache, result.parameters);
0467         break;
0468       }
0469       default:
0470         ACTS_WARNING(__func__ << "() " << __LINE__
0471                               << ": Invalid parameter size "
0472                               << pullCalculator.config().parsToUse.size());
0473         return result;
0474     }
0475     /// Check whether the fit parameters are within range
0476     for (const FitParIndex par : pullCalculator.config().parsToUse) {
0477       const auto p = toUnderlying(par);
0478       if (m_cfg.ranges[p][0] < m_cfg.ranges[p][1] &&
0479           (result.parameters[p] < m_cfg.ranges[p][0] ||
0480            result.parameters[p] > m_cfg.ranges[p][1])) {
0481         ACTS_VERBOSE(__func__ << "() " << __LINE__ << ": The parameter "
0482                               << pullCalculator.parName(par) << " "
0483                               << result.parameters[p] << " is out range ["
0484                               << m_cfg.ranges[p][0] << ";" << m_cfg.ranges[p][1]
0485                               << "]");
0486         update = UpdateStep::outOfBounds;
0487       }
0488     }
0489     switch (update) {
0490       using enum UpdateStep;
0491       case converged: {
0492         result.converged = true;
0493         break;
0494       }
0495       case goodStep: {
0496         break;
0497       }
0498       case outOfBounds: {
0499         return result;
0500       }
0501     }
0502   }
0503   // Parameters converged
0504   if (result.converged) {
0505     const DoFcounts doF = countDoF(result.measurements, fitOpts.selector);
0506     result.nDoF = doF.nonBending + doF.bending + doF.time -
0507                   pullCalculator.config().parsToUse.size();
0508     line.updateParameters(result.parameters);
0509     const double t0 = result.parameters[toUnderlying(FitParIndex::t0)];
0510     // Recalibrate the measurements before returning
0511     result.measurements = fitOpts.calibrator->calibrate(
0512         fitOpts.calibContext, line.position(), line.direction(), t0,
0513         result.measurements);
0514 
0515     // Assign the Hessian
0516     switch (pullCalculator.config().parsToUse.size()) {
0517       // 2D fit (intercept + inclination angle)
0518       case 2: {
0519         fillCovariance<2>(pullCalculator.config().parsToUse.front(),
0520                           cache.hessian, result.covariance);
0521         break;
0522       }
0523       // 2D fit + time
0524       case 3: {
0525         fillCovariance<3>(pullCalculator.config().parsToUse.front(),
0526                           cache.hessian, result.covariance);
0527         break;
0528       }
0529       // 3D spatial fit (x0, y0, theta, phi)
0530       case 4: {
0531         fillCovariance<4>(pullCalculator.config().parsToUse.front(),
0532                           cache.hessian, result.covariance);
0533         break;
0534       }
0535       // full fit
0536       case 5: {
0537         fillCovariance<5>(pullCalculator.config().parsToUse.front(),
0538                           cache.hessian, result.covariance);
0539         break;
0540       }
0541       // No need to warn here -> captured by the fit iterations
0542       default:
0543         break;
0544     }
0545   }
0546   return result;
0547 }
0548 
0549 template <unsigned N>
0550 CompositeSpacePointLineFitter::UpdateStep
0551 CompositeSpacePointLineFitter::updateParameters(const FitParIndex firstPar,
0552                                                 ChiSqCache& cache,
0553                                                 ParamVec_t& currentPars) const
0554   requires(N >= 2 && N <= s_nPars)
0555 {
0556   auto firstIdx = toUnderlying(firstPar);
0557   assert(firstIdx + N <= s_nPars);
0558   // Current parameters mapped to an Eigen interface
0559   if constexpr (N == 3) {
0560     auto t0Idx = toUnderlying(FitParIndex::t0);
0561     assert(firstIdx + N - 1 == 2);
0562     std::swap(currentPars[2], currentPars[t0Idx]);
0563     std::swap(cache.gradient[2], cache.gradient[t0Idx]);
0564     cache.hessian.row(2).swap(cache.hessian.row(t0Idx));
0565     cache.hessian.col(2).swap(cache.hessian.col(t0Idx));
0566   }
0567   Eigen::Map<ActsVector<N>> miniPars{currentPars.data() + firstIdx};
0568   ACTS_VERBOSE(__func__ << "<" << N << ">() - " << __LINE__
0569                         << ": Current parameters " << toString(miniPars)
0570                         << " with chi2: " << cache.chi2 << ",  gradient: "
0571                         << toString(cache.gradient) << ", hessian: \n"
0572                         << cache.hessian);
0573 
0574   // Take out the filled block from the gradient
0575   Eigen::Map<const ActsVector<N>> miniGradient{cache.gradient.data() +
0576                                                firstIdx};
0577   // The gradient is already small enough
0578   if (miniGradient.norm() < m_cfg.precCutOff) {
0579     ACTS_DEBUG(__func__ << "<" << N << ">() - " << __LINE__
0580                         << ": Gradient is small enough");
0581     return UpdateStep::converged;
0582   }
0583   // Take out the filled block from the hessian
0584   Acts::ActsSquareMatrix<N> miniHessian{
0585       cache.hessian.block<N, N>(firstIdx, firstIdx)};
0586   ACTS_VERBOSE(__func__ << "<" << N << ">() - " << __LINE__
0587                         << ": Projected parameters: " << toString(miniPars)
0588                         << " gradient: " << toString(miniGradient)
0589                         << ", hessian: \n"
0590                         << miniHessian
0591                         << ", determinant: " << miniHessian.determinant());
0592 
0593   auto inverseH = safeInverse(miniHessian);
0594   // The Hessian can safely be inverted
0595   if (inverseH) {
0596     const ActsVector<N> update{(*inverseH) * miniGradient};
0597 
0598     if (update.norm() < m_cfg.precCutOff) {
0599       ACTS_DEBUG(__func__ << "<" << N << ">() - " << __LINE__ << ": Update "
0600                           << toString(update) << " is negligible small.");
0601       return UpdateStep::converged;
0602     }
0603     ACTS_VERBOSE(__func__ << "<" << N << ">() - " << __LINE__
0604                           << ": Inverted Hessian \n"
0605                           << (*inverseH) << "\n-> Update parameters by "
0606                           << toString(update));
0607     miniPars -= update;
0608 
0609   } else {
0610     // Fall back to gradient decent with a fixed damping factor
0611     const ActsVector<N> update{
0612         std::min(m_cfg.gradientStep, miniGradient.norm()) *
0613         miniGradient.normalized()};
0614 
0615     ACTS_VERBOSE(__func__ << "<" << N << ">() - " << __LINE__
0616                           << ": Update parameters by " << toString(update));
0617     miniPars -= update;
0618   }
0619   // swap back t0 component if needed
0620   if constexpr (N == 3) {
0621     std::swap(currentPars[2], currentPars[toUnderlying(FitParIndex::t0)]);
0622   }
0623   // Check parameter ranges
0624 
0625   return UpdateStep::goodStep;
0626 }
0627 
0628 template <unsigned N>
0629 void CompositeSpacePointLineFitter::fillCovariance(const FitParIndex firstPar,
0630                                                    const CovMat_t& hessian,
0631                                                    CovMat_t& covariance) const
0632   requires(N >= 2 && N <= s_nPars)
0633 {
0634   auto firstIdx = toUnderlying(firstPar);
0635   assert(firstIdx + N <= s_nPars);
0636 
0637   Acts::ActsSquareMatrix<N> miniHessian =
0638       hessian.block<N, N>(firstIdx, firstIdx).eval();
0639 
0640   auto inverseH = safeInverse(miniHessian);
0641   // The Hessian can safely be inverted
0642   if (inverseH) {
0643     auto covBlock = covariance.block<N, N>(firstIdx, firstIdx);
0644     covBlock = *inverseH;
0645 
0646     // swap back t0 component if needed
0647     if constexpr (N == 3) {
0648       auto t0Idx = toUnderlying(FitParIndex::t0);
0649       covariance(t0Idx, t0Idx) = 1.;
0650       covariance.row(2).swap(covariance.row(t0Idx));
0651       covariance.col(2).swap(covariance.col(t0Idx));
0652     }
0653   }
0654 
0655   ACTS_DEBUG(__func__ << "<" << N << ">() - " << __LINE__
0656                       << ": Evaluated covariance: \n"
0657                       << covariance);
0658 }
0659 
0660 }  // namespace Acts::Experimental
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!