File indexing completed on 2026-03-28 07:45:41
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Utilities/MathHelpers.hpp"
0013
0014 Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0,
0015 const Vector3& sp1,
0016 const Vector3& sp2,
0017 const Vector3& bField) {
0018 return estimateTrackParamsFromSeed(sp0, 0, sp1, sp2, bField);
0019 }
0020
0021 Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0,
0022 const double t0,
0023 const Vector3& sp1,
0024 const Vector3& sp2,
0025 const Vector3& bField) {
0026
0027
0028
0029
0030
0031 const Vector3 relVec = sp1 - sp0;
0032 const Vector3 newZAxis = bField.normalized();
0033 const Vector3 newYAxis = newZAxis.cross(relVec).normalized();
0034 const Vector3 newXAxis = newYAxis.cross(newZAxis);
0035 RotationMatrix3 rotation;
0036 rotation.col(0) = newXAxis;
0037 rotation.col(1) = newYAxis;
0038 rotation.col(2) = newZAxis;
0039
0040 const Translation3 trans(sp0);
0041
0042 const Transform3 transform(trans * rotation);
0043
0044
0045 const Vector3 local1 = transform.inverse() * sp1;
0046 const Vector3 local2 = transform.inverse() * sp2;
0047
0048
0049 const Vector2 uv1 = local1.head<2>() / local1.head<2>().squaredNorm();
0050 const Vector2 uv2 = local2.head<2>() / local2.head<2>().squaredNorm();
0051 const Vector2 deltaUV2 = uv2 - uv1;
0052 const double A = deltaUV2.y() / deltaUV2.x();
0053 const double bOverS =
0054 (uv1.y() * uv2.x() - uv2.y() * uv1.x()) / deltaUV2.norm();
0055
0056 const double invTanTheta = local2.z() / local2.head<2>().norm();
0057 const Vector3 transDirection(1, A, fastHypot(1, A) * invTanTheta);
0058
0059
0060 const Vector3 direction = rotation * transDirection.normalized();
0061
0062
0063 FreeVector params = FreeVector::Zero();
0064
0065
0066 params.segment<3>(eFreePos0) = sp0;
0067
0068
0069 params.segment<3>(eFreeDir0) = direction;
0070
0071
0072
0073 const double qOverPt = 2 * bOverS / bField.norm();
0074
0075 params[eFreeQOverP] = qOverPt / fastHypot(1, invTanTheta);
0076
0077
0078 params[eFreeTime] = t0;
0079
0080 return params;
0081 }
0082
0083 Acts::Result<Acts::BoundVector> Acts::estimateTrackParamsFromSeed(
0084 const GeometryContext& gctx, const Surface& surface, const Vector3& sp0,
0085 const double t0, const Vector3& sp1, const Vector3& sp2,
0086 const Vector3& bField) {
0087 const FreeVector freeParams =
0088 estimateTrackParamsFromSeed(sp0, t0, sp1, sp2, bField);
0089 return transformFreeToBoundParameters(freeParams, surface, gctx);
0090 }
0091
0092 Acts::BoundMatrix Acts::estimateTrackParamCovariance(
0093 const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
0094 bool hasTime) {
0095 assert((params[eBoundTheta] > 0 && params[eBoundTheta] < std::numbers::pi) &&
0096 "Theta must be in the range (0, pi)");
0097
0098 BoundMatrix result = BoundMatrix::Zero();
0099
0100 for (std::size_t i = eBoundLoc0; i < eBoundSize; ++i) {
0101 double sigma = config.initialSigmas[i];
0102 double variance = sigma * sigma;
0103
0104 if (i == eBoundQOverP) {
0105
0106 double varianceTheta = result(eBoundTheta, eBoundTheta);
0107
0108
0109 variance += std::pow(
0110 config.initialSigmaQoverPt * std::sin(params[eBoundTheta]), 2);
0111
0112
0113 variance += std::pow(config.initialSigmaPtRel * params[eBoundQOverP], 2);
0114
0115
0116 variance +=
0117 varianceTheta *
0118 std::pow(params[eBoundQOverP] / std::tan(params[eBoundTheta]), 2);
0119 }
0120
0121 if (i == eBoundTime && !hasTime) {
0122
0123 variance *= config.noTimeVarInflation;
0124 }
0125
0126
0127 variance *= config.initialVarInflation[i];
0128
0129 result(i, i) = variance;
0130 }
0131
0132 return result;
0133 }