Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:09:43

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023-2024 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 http://mozilla.org/MPL/2.0/.
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& /*fieldCache*/) const {
0021   // Create propagator options
0022   PropagatorOptions<> pOptions(gctx, mctx);
0023 
0024   // Length scale at which we consider to be sufficiently close to the Perigee
0025   // surface to skip the propagation.
0026   pOptions.surfaceTolerance = m_cfg.targetTolerance;
0027 
0028   // Get intersection of the track with the Perigee if the particle would
0029   // move on a straight line.
0030   // This allows us to determine whether we need to propagate the track
0031   // forward or backward to arrive at the PCA.
0032   auto intersection = perigeeSurface
0033                           .intersect(gctx, params.position(gctx),
0034                                      params.direction(), BoundaryCheck(false))
0035                           .closest();
0036 
0037   // Setting the propagation direction using the intersection length from
0038   // above.
0039   // We handle zero path length as forward propagation, but we could actually
0040   // skip the whole propagation in this case.
0041   pOptions.direction =
0042       Direction::fromScalarZeroAsPositive(intersection.pathLength());
0043 
0044   // Propagate to the PCA of the reference point
0045   auto result =
0046       m_cfg.propagator->propagateToSurface(params, perigeeSurface, pOptions);
0047   if (!result.ok()) {
0048     return result.error();
0049   }
0050 
0051   // Extracting the Perigee representation of the track wrt the reference point
0052   auto endParams = *result;
0053   BoundVector perigeeParams = endParams.parameters();
0054 
0055   // Covariance and weight matrix at the PCA to the reference point
0056   BoundSquareMatrix parCovarianceAtPCA = endParams.covariance().value();
0057   BoundSquareMatrix weightAtPCA = parCovarianceAtPCA.inverse();
0058 
0059   // Vector containing the track parameters at the PCA
0060   // Note that we parametrize the track using the following parameters:
0061   // (x, y, z, t, phi, theta, q/p),
0062   // where
0063   // -) (x, y, z, t) is the global 4D position of the PCA
0064   // -) phi and theta are the global angles of the momentum at the PCA
0065   // -) q/p is the charge divided by the total momentum at the PCA
0066   Acts::ActsVector<eLinSize> paramVec;
0067 
0068   // 4D PCA and the momentum of the track at the PCA
0069   // These quantities will be used in the computation of the constant term in
0070   // the Taylor expansion
0071   Vector4 pca;
0072   Vector3 momentumAtPCA;
0073 
0074   // Fill "paramVec", "pca", and "momentumAtPCA"
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   // Complete Jacobian (consists of positionJacobian and momentumJacobian)
0088   ActsMatrix<eBoundSize, eLinSize> completeJacobian =
0089       ActsMatrix<eBoundSize, eLinSize>::Zero(eBoundSize, eLinSize);
0090 
0091   // Perigee parameters wrt the reference point after wiggling
0092   BoundVector newPerigeeParams;
0093 
0094   // Check if wiggled angle theta are within definition range [0, pi]
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   // Wiggling each of the parameters at the PCA and computing the Perigee
0102   // parametrization of the resulting new track. This allows us to approximate
0103   // the numerical derivatives.
0104   for (unsigned int i = 0; i < eLinSize; i++) {
0105     Acts::ActsVector<eLinSize> paramVecCopy = paramVec;
0106     // Wiggle
0107     paramVecCopy(i) += m_cfg.delta;
0108 
0109     // Create curvilinear track object from our parameters. This is needed for
0110     // the propagation. Note that we work without covariance since we don't need
0111     // it to compute the derivative.
0112     Vector3 wiggledDir = makeDirectionFromPhiTheta(paramVecCopy(eLinPhi),
0113                                                    paramVecCopy(eLinTheta));
0114     // Since we work in 4D we have eLinPosSize = 4
0115     CurvilinearTrackParameters wiggledCurvilinearParams(
0116         paramVecCopy.template head<eLinPosSize>(), wiggledDir,
0117         paramVecCopy(eLinQOverP), std::nullopt, ParticleHypothesis::pion());
0118 
0119     // Obtain propagation direction
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     // Propagate to the new PCA and extract Perigee parameters
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     // Computing the numerical derivatives and filling the Jacobian
0136     completeJacobian.array().col(i) =
0137         (newPerigeeParams - perigeeParams) / m_cfg.delta;
0138     // We need to account for the periodicity of phi. We overwrite the
0139     // previously computed value for better readability.
0140     completeJacobian(eLinPhi, i) =
0141         Acts::detail::difference_periodic(newPerigeeParams(eLinPhi),
0142                                           perigeeParams(eLinPhi), 2 * M_PI) /
0143         m_cfg.delta;
0144   }
0145 
0146   // Extracting positionJacobian and momentumJacobian from the complete Jacobian
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   // Constant term of Taylor expansion (Eq. 5.38 in Ref. (1))
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 }