Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-07-16 08:07:35 

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/detail/FastStrawLineFitter.hpp"
0012 
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Utilities/Enumerate.hpp"
0015 
0016 namespace Acts::Experimental::detail {
0017 
0018 template <CompositeSpacePointContainer StrawCont_t>
0019 std::optional<FastStrawLineFitter::FitResult> FastStrawLineFitter::fit(
0020     const StrawCont_t& measurements, const std::vector<int>& signs) const {
0021   if (measurements.size() != signs.size()) {
0022     ACTS_WARNING(
0023         __func__ << "() - " << __LINE__
0024                  << ": Not all measurements are associated with a drift sign");
0025     return std::nullopt;
0026   }
0027 
0028   auto result = fit(fillAuxiliaries(measurements, signs));
0029   if (!result) {
0030     return std::nullopt;
0031   }
0032   // Calculate the chi2
0033   calcPostFitChi2(measurements, *result);
0034   return result;
0035 }
0036 
0037 template <CompositeSpacePointContainer StrawCont_t>
0038 void FastStrawLineFitter::calcPostFitChi2(const StrawCont_t& measurements,
0039                                           FitResult& result) const {
0040   const TrigonomHelper angles{result.theta};
0041   result.chi2 = 0.;
0042   for (const auto& strawMeas : measurements) {
0043     result.chi2 += chi2Term(angles, result.y0, *strawMeas);
0044   }
0045   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Overall chi2: "
0046                       << result.chi2 << ", nDoF: " << result.nDoF
0047                       << ", redChi2: " << (result.chi2 / result.nDoF));
0048 }
0049 
0050 template <CompositeSpacePoint Point_t>
0051 double FastStrawLineFitter::chi2Term(const TrigonomHelper& angle,
0052                                      const double y0, const Point_t& strawMeas,
0053                                      std::optional<double> r) const {
0054   if (!strawMeas.isStraw()) {
0055     return 0.;
0056   }
0057   const double cov = strawMeas.covariance()[s_covIdx];
0058   if (cov < std::numeric_limits<double>::epsilon()) {
0059     return 0.;
0060   }
0061   const Vector& pos = strawMeas.localPosition();
0062   const double y = pos.dot(strawMeas.toNextSensor());
0063   const double z = pos.dot(strawMeas.planeNormal());
0064   const double dist = Acts::abs((y - y0) * angle.cosTheta - z * angle.sinTheta);
0065   ACTS_VERBOSE(__func__ << "() - " << __LINE__ << ": Distance straw (" << y
0066                         << ", " << z
0067                         << "), r: " << r.value_or(strawMeas.driftRadius())
0068                         << " - track: " << dist);
0069   return Acts::pow(dist - r.value_or(strawMeas.driftRadius()), 2) / cov;
0070 }
0071 
0072 template <CompositeSpacePointContainer StripCont_t>
0073 std::optional<FastStrawLineFitter::FitResult> FastStrawLineFitter::fit(
0074     const StripCont_t& measurements, const ResidualIdx projection) const {
0075   if (projection == ResidualIdx::time) {
0076     ACTS_WARNING(__func__ << "() - " << __LINE__
0077                           << ": Only spatial projections, "
0078                           << "i.e. nonBending / bending are sensible");
0079     return std::nullopt;
0080   }
0081 
0082   FitAuxiliaries auxVars{};
0083   Vector centerOfGravity{Vector::Zero()};
0084   // @brief Selector function to check that the current strip provides a constraint in the projetor direction
0085   auto select = [&projection](const auto& strip) -> bool {
0086     // Skip straw measurements that are not twins &
0087     // don't measure non-bending coordinate
0088     if (strip->isStraw()) {
0089       return strip->measuresLoc0() && projection == ResidualIdx::nonBending;
0090     }
0091     // Check that the strip is actually measuring the projection
0092     return (strip->measuresLoc0() && projection == ResidualIdx::nonBending) ||
0093            (strip->measuresLoc1() && projection == ResidualIdx::bending);
0094   };
0095 
0096   auxVars.invCovs.resize(measurements.size());
0097   for (const auto& [sIdx, strip] : enumerate(measurements)) {
0098     if (!select(strip)) {
0099       ACTS_VERBOSE(__func__ << "() - " << __LINE__
0100                             << ": Skip strip measurement " << toString(*strip));
0101       continue;
0102     }
0103     const auto& invCov =
0104         (auxVars.invCovs[sIdx] =
0105              1. / strip->covariance()[toUnderlying(projection)]);
0106     auxVars.covNorm += invCov;
0107     centerOfGravity += invCov * strip->localPosition();
0108     ++auxVars.nDoF;
0109   }
0110   // To little information provided
0111   if (auxVars.nDoF < 3) {
0112     return std::nullopt;
0113   }
0114   // Reduce the number of degrees of freedom by 2 to account
0115   // for the two free parameters to fit
0116   auxVars.nDoF -= 2u;
0117   auxVars.covNorm = 1. / auxVars.covNorm;
0118   centerOfGravity *= auxVars.covNorm;
0119 
0120   bool centerSet{false};
0121   for (const auto& [sIdx, strip] : enumerate(measurements)) {
0122     if (!select(strip)) {
0123       continue;
0124     }
0125     const Vector pos = strip->localPosition() - centerOfGravity;
0126     const Vector& measDir{
0127         (projection == ResidualIdx::nonBending && strip->measuresLoc1()) ||
0128                 strip->isStraw()
0129             ? strip->sensorDirection()
0130             : strip->toNextSensor()};
0131 
0132     if (!centerSet) {
0133       auxVars.centerY = centerOfGravity.dot(measDir);
0134       auxVars.centerZ = centerOfGravity.dot(strip->planeNormal());
0135       centerSet = true;
0136     }
0137 
0138     const double y = pos.dot(measDir);
0139     const double z = pos.dot(strip->planeNormal());
0140 
0141     const auto& invCov = auxVars.invCovs[sIdx];
0142     auxVars.T_zzyy += invCov * (Acts::square(z) - Acts::square(y));
0143     auxVars.T_yz += invCov * z * y;
0144   }
0145   return fit(auxVars);
0146 }
0147 
0148 template <CompositeSpacePointContainer StrawCont_t,
0149           CompositeSpacePointFastCalibrator<
0150               Acts::RemovePointer_t<typename StrawCont_t::value_type>>
0151               Calibrator_t>
0152 void FastStrawLineFitter::calcPostFitChi2(const Acts::CalibrationContext& ctx,
0153                                           const StrawCont_t& measurements,
0154                                           const Calibrator_t& calibrator,
0155                                           FitResultT0& result) const {
0156   const TrigonomHelper angles{result.theta};
0157   result.chi2 = 0.;
0158   for (const auto& strawMeas : measurements) {
0159     result.chi2 += chi2Term(angles, result.y0, *strawMeas,
0160                             calibrator.driftRadius(ctx, *strawMeas, result.t0));
0161   }
0162   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Overall chi2: "
0163                       << result.chi2 << ", nDoF: " << result.nDoF
0164                       << ", redChi2: " << (result.chi2 / result.nDoF));
0165 }
0166 
0167 template <CompositeSpacePointContainer StrawCont_t>
0168 FastStrawLineFitter::FitAuxiliaries FastStrawLineFitter::fillAuxiliaries(
0169     const StrawCont_t& measurements, const std::vector<int>& signs) const {
0170   FitAuxiliaries auxVars{};
0171   auxVars.invCovs.resize(signs.size(), -1.);
0172   Vector centerOfGravity{Vector::Zero()};
0173 
0174   // Calculate first the center of gravity
0175   for (const auto& [sIdx, strawMeas] : enumerate(measurements)) {
0176     if (!strawMeas->isStraw()) {
0177       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": The measurement "
0178                           << toString(*strawMeas) << " is not a straw");
0179       continue;
0180     }
0181     const double cov = strawMeas->covariance()[s_covIdx];
0182     if (cov < std::numeric_limits<double>::epsilon()) {
0183       ACTS_WARNING(__func__ << "() - " << __LINE__ << ": The covariance ("
0184                             << cov << ") of the measurement "
0185                             << toString(*strawMeas) << " is invalid.");
0186       continue;
0187     }
0188     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << ": Fill "
0189                           << toString(*strawMeas) << ".");
0190 
0191     auto& invCov = (auxVars.invCovs[sIdx] = 1. / cov);
0192     auxVars.covNorm += invCov;
0193     centerOfGravity += invCov * strawMeas->localPosition();
0194     ++auxVars.nDoF;
0195   }
0196   if (auxVars.nDoF < 3) {
0197     std::stringstream sstr{};
0198     for (const auto& [sIdx, strawMeas] : enumerate(measurements)) {
0199       sstr << " --- " << (sIdx + 1) << ") " << toString(*strawMeas)
0200            << ", weight: " << auxVars.invCovs[sIdx] << std::endl;
0201     }
0202     ACTS_WARNING(__func__ << "() - " << __LINE__
0203                           << ": At least 3 measurements are required to "
0204                              "perform the straw line fit\n"
0205                           << sstr.str());
0206     auxVars.nDoF = 0u;
0207     return auxVars;
0208   }
0209   // Reduce the number of degrees of freedom by 2 to account
0210   // for the two free parameters to fit
0211   auxVars.nDoF -= 2u;
0212   auxVars.covNorm = 1. / auxVars.covNorm;
0213   centerOfGravity *= auxVars.covNorm;
0214 
0215   // Now calculate the fit constants
0216   bool centerSet{false};
0217   for (const auto& [sIdx, strawMeas] : enumerate(measurements)) {
0218     const auto& invCov = auxVars.invCovs[sIdx];
0219     // Invalid measurements were marked
0220     if (invCov < 0.) {
0221       continue;
0222     }
0223     if (!centerSet) {
0224       auxVars.centerY = centerOfGravity.dot(strawMeas->toNextSensor());
0225       auxVars.centerZ = centerOfGravity.dot(strawMeas->planeNormal());
0226       centerSet = true;
0227     }
0228     const Vector pos = strawMeas->localPosition() - centerOfGravity;
0229     const double y = pos.dot(strawMeas->toNextSensor());
0230     const double z = pos.dot(strawMeas->planeNormal());
0231     const double r = strawMeas->driftRadius();
0232 
0233     auxVars.T_zzyy += invCov * (Acts::square(z) - Acts::square(y));
0234     auxVars.T_yz += invCov * z * y;
0235     const double sInvCov = -invCov * signs[sIdx];
0236     auxVars.T_rz += sInvCov * z * r;
0237     auxVars.T_ry += sInvCov * y * r;
0238     auxVars.fitY0 += sInvCov * r;
0239   }
0240   auxVars.fitY0 *= auxVars.covNorm;
0241 
0242   return auxVars;
0243 }
0244 
0245 template <CompositeSpacePointContainer StrawCont_t,
0246           CompositeSpacePointFastCalibrator<
0247               Acts::RemovePointer_t<typename StrawCont_t::value_type>>
0248               Calibrator_t>
0249 std::optional<FastStrawLineFitter::FitResultT0> FastStrawLineFitter::fit(
0250     const Acts::CalibrationContext& ctx, const Calibrator_t& calibrator,
0251     const StrawCont_t& measurements, const std::vector<int>& signs,
0252     std::optional<double> startT0) const {
0253   using namespace Acts::UnitLiterals;
0254   if (measurements.size() != signs.size()) {
0255     ACTS_WARNING(
0256         __func__ << "() - " << __LINE__
0257                  << ": Not all measurements are associated with a drift sign");
0258     return std::nullopt;
0259   }
0260 
0261   FitResultT0 result{};
0262   result.t0 = startT0.value_or(0.);
0263 
0264   FitAuxiliariesWithT0 fitPars{
0265       fillAuxiliaries(ctx, calibrator, measurements, signs, result.t0)};
0266   result.theta = startTheta(fitPars);
0267   result.nDoF = fitPars.nDoF;
0268   ACTS_DEBUG(__func__ << "() - " << __LINE__
0269                       << ": Initial fit parameters: " << result);
0270   UpdateStatus iterStatus{UpdateStatus::GoodStep};
0271 
0272   while ((iterStatus = updateIteration(fitPars, result)) !=
0273          UpdateStatus::Exceeded) {
0274     if (iterStatus == UpdateStatus::Converged) {
0275       calcPostFitChi2(ctx, measurements, calibrator, result);
0276       return result;
0277     }
0278     fitPars = fillAuxiliaries(ctx, calibrator, measurements, signs, result.t0);
0279   }
0280   if (logger().doPrint(Logging::VERBOSE)) {
0281     ACTS_VERBOSE("Fit failed, printing all measurements:");
0282     for (const auto& meas : measurements) {
0283       ACTS_VERBOSE(toString(*meas)
0284                    << ", t0: " << result.t0 / 1._ns
0285                    << ", truthR, RecoR: " << meas->driftRadius() << ", "
0286                    << calibrator.driftRadius(ctx, *meas, result.t0)
0287                    << ", velocity: "
0288                    << calibrator.driftVelocity(ctx, *meas, result.t0) * 1._ns
0289                    << ", acceleration: "
0290                    << calibrator.driftAcceleration(ctx, *meas, result.t0) *
0291                           Acts::square(1._ns));
0292     }
0293     ACTS_VERBOSE("Result: " << result);
0294   }
0295   return std::nullopt;
0296 }
0297 
0298 template <CompositeSpacePointContainer StrawCont_t,
0299           CompositeSpacePointFastCalibrator<
0300               Acts::RemovePointer_t<typename StrawCont_t::value_type>>
0301               Calibrator_t>
0302 FastStrawLineFitter::FitAuxiliariesWithT0 FastStrawLineFitter::fillAuxiliaries(
0303     const CalibrationContext& ctx, const Calibrator_t& calibrator,
0304     const StrawCont_t& measurements, const std::vector<int>& signs,
0305     const double t0) const {
0306   using namespace Acts::UnitLiterals;
0307   FitAuxiliariesWithT0 auxVars{fillAuxiliaries(measurements, signs)};
0308   if (auxVars.nDoF <= 1) {
0309     auxVars.nDoF = 0;
0310     return auxVars;
0311   }
0312   // Account for the time offset as extra degree of freedom
0313   --auxVars.nDoF;
0314   // Fill the new extra variables
0315   auxVars.T_rz = 0.;
0316   auxVars.T_ry = 0.;
0317   auxVars.fitY0 = 0.;
0318   for (const auto& [spIdx, strawMeas] : enumerate(measurements)) {
0319     const double& invCov = auxVars.invCovs[spIdx];
0320     // Invalid (non)-straw measurements
0321     if (invCov < 0.) {
0322       continue;
0323     }
0324     const double sInvCov = -invCov * signs[spIdx];
0325     const double r = calibrator.driftRadius(ctx, *strawMeas, t0);
0326     const double v = calibrator.driftVelocity(ctx, *strawMeas, t0);
0327     const double a = calibrator.driftAcceleration(ctx, *strawMeas, t0);
0328     const double y = strawMeas->localPosition().dot(strawMeas->toNextSensor()) -
0329                      auxVars.centerY;
0330     const double z = strawMeas->localPosition().dot(strawMeas->planeNormal()) -
0331                      auxVars.centerZ;
0332 
0333     ACTS_VERBOSE(__func__ << "() - " << __LINE__ << ": # " << (spIdx + 1)
0334                           << ") " << toString(*strawMeas) << ", t0: "
0335                           << t0 / 1._ns << " r: " << r << ", v: " << v * 1._ns
0336                           << ", a: " << a * Acts::square(1._ns));
0337     auxVars.fitY0 += sInvCov * r;
0338     auxVars.R_v += sInvCov * v;
0339     auxVars.R_a += sInvCov * a;
0340 
0341     auxVars.T_rz += sInvCov * z * r;
0342     auxVars.T_ry += sInvCov * y * r;
0343 
0344     auxVars.T_vy += sInvCov * v * y;
0345     auxVars.T_vz += sInvCov * v * z;
0346 
0347     auxVars.R_vr += invCov * r * v;
0348     auxVars.R_vv += invCov * v * v;
0349 
0350     auxVars.T_ay += sInvCov * a * y;
0351     auxVars.T_az += sInvCov * a * z;
0352 
0353     auxVars.R_ar += invCov * a * r;
0354   }
0355   auxVars.fitY0 *= auxVars.covNorm;
0356   ACTS_DEBUG(__func__ << "() - " << __LINE__ << " Fit constants calculated \n"
0357                       << auxVars);
0358   return auxVars;
0359 }
0360 
0361 }  // 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!