src/Seeding/EstimateTrackParamsFromSeed.cpp

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 #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   // Define a new coordinate frame with its origin at the bottom space point, z
0019   // axis long the magnetic field direction and y axis perpendicular to vector
0020   // from the bottom to middle space point. Hence, the projection of the middle
0021   // space point on the transverse plane will be located at the x axis of the
0022   // new frame.
0023   const Vector3 relVec = sp1 - sp0;
0024   const Vector3 newZAxis = bField.normalized();
0025   const Vector3 newYAxis = newZAxis.cross(relVec).normalized();
0026   const Vector3 newXAxis = newYAxis.cross(newZAxis);
0027   RotationMatrix3 rotation;
0028   rotation.col(0) = newXAxis;
0029   rotation.col(1) = newYAxis;
0030   rotation.col(2) = newZAxis;
0031   // The center of the new frame is at the bottom space point
0032   const Translation3 trans(sp0);
0033   // The transform which constructs the new frame
0034   const Transform3 transform(trans * rotation);
0035
0036   // The coordinate of the middle and top space point in the new frame
0037   const Vector3 local1 = transform.inverse() * sp1;
0038   const Vector3 local2 = transform.inverse() * sp2;
0039
0040   // Use the uv-plane to estimate the circle parameters
0041   const Vector2 uv1 = local1.head<2>() / local1.head<2>().squaredNorm();
0042   const Vector2 uv2 = local2.head<2>() / local2.head<2>().squaredNorm();
0043   const Vector2 deltaUV2 = uv2 - uv1;
0044   const double A = deltaUV2.y() / deltaUV2.x();
0045   const double bOverS =
0046       (uv1.y() * uv2.x() - uv2.y() * uv1.x()) / deltaUV2.norm();
0047
0048   const double invTanTheta = local2.z() / local2.head<2>().norm();
0049   const Vector3 transDirection(1, A, fastHypot(1, A) * invTanTheta);
0050
0051   // Transform it back to the original frame
0052   const Vector3 direction = rotation * transDirection.normalized();
0053
0054   // Initialize the free parameters vector
0055   FreeVector params = FreeVector::Zero();
0056
0057   // The bottom space point position
0058   params.segment<3>(eFreePos0) = sp0;
0059
0060   // The estimated direction
0061   params.segment<3>(eFreeDir0) = direction;
0062
0063   // The estimated q/pt in [GeV/c]^-1 (note that the pt is the projection of
0064   // momentum on the transverse plane of the new frame)
0065   const double qOverPt = 2 * bOverS / bField.norm();
0066   // The estimated q/p in [GeV/c]^-1
0067   params[eFreeQOverP] = qOverPt / fastHypot(1, invTanTheta);
0068
0069   return params;
0070 }
0071
0072 Acts::BoundMatrix Acts::estimateTrackParamCovariance(
0073     const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
0074     bool hasTime) {
0075   assert((params[eBoundTheta] > 0 && params[eBoundTheta] < std::numbers::pi) &&
0076          "Theta must be in the range (0, pi)");
0077
0078   BoundSquareMatrix result = BoundSquareMatrix::Zero();
0079
0080   for (std::size_t i = eBoundLoc0; i < eBoundSize; ++i) {
0081     double sigma = config.initialSigmas[i];
0082     double variance = sigma * sigma;
0083
0084     if (i == eBoundQOverP) {
0085       // note that we rely on the fact that sigma theta is already computed
0086       double varianceTheta = result(eBoundTheta, eBoundTheta);
0087
0088       // contribution from sigma(q/pt)
0089       variance += std::pow(
0090           config.initialSigmaQoverPt * std::sin(params[eBoundTheta]), 2);
0091
0092       // contribution from sigma(pt)/pt
0093       variance += std::pow(config.initialSigmaPtRel * params[eBoundQOverP], 2);
0094
0095       // contribution from sigma(theta)
0096       variance +=
0097           varianceTheta *
0098           std::pow(params[eBoundQOverP] / std::tan(params[eBoundTheta]), 2);
0099     }
0100
0101     if (i == eBoundTime && !hasTime) {
0102       // Inflate the time uncertainty if no time measurement is available
0103       variance *= config.noTimeVarInflation;
0104     }
0105
0106     // Inflate the initial covariance
0107     variance *= config.initialVarInflation[i];
0108
0109     result(i, i) = variance;
0110   }
0111
0112   return result;
0113 }