Back to home page

EIC code displayed by LXR

 
 

    


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