Seeding/detail/CompSpacePointAuxiliaries.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 "Acts/Seeding/detail/CompSpacePointAuxiliaries.hpp"
0012
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Surfaces/detail/LineHelper.hpp"
0015 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0016 #include "Acts/Utilities/StringHelpers.hpp"
0017
0018 #include <format>
0019
0020 namespace Acts {
0021 template <Experimental::CompositeSpacePoint SpacePoint_t>
0022 std::string toString(const SpacePoint_t& measurement) {
0023   if constexpr (hasPrintOperator<SpacePoint_t>) {
0024     std::ostringstream sstr{};
0025     sstr << measurement;
0026     return sstr.str();
0027   } else {
0028     using ResidualIdx =
0029         Experimental::detail::CompSpacePointAuxiliaries::ResidualIdx;
0030     if (measurement.isStraw()) {
0031       return std::format(
0032           "straw SP @ {:} with r: {:.3f}+-{:.3f} & wire : {:} ",
0033           toString(measurement.localPosition()), measurement.driftRadius(),
0034           std::sqrt(
0035               measurement.covariance()[toUnderlying(ResidualIdx::bending)]),
0036           toString(measurement.sensorDirection()));
0037     } else {
0038       return std::format(
0039           "strip SP @ {:} with normal: {:}, strip dir: {:}, to next {:}, "
0040           "measures loc0/loc1/time: {:}/{:}/{:}",
0041           toString(measurement.localPosition()),
0042           toString(measurement.planeNormal()),
0043           toString(measurement.sensorDirection()),
0044           toString(measurement.toNextSensor()),
0045           measurement.measuresLoc0() ? "yes" : "no",
0046           measurement.measuresLoc1() ? "yes" : "no",
0047           measurement.hasTime() ? "yes" : "no");
0048     }
0049   }
0050 }
0051 }  // namespace Acts
0052
0053 namespace Acts::Experimental::detail {
0054
0055 template <CompositeSpacePoint SpacePoint_t>
0056 CompSpacePointAuxiliaries::Vector CompSpacePointAuxiliaries::extrapolateToPlane(
0057     const Line_t& line, const SpacePoint_t& hit) {
0058   return extrapolateToPlane(line.position(), line.direction(), hit);
0059 }
0060 template <CompositeSpacePoint SpacePoint_t>
0061 CompSpacePointAuxiliaries::Vector CompSpacePointAuxiliaries::extrapolateToPlane(
0062     const Vector& pos, const Vector& dir, const SpacePoint_t& hit) {
0063   using namespace Acts::PlanarHelper;
0064   const auto planeIsect =
0065       intersectPlane(pos, dir, hit.planeNormal(), hit.localPosition());
0066   return planeIsect.position();
0067 }
0068 template <CompositeSpacePoint SpacePoint_t>
0069 double CompSpacePointAuxiliaries::chi2Term(const Line_t& line,
0070                                            const SpacePoint_t& hit) {
0071   return chi2Term(line.position(), line.direction(), hit);
0072 }
0073 template <CompositeSpacePoint SpacePoint_t>
0074 double CompSpacePointAuxiliaries::chi2Term(const Vector& pos, const Vector& dir,
0075                                            const SpacePoint_t& hit) {
0076   double chiSq{0.};
0077   using namespace Acts::detail::LineHelper;
0078   constexpr auto bendIdx = toUnderlying(ResidualIdx::bending);
0079   constexpr auto nonBendIdx = toUnderlying(ResidualIdx::nonBending);
0080   if (hit.isStraw()) {
0081     const double dist = Acts::abs(
0082         signedDistance(pos, dir, hit.localPosition(), hit.sensorDirection()));
0083     chiSq = Acts::square(dist - hit.driftRadius()) / hit.covariance()[bendIdx];
0084     if (hit.measuresLoc0()) {
0085       auto closePointOnStraw =
0086           lineIntersect(pos, dir, hit.localPosition(), hit.sensorDirection());
0087       chiSq += Acts::square(closePointOnStraw.pathLength()) /
0088                hit.covariance()[nonBendIdx];
0089     }
0090
0091   } else {
0092     const Vector distOnPlane =
0093         (extrapolateToPlane(pos, dir, hit) - hit.localPosition());
0094     const Vector& b1{hit.toNextSensor()};
0095     const Vector& b2{hit.sensorDirection()};
0096     Vector2 dist{distOnPlane.dot(b1), distOnPlane.dot(b2)};
0097     if (hit.measuresLoc0() && hit.measuresLoc1()) {
0098       /// Check whether the two vectors are orthogonal
0099       constexpr double tolerance = 1.e-8;
0100       const double stripAngle = b1.dot(b2);
0101       if (stripAngle > tolerance) {
0102         const double invDist = 1. / (1. - square(stripAngle));
0103         ActsSquareMatrix<2> stereoDecomp{invDist *
0104                                          ActsSquareMatrix<2>::Identity()};
0105         stereoDecomp(1, 0) = stereoDecomp(0, 1) = -stripAngle * invDist;
0106         dist = stereoDecomp * dist;
0107       }
0108       chiSq = Acts::square(dist[0]) / hit.covariance()[bendIdx] +
0109               Acts::square(dist[1]) / hit.covariance()[nonBendIdx];
0110     } else if (hit.measuresLoc0()) {
0111       chiSq = Acts::square(dist[0]) / hit.covariance()[nonBendIdx];
0112     } else {
0113       chiSq = Acts::square(dist[0]) / hit.covariance()[bendIdx];
0114       if (hit.covariance()[nonBendIdx] >
0115           std::numeric_limits<double>::epsilon()) {
0116         chiSq += square(dist[1]) / hit.covariance()[nonBendIdx];
0117       }
0118     }
0119   }
0120   return chiSq;
0121 }
0122 template <CompositeSpacePoint SpacePoint_t>
0123 double CompSpacePointAuxiliaries::chi2Term(
0124     const Line_t& line, const Acts::Transform3& localToGlobal, const double t0,
0125     const SpacePoint_t& hit) {
0126   return chi2Term(line.position(), line.direction(), localToGlobal, t0, hit);
0127 }
0128 template <CompositeSpacePoint SpacePoint_t>
0129 double CompSpacePointAuxiliaries::chi2Term(const Vector& pos, const Vector& dir,
0130                                            const Transform3& localToGlobal,
0131                                            const double t0,
0132                                            const SpacePoint_t& hit) {
0133   using namespace Acts::UnitLiterals;
0134   return chi2Term(
0135       pos, dir,
0136       t0 + (localToGlobal * extrapolateToPlane(pos, dir, hit)).norm() /
0137                PhysicalConstants::c,
0138       hit);
0139 }
0140
0141 template <CompositeSpacePoint SpacePoint_t>
0142 double CompSpacePointAuxiliaries::chi2Term(const Line_t& line, const double t0,
0143                                            const SpacePoint_t& hit) {
0144   return chi2Term(line.position(), line.direction(), t0, hit);
0145 }
0146 template <CompositeSpacePoint SpacePoint_t>
0147 double CompSpacePointAuxiliaries::chi2Term(const Vector& pos, const Vector& dir,
0148                                            const double t0,
0149                                            const SpacePoint_t& hit) {
0150   double chiSq = chi2Term(pos, dir, hit);
0151   if (!hit.hasTime() || hit.isStraw()) {
0152     return chiSq;
0153   }
0154   chiSq += Acts::square(hit.time() - t0) /
0155            hit.covariance()[toUnderlying(ResidualIdx::time)];
0156   return chiSq;
0157 }
0158
0159 template <CompositeSpacePoint Point_t>
0160 int CompSpacePointAuxiliaries::strawSign(const Line_t& line,
0161                                          const Point_t& strawSp) {
0162   return strawSign(line.position(), line.direction(), strawSp);
0163 }
0164
0165 template <CompositeSpacePoint Point_t>
0166 int CompSpacePointAuxiliaries::strawSign(const Vector& pos, const Vector& dir,
0167                                          const Point_t& strawSp) {
0168   if (!strawSp.isStraw()) {
0169     return 0;
0170   }
0171   const double dist = Acts::detail::LineHelper::signedDistance(
0172       pos, dir, strawSp.localPosition(), strawSp.sensorDirection());
0173   return dist > 0. ? 1 : -1;
0174 }
0175 template <CompositeSpacePointContainer StrawCont_t>
0176 std::vector<int> CompSpacePointAuxiliaries::strawSigns(
0177     const Line_t& line, const StrawCont_t& measurements) {
0178   return strawSigns(line.position(), line.direction(), measurements);
0179 }
0180 template <CompositeSpacePointContainer StrawCont_t>
0181 std::vector<int> CompSpacePointAuxiliaries::strawSigns(
0182     const Vector& pos, const Vector& dir, const StrawCont_t& measurements) {
0183   std::vector<int> signs{};
0184   signs.reserve(measurements.size());
0185   for (const auto& strawSp : measurements) {
0186     signs.push_back(strawSign(pos, dir, *strawSp));
0187   }
0188   return signs;
0189 }
0190 template <CompositeSpacePoint Point_t>
0191 void CompSpacePointAuxiliaries::updateSpatialResidual(
0192     const Line_t& line, const Point_t& spacePoint) {
0193   ACTS_DEBUG(__func__ << "() " << __LINE__ << " Update residual of "
0194                       << toString(line.position()) << " + "
0195                       << toString(line.direction()) << " w.r.t\n"
0196                       << toString(spacePoint));
0197   if (spacePoint.isStraw()) {
0198     /// Fetch the hit position & direction
0199     const auto& wireDir{spacePoint.sensorDirection()};
0200     /// Calculate the distance from the two reference points
0201     const Vector hitMinSeg = spacePoint.localPosition() - line.position();
0202
0203     if (!updateStrawResidual(line, hitMinSeg, wireDir,
0204                              spacePoint.driftRadius())) {
0205       return;
0206     }
0207
0208     if (m_cfg.calcAlongStraw && spacePoint.measuresLoc0()) {
0209       /// If the tube is a twin-tube, the hit position is no longer arbitrary
0210       /// along the wire. Calculate the distance along the wire towards the
0211       /// point of closest approach.
0212       updateAlongTheStraw(line, hitMinSeg, wireDir);
0213     }
0214   } else {
0215     updateStripResidual(line, spacePoint.planeNormal(),
0216                         spacePoint.toNextSensor(), spacePoint.sensorDirection(),
0217                         spacePoint.localPosition(), spacePoint.measuresLoc1(),
0218                         spacePoint.measuresLoc0());
0219   }
0220 }
0221
0222 template <CompositeSpacePoint Point_t>
0223 void CompSpacePointAuxiliaries::updateFullResidual(const Line_t& line,
0224                                                    const double timeOffset,
0225                                                    const Point_t& spacePoint,
0226                                                    const double driftV,
0227                                                    const double driftA) {
0228   /// Calculate first the spatial residual
0229   updateSpatialResidual(line, spacePoint);
0230
0231   /// Calculate the time residual for strip-like measurements
0232   if (!spacePoint.isStraw()) {
0233     /// If the measurement does not provide time, then simply reset the time
0234     /// partial components
0235     if (!spacePoint.hasTime()) {
0236       resetTime();
0237       return;
0238     }
0239     updateTimeStripRes(spacePoint.toNextSensor(), spacePoint.sensorDirection(),
0240                        spacePoint.localPosition(), spacePoint.measuresLoc1(),
0241                        spacePoint.time(), timeOffset);
0242   } else {
0243     updateTimeStrawRes(line, spacePoint.localPosition() - line.position(),
0244                        spacePoint.sensorDirection(), spacePoint.driftRadius(),
0245                        driftV, driftA);
0246   }
0247 }
0248
0249 }  // namespace Acts::Experimental::detail