File indexing completed on 2025-08-05 08:09:43
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Vertexing/NumericalTrackLinearizer.hpp"
0010
0011 #include "Acts/Surfaces/PerigeeSurface.hpp"
0012 #include "Acts/Utilities/UnitVectors.hpp"
0013 #include "Acts/Vertexing/LinearizerTrackParameters.hpp"
0014
0015 Acts::Result<Acts::LinearizedTrack>
0016 Acts::NumericalTrackLinearizer::linearizeTrack(
0017 const BoundTrackParameters& params, double linPointTime,
0018 const Surface& perigeeSurface, const Acts::GeometryContext& gctx,
0019 const Acts::MagneticFieldContext& mctx,
0020 MagneticFieldProvider::Cache& ) const {
0021
0022 PropagatorOptions<> pOptions(gctx, mctx);
0023
0024
0025
0026 pOptions.surfaceTolerance = m_cfg.targetTolerance;
0027
0028
0029
0030
0031
0032 auto intersection = perigeeSurface
0033 .intersect(gctx, params.position(gctx),
0034 params.direction(), BoundaryCheck(false))
0035 .closest();
0036
0037
0038
0039
0040
0041 pOptions.direction =
0042 Direction::fromScalarZeroAsPositive(intersection.pathLength());
0043
0044
0045 auto result =
0046 m_cfg.propagator->propagateToSurface(params, perigeeSurface, pOptions);
0047 if (!result.ok()) {
0048 return result.error();
0049 }
0050
0051
0052 auto endParams = *result;
0053 BoundVector perigeeParams = endParams.parameters();
0054
0055
0056 BoundSquareMatrix parCovarianceAtPCA = endParams.covariance().value();
0057 BoundSquareMatrix weightAtPCA = parCovarianceAtPCA.inverse();
0058
0059
0060
0061
0062
0063
0064
0065
0066 Acts::ActsVector<eLinSize> paramVec;
0067
0068
0069
0070
0071 Vector4 pca;
0072 Vector3 momentumAtPCA;
0073
0074
0075 {
0076 Vector3 globalCoords = endParams.position(gctx);
0077 ActsScalar globalTime = endParams.time();
0078 ActsScalar phi = perigeeParams(BoundIndices::eBoundPhi);
0079 ActsScalar theta = perigeeParams(BoundIndices::eBoundTheta);
0080 ActsScalar qOvP = perigeeParams(BoundIndices::eBoundQOverP);
0081
0082 paramVec << globalCoords, globalTime, phi, theta, qOvP;
0083 pca << globalCoords, globalTime;
0084 momentumAtPCA << phi, theta, qOvP;
0085 }
0086
0087
0088 ActsMatrix<eBoundSize, eLinSize> completeJacobian =
0089 ActsMatrix<eBoundSize, eLinSize>::Zero(eBoundSize, eLinSize);
0090
0091
0092 BoundVector newPerigeeParams;
0093
0094
0095 if (paramVec(eLinTheta) + m_cfg.delta > M_PI) {
0096 ACTS_ERROR(
0097 "Wiggled theta outside range, choose a smaller wiggle (i.e., delta)! "
0098 "You might need to decrease targetTolerance as well.")
0099 }
0100
0101
0102
0103
0104 for (unsigned int i = 0; i < eLinSize; i++) {
0105 Acts::ActsVector<eLinSize> paramVecCopy = paramVec;
0106
0107 paramVecCopy(i) += m_cfg.delta;
0108
0109
0110
0111
0112 Vector3 wiggledDir = makeDirectionFromPhiTheta(paramVecCopy(eLinPhi),
0113 paramVecCopy(eLinTheta));
0114
0115 CurvilinearTrackParameters wiggledCurvilinearParams(
0116 paramVecCopy.template head<eLinPosSize>(), wiggledDir,
0117 paramVecCopy(eLinQOverP), std::nullopt, ParticleHypothesis::pion());
0118
0119
0120 intersection = perigeeSurface
0121 .intersect(gctx, paramVecCopy.template head<3>(),
0122 wiggledDir, BoundaryCheck(false))
0123 .closest();
0124 pOptions.direction =
0125 Direction::fromScalarZeroAsPositive(intersection.pathLength());
0126
0127
0128 auto newResult = m_cfg.propagator->propagateToSurface(
0129 wiggledCurvilinearParams, perigeeSurface, pOptions);
0130 if (!newResult.ok()) {
0131 return newResult.error();
0132 }
0133 newPerigeeParams = newResult->parameters();
0134
0135
0136 completeJacobian.array().col(i) =
0137 (newPerigeeParams - perigeeParams) / m_cfg.delta;
0138
0139
0140 completeJacobian(eLinPhi, i) =
0141 Acts::detail::difference_periodic(newPerigeeParams(eLinPhi),
0142 perigeeParams(eLinPhi), 2 * M_PI) /
0143 m_cfg.delta;
0144 }
0145
0146
0147 ActsMatrix<eBoundSize, eLinPosSize> positionJacobian =
0148 completeJacobian.block<eBoundSize, eLinPosSize>(0, 0);
0149 ActsMatrix<eBoundSize, eLinMomSize> momentumJacobian =
0150 completeJacobian.block<eBoundSize, eLinMomSize>(0, eLinPosSize);
0151
0152
0153 BoundVector constTerm =
0154 perigeeParams - positionJacobian * pca - momentumJacobian * momentumAtPCA;
0155
0156 Vector4 linPoint;
0157 linPoint.head<3>() = perigeeSurface.center(gctx);
0158 linPoint[3] = linPointTime;
0159
0160 return LinearizedTrack(perigeeParams, parCovarianceAtPCA, weightAtPCA,
0161 linPoint, positionJacobian, momentumJacobian, pca,
0162 momentumAtPCA, constTerm);
0163 }