Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:23:17 

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 #include "Acts/Seeding/detail/FastStrawLineFitter.hpp"
0010 
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Utilities/StringHelpers.hpp"
0013 #include "Acts/Utilities/detail/periodic.hpp"
0014 
0015 #include <format>
0016 
0017 using namespace Acts::UnitLiterals;
0018 
0019 namespace {
0020 /// @brief Expresses an angle in degree
0021 /// @param x: Angle in radians
0022 constexpr double inDeg(const double x) {
0023   return x / 1._degree;
0024 }
0025 /// @brief Express a time in units of nanoseconds
0026 constexpr double inNanoS(const double x) {
0027   return x / 1._ns;
0028 }
0029 std::string printThetaStep(const double theta, const double thetaPrime,
0030                            const double thetaTwoPrime) {
0031   return std::format(
0032       "-- theta: {:.3f}, thetaPrime: {:.3f}, thetaTwoPrime: {:.3f}  --> "
0033       "step-size: {:.8f}",
0034       inDeg(theta), inDeg(thetaPrime), inDeg(thetaTwoPrime),
0035       inDeg(thetaPrime / thetaTwoPrime));
0036 }
0037 }  // namespace
0038 namespace Acts::Experimental::detail {
0039 
0040 void FastStrawLineFitter::FitAuxiliaries::print(std::ostream& ostr) const {
0041   ostr << std::format("T_zzyy: {:.3f}, ", T_zzyy)
0042        << std::format("T_yz: {:.3f}, ", T_yz)
0043        << std::format("T_rz: {:.3f}, ", T_rz)
0044        << std::format("T_ry: {:.3f}, ", T_ry)
0045        << std::format("centre ( {:.3f}, {:3f}), ", centerY, centerZ)
0046        << std::format("y0: {:.3f}", fitY0)
0047        << std::format(", norm: {:.3f}", covNorm) << ", nDoF: " << nDoF;
0048 }
0049 void FastStrawLineFitter::FitAuxiliariesWithT0 ::print(
0050     std::ostream& ostr) const {
0051   ostr << std::format("R_vr: {:.3f}, ", R_vr);
0052   ostr << std::format("R_ar: {:.3f}, ", R_ar);
0053   ostr << std::format("R_vv: {:.3f}, ", R_vv);
0054   ostr << std::format("R_a: {:.3f}, ", R_a);
0055   ostr << std::format("R_v: {:.3f}, ", R_v);
0056   ostr << std::format("T_vz: {:.3f}, ", T_vz);
0057   ostr << std::format("T_vy: {:.3f}, ", T_vy);
0058   ostr << std::format("T_az: {:.3f}, ", T_az);
0059   ostr << std::format("T_ay: {:.3f}, ", T_ay);
0060   FitAuxiliaries::print(ostr);
0061 }
0062 
0063 void FastStrawLineFitter::FitResult::print(std::ostream& ostr) const {
0064   ostr << "# iteration: " << nIter << ", nDoF: " << nDoF << ", chi2: " << chi2
0065        << ", chi2 / nDoF: " << (chi2 / static_cast<double>(nDoF)) << ",\n";
0066   ostr << std::format("theta: {:.3f} pm {:.3f}, ", inDeg(theta), inDeg(dTheta))
0067        << std::format("y0: {:.3f}  pm {:.3f}", y0, dY0);
0068 }
0069 void FastStrawLineFitter::FitResultT0::print(std::ostream& ostr) const {
0070   FitResult::print(ostr);
0071   ostr << std::format(", t0: {:.3f} pm {:.3f}", inNanoS(t0), inNanoS(dT0));
0072 }
0073 
0074 using Vector = FastStrawLineFitter::Vector;
0075 const Acts::Logger& FastStrawLineFitter::logger() const {
0076   assert(m_logger != nullptr);
0077   return *m_logger;
0078 }
0079 FastStrawLineFitter::FastStrawLineFitter(const Config& cfg,
0080                                          std::unique_ptr<const Logger> logger)
0081     : m_cfg{cfg}, m_logger{std::move(logger)} {}
0082 
0083 inline void FastStrawLineFitter::calcAngularDerivatives(
0084     const TrigonomHelper& angles, const FitAuxiliaries& fitPars,
0085     double& thetaPrime, double& thetaTwoPrime) const {
0086   thetaPrime = 0.5 * fitPars.T_zzyy * angles.sinTwoTheta -
0087                fitPars.T_yz * angles.cosTwoTheta -
0088                fitPars.T_rz * angles.cosTheta - fitPars.T_ry * angles.sinTheta;
0089   thetaTwoPrime = fitPars.T_zzyy * angles.cosTwoTheta +
0090                   2. * fitPars.T_yz * angles.sinTwoTheta +
0091                   fitPars.T_rz * angles.sinTheta -
0092                   fitPars.T_ry * angles.cosTheta;
0093 }
0094 double FastStrawLineFitter::startTheta(const FitAuxiliaries& fitPars) {
0095   double thetaGuess =
0096       std::atan2(2. * (fitPars.T_yz - fitPars.T_ry), fitPars.T_zzyy) / 2.;
0097   return Acts::detail::wrap_periodic(thetaGuess, 0., std::numbers::pi);
0098 }
0099 double FastStrawLineFitter::calcTimeGrad(const TrigonomHelper& angles,
0100                                          const FitAuxiliariesWithT0& pars) {
0101   return pars.fitY0 * pars.R_v - pars.R_vr + pars.T_vz * angles.sinTheta -
0102          pars.T_vy * angles.cosTheta;
0103 }
0104 
0105 void FastStrawLineFitter::completeResult(const FitAuxiliaries& fitPars,
0106                                          const double thetaTwoPrime,
0107                                          FitResult& result) const {
0108   result.dTheta = std::sqrt(1. / Acts::abs(thetaTwoPrime));
0109   const double tanTheta = std::tan(result.theta);
0110   const double secTheta = 1. / std::cos(result.theta);
0111   result.y0 =
0112       fitPars.centerY - fitPars.centerZ * tanTheta + fitPars.fitY0 * secTheta;
0113   result.dY0 = Acts::fastHypot(-fitPars.centerZ * Acts::square(secTheta) +
0114                                    fitPars.fitY0 * secTheta * tanTheta,
0115                                secTheta);
0116 
0117   result.dY0 *= result.dTheta;
0118   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Fit succeeded " << result);
0119 }
0120 
0121 std::optional<FastStrawLineFitter::FitResult> FastStrawLineFitter::fit(
0122     const FitAuxiliaries& fitPars) const {
0123   /// No degrees of freedom -> no valid parameters
0124   if (fitPars.nDoF == 0u) {
0125     return std::nullopt;
0126   }
0127   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Input fit parameters "
0128                       << fitPars);
0129 
0130   const double thetaGuess = startTheta(fitPars);
0131   ACTS_DEBUG(__func__ << "() - " << __LINE__
0132                       << ": Start fast fit seed  theta: " << inDeg(thetaGuess));
0133   ////
0134   FitResult result{};
0135   result.theta = thetaGuess;
0136   result.nDoF = fitPars.nDoF;
0137 
0138   double thetaPrime{0.};
0139   double thetaTwoPrime{0.};
0140   while (result.nIter <= m_cfg.maxIter) {
0141     ++result.nIter;
0142     const TrigonomHelper angles{result.theta};
0143     calcAngularDerivatives(angles, fitPars, thetaPrime, thetaTwoPrime);
0144     if (std::abs(thetaTwoPrime) < 2. * std::numeric_limits<double>::epsilon()) {
0145       ACTS_WARNING(__func__
0146                    << "() - " << __LINE__
0147                    << ": The fast straw fit encountered a singular "
0148                       "second derivative "
0149                    << fitPars << ", \n"
0150                    << printThetaStep(result.theta, thetaPrime, thetaTwoPrime)
0151                    << "\n"
0152                    << result);
0153       return std::nullopt;
0154     }
0155     const double update = thetaPrime / thetaTwoPrime;
0156     ACTS_VERBOSE(
0157         __func__ << "() - " << __LINE__ << ": Fit iteration #" << result.nIter
0158                  << " "
0159                  << printThetaStep(result.theta, thetaPrime, thetaTwoPrime));
0160 
0161     if (std::abs(update) < m_cfg.precCutOff) {
0162       completeResult(fitPars, thetaTwoPrime, result);
0163       return result;
0164     }
0165     result.theta -= update;
0166   }
0167 
0168   ACTS_WARNING(
0169       __func__ << "() - " << __LINE__
0170                << ": The fast straw fit did not converge " << fitPars << ", \n"
0171                << printThetaStep(result.theta, thetaPrime, thetaTwoPrime)
0172                << "\n"
0173                << result);
0174   return std::nullopt;
0175 }
0176 
0177 FastStrawLineFitter::UpdateStatus FastStrawLineFitter::updateIteration(
0178     const FitAuxiliariesWithT0& fitPars, FitResultT0& fitResult) const {
0179   ++fitResult.nIter;
0180   if (fitResult.nIter > m_cfg.maxIter) {
0181     ACTS_WARNING(__func__ << "() - " << __LINE__
0182                           << ": The fast straw fit did not converge " << fitPars
0183                           << "\n"
0184                           << fitResult);
0185     return UpdateStatus::Exceeded;
0186   }
0187 
0188   ActsSquareMatrix<2> cov{ActsSquareMatrix<2>::Zero()};
0189   Vector2 grad{Vector2::Zero()};
0190 
0191   const TrigonomHelper angles{fitResult.theta};
0192 
0193   calcAngularDerivatives(angles, fitPars, grad[0], cov(0, 0));
0194   grad[1] = calcTimeGrad(angles, fitPars);
0195   if (grad.norm() < m_cfg.precCutOff) {
0196     completeResult(fitPars, cov(0, 0), fitResult);
0197     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Fit converged "
0198                         << fitResult);
0199     return UpdateStatus::Converged;
0200   }
0201 
0202   cov(1, 0) = cov(0, 1) =
0203       fitPars.T_vz * angles.cosTheta + fitPars.T_vy * angles.sinTheta;
0204 
0205   cov(1, 1) = -fitPars.R_v * fitPars.R_v * fitPars.covNorm -
0206               fitPars.fitY0 * fitPars.R_a + fitPars.R_vv + fitPars.R_ar -
0207               (fitPars.T_az * angles.sinTheta - fitPars.T_ay * angles.cosTheta);
0208   ACTS_VERBOSE(
0209       __func__
0210       << "() - " << __LINE__
0211       << ": -fitPars.R_v * fitPars.R_v * fitPars.covNorm + fitPars.R_vv: "
0212       << (-fitPars.R_v * fitPars.R_v * fitPars.covNorm + fitPars.R_vv)
0213       << ", fitPars.fitY0 * fitPars.R_a: " << (fitPars.fitY0 * fitPars.R_a)
0214       << ", fitPars.R_ar: " << fitPars.R_ar
0215       << ", fitPars.T_az * angles.sinTheta - fitPars.T_ay * angles.cosTheta: "
0216       << (fitPars.T_az * angles.sinTheta - fitPars.T_ay * angles.cosTheta));
0217   const auto invCov = cov.inverse();
0218   if (invCov(1, 1) < 0) {
0219     ACTS_WARNING("Invalid covariance\n"
0220                  << cov << ", determinant: " << cov.determinant() << ", "
0221                  << fitPars);
0222     return UpdateStatus::Exceeded;
0223   }
0224   const Vector2 update = invCov * grad;
0225 
0226   ACTS_VERBOSE(__func__ << "() - " << __LINE__ << " intermediate result "
0227                         << fitResult << "\n"
0228                         << std::format("gradient: ({:.3f}, {:.3f})",
0229                                        inDeg(grad[0]), inNanoS(grad[1]))
0230                         << ", covariance:" << std::endl
0231                         << toString(cov) << std::endl
0232                         << cov.determinant()
0233                         << std::format(" update: ({:.3f}, {:.3f}).",
0234                                        inDeg(update[0]), inNanoS(update[1])));
0235 
0236   if (update.norm() < m_cfg.precCutOff) {
0237     completeResult(fitPars, cov(0, 0), fitResult);
0238     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Fit converged "
0239                         << fitResult);
0240     return UpdateStatus::Converged;
0241   }
0242   fitResult.t0 -= update[1];
0243   fitResult.theta -= update[0];
0244   fitResult.dT0 = std::sqrt(1. / cov(1, 1));
0245   completeResult(fitPars, cov(0, 0), fitResult);
0246 
0247   return UpdateStatus::GoodStep;
0248 }
0249 }  // namespace Acts::Experimental::detail
[ 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!