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