src/Seeding/FastStrawLineFitter.cpp

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 std::string printThetaStep(const double theta, const double thetaPrime,
0026                            const double thetaTwoPrime) {
0027   return std::format(
0028       "-- theta: {:.3f}, thetaPrime: {:.3f}, thetaTwoPrime: {:.3f}  --> "
0029       "step-size: {:.8f}",
0030       inDeg(theta), inDeg(thetaPrime), inDeg(thetaTwoPrime),
0031       inDeg(thetaPrime / thetaTwoPrime));
0032 }
0033 }  // namespace
0034 namespace Acts::Experimental::detail {
0035
0036 void FastStrawLineFitter::FitAuxiliaries::print(std::ostream& ostr) const {
0037   ostr << std::format("T_zzyy: {:.3f}, ", T_zzyy)
0038        << std::format("T_yz: {:.3f}, ", T_yz)
0039        << std::format("T_rz: {:.3f}, ", T_rz)
0040        << std::format("T_ry: {:.3f}, ", T_ry)
0041        << std::format("centre ( {:.3f}, {:3f}), ", centerY, centerZ)
0042        << std::format("y0: {:.3f}", fitY0)
0043        << std::format(", norm: {:.3f}", covNorm) << ", nDoF: " << nDoF;
0044 }
0045 void FastStrawLineFitter::FitResult::print(std::ostream& ostr) const {
0046   ostr << "# iteration: " << nIter << ", nDoF: " << nDoF << ", chi2: " << chi2
0047        << ", chi2 / nDoF: " << (chi2 / static_cast<double>(nDoF)) << ",\n";
0048   ostr << std::format("theta: {:.3f} pm {:.3f}, ", inDeg(theta), inDeg(dTheta))
0049        << std::format("y0: {:.3f}  pm {:.3f}", y0, dY0);
0050 }
0051
0052 using Vector = FastStrawLineFitter::Vector;
0053 const Acts::Logger& FastStrawLineFitter::logger() const {
0054   assert(m_logger != nullptr);
0055   return *m_logger;
0056 }
0057 FastStrawLineFitter::FastStrawLineFitter(const Config& cfg,
0058                                          std::unique_ptr<const Logger> logger)
0059     : m_cfg{cfg}, m_logger{std::move(logger)} {}
0060
0061 inline void FastStrawLineFitter::calcAngularDerivatives(
0062     const TrigonomHelper& angles, const FitAuxiliaries& fitPars,
0063     double& thetaPrime, double& thetaTwoPrime) const {
0064   thetaPrime = 0.5 * fitPars.T_zzyy * angles.sinTwoTheta -
0065                fitPars.T_yz * angles.cosTwoTheta -
0066                fitPars.T_rz * angles.cosTheta - fitPars.T_ry * angles.sinTheta;
0067   thetaTwoPrime = fitPars.T_zzyy * angles.cosTwoTheta +
0068                   2. * fitPars.T_yz * angles.sinTwoTheta +
0069                   fitPars.T_rz * angles.sinTheta -
0070                   fitPars.T_ry * angles.cosTheta;
0071 }
0072 double FastStrawLineFitter::startTheta(const FitAuxiliaries& fitPars) {
0073   double thetaGuess =
0074       std::atan2(2. * (fitPars.T_yz - fitPars.T_ry), fitPars.T_zzyy) / 2.;
0075   return Acts::detail::wrap_periodic(thetaGuess, 0., std::numbers::pi);
0076 }
0077 void FastStrawLineFitter::completeResult(const FitAuxiliaries& fitPars,
0078                                          const double thetaTwoPrime,
0079                                          FitResult& result) const {
0080   result.dTheta = std::sqrt(1. / thetaTwoPrime);
0081   const double tanTheta = std::tan(result.theta);
0082   const double secTheta = 1. / std::cos(result.theta);
0083   result.y0 =
0084       fitPars.centerY - fitPars.centerZ * tanTheta + fitPars.fitY0 * secTheta;
0085   result.dY0 = Acts::fastHypot(-fitPars.centerZ * Acts::square(secTheta) +
0086                                    fitPars.fitY0 * secTheta * tanTheta,
0087                                secTheta);
0088
0089   result.dY0 *= result.dTheta;
0090   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Fit succeeded " << result);
0091 }
0092
0093 std::optional<FastStrawLineFitter::FitResult> FastStrawLineFitter::fit(
0094     const FitAuxiliaries& fitPars) const {
0095   /// No degrees of freedom -> no valid parameters
0096   if (fitPars.nDoF == 0u) {
0097     return std::nullopt;
0098   }
0099   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Input fit parameters "
0100                       << fitPars);
0101
0102   const double thetaGuess = startTheta(fitPars);
0103   ACTS_DEBUG(__func__ << "() - " << __LINE__
0104                       << ": Start fast fit seed  theta: " << inDeg(thetaGuess));
0105   ////
0106   FitResult result{};
0107   result.theta = thetaGuess;
0108   result.nDoF = fitPars.nDoF;
0109
0110   double thetaPrime{0.};
0111   double thetaTwoPrime{0.};
0112   while (result.nIter <= m_cfg.maxIter) {
0113     ++result.nIter;
0114     const TrigonomHelper angles{result.theta};
0115     calcAngularDerivatives(angles, fitPars, thetaPrime, thetaTwoPrime);
0116     const double update = thetaPrime / thetaTwoPrime;
0117     ACTS_VERBOSE(
0118         __func__ << "() - " << __LINE__ << ": Fit iteration #" << result.nIter
0119                  << " "
0120                  << printThetaStep(result.theta, thetaPrime, thetaTwoPrime));
0121
0122     if (std::abs(update) < m_cfg.precCutOff) {
0123       completeResult(fitPars, thetaTwoPrime, result);
0124       return result;
0125     }
0126     result.theta -= update;
0127   }
0128
0129   ACTS_WARNING(
0130       __func__ << "() - " << __LINE__
0131                << ": The fast straw fit did not converge " << fitPars << ", \n"
0132                << printThetaStep(result.theta, thetaPrime, thetaTwoPrime)
0133                << "\n"
0134                << result);
0135   return std::nullopt;
0136 }
0137
0138 }  // namespace Acts::Experimental::detail