Core/Vertexing/LinearizedTrackFactoryTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019-2023 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 <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/output_test_stream.hpp>
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/Algebra.hpp"
0014 #include "Acts/Definitions/Direction.hpp"
0015 #include "Acts/Definitions/TrackParametrization.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/Geometry/GeometryIdentifier.hpp"
0020 #include "Acts/MagneticField/ConstantBField.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0023 #include "Acts/MagneticField/NullBField.hpp"
0024 #include "Acts/Propagator/EigenStepper.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/Propagator/StraightLineStepper.hpp"
0027 #include "Acts/Surfaces/PerigeeSurface.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0030 #include "Acts/Utilities/Result.hpp"
0031 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0032 #include "Acts/Vertexing/LinearizedTrack.hpp"
0033 #include "Acts/Vertexing/NumericalTrackLinearizer.hpp"
0034
0035 #include <algorithm>
0036 #include <array>
0037 #include <cmath>
0038 #include <memory>
0039 #include <random>
0040 #include <tuple>
0041 #include <utility>
0042 #include <vector>
0043
0044 namespace bdata = boost::unit_test::data;
0045 using namespace Acts::UnitLiterals;
0046
0047 namespace Acts::Test {
0048
0049 using Covariance = BoundSquareMatrix;
0050 // We will compare analytical and numerical computations in the case of a
0051 // (non-zero) constant B-field and a zero B-field.
0052 using HelicalPropagator = Propagator<EigenStepper<>>;
0053 using StraightPropagator = Propagator<StraightLineStepper>;
0054 using AnalyticalLinearizer = HelicalTrackLinearizer;
0055 using StraightAnalyticalLinearizer = HelicalTrackLinearizer;
0056 using NumericalLinearizer = NumericalTrackLinearizer;
0057 using StraightNumericalLinearizer = NumericalTrackLinearizer;
0058
0059 // Create a test context
0060 GeometryContext geoContext = GeometryContext();
0061 MagneticFieldContext magFieldContext = MagneticFieldContext();
0062
0063 // Vertex x/y position distribution
0064 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0065 // Vertex z position distribution
0066 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0067 // Vertex time distribution
0068 std::uniform_real_distribution<double> vTDist(-1_ns, 1_ns);
0069 // Track d0 distribution
0070 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0071 // Track z0 distribution
0072 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0073 // Track pT distribution
0074 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0075 // Track phi distribution
0076 std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
0077 // Track theta distribution
0078 std::uniform_real_distribution<double> thetaDist(1.0, M_PI - 1.0);
0079 // Track charge helper distribution
0080 std::uniform_real_distribution<double> qDist(-1, 1);
0081 // Track time distribution
0082 std::uniform_real_distribution<double> tDist(-0.002_ns, 0.002_ns);
0083 // Track IP resolution distribution
0084 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0085 // Track angular distribution
0086 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0087 // Track q/p resolution distribution
0088 std::uniform_real_distribution<double> resQoPDist(0.0, 0.1);
0089 // Track time resolution distribution
0090 std::uniform_real_distribution<double> resTDist(0.1_ns, 0.2_ns);
0091
0092 ///
0093 /// @brief Test HelicalTrackLinearizer by comparing it to NumericalTrackLinearizer.
0094 ///
0095 /// @note While HelicalTrackLinearizer computes the Jacobians using analytically derived formulae,
0096 /// NumericalTrackLinearizer uses numerical differentiation:
0097 /// f'(x) ~ (f(x+dx) - f(x)) / dx).
0098 ///
0099 BOOST_AUTO_TEST_CASE(linearized_track_factory_test) {
0100   // Number of tracks to linearize
0101   unsigned int nTracks = 100;
0102
0103   // Set up RNG
0104   int seed = 31415;
0105   std::mt19937 gen(seed);
0106
0107   // Constant B-Field and 0 B-field
0108   auto constField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 2_T});
0109   auto zeroField = std::make_shared<NullBField>();
0110
0111   // Set up stepper and propagator for constant B-field
0112   EigenStepper<> stepper(constField);
0113   auto propagator = std::make_shared<HelicalPropagator>(stepper);
0114
0115   // Set up stepper and propagator for 0 B-field
0116   StraightLineStepper straightStepper;
0117   auto straightPropagator =
0118       std::make_shared<StraightPropagator>(straightStepper);
0119
0120   // Create perigee surface, initial track parameters will be relative to it
0121   std::shared_ptr<PerigeeSurface> perigeeSurface{
0122       Surface::makeShared<PerigeeSurface>(Vector3{0., 0., 0.})};
0123
0124   // Vertex position and corresponding d0 and z0
0125   Vector4 vtxPos;
0126   double d0v{};
0127   double z0v{};
0128   double t0v{};
0129   {
0130     double x = vXYDist(gen);
0131     double y = vXYDist(gen);
0132     double z = vZDist(gen);
0133     double t = vTDist(gen);
0134     d0v = std::hypot(x, y);
0135     z0v = z;
0136     t0v = t;
0137     vtxPos << x, y, z, t;
0138   }
0139
0140   // Vector storing the tracks that we linearize
0141   std::vector<BoundTrackParameters> tracks;
0142
0143   // Construct random track emerging from vicinity of vertex position
0144   for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
0145     // Random charge
0146     double q = qDist(gen) < 0 ? -1. : 1.;
0147
0148     // Random track parameters
0149     BoundVector paramVec;
0150     paramVec << d0v + d0Dist(gen), z0v + z0Dist(gen), phiDist(gen),
0151         thetaDist(gen), q / pTDist(gen), t0v + tDist(gen);
0152
0153     // Resolutions
0154     double resD0 = resIPDist(gen);
0155     double resZ0 = resIPDist(gen);
0156     double resPh = resAngDist(gen);
0157     double resTh = resAngDist(gen);
0158     double resQp = resQoPDist(gen);
0159     double resT = resTDist(gen);
0160
0161     // Fill vector of track objects with simple covariance matrix
0162     Covariance covMat;
0163
0164     covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
0165         0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
0166         0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., resT * resT;
0167     tracks.emplace_back(perigeeSurface, paramVec, std::move(covMat),
0168                         ParticleHypothesis::pion());
0169   }
0170
0171   // Linearizer for constant field and corresponding state
0172   AnalyticalLinearizer::Config linConfig;
0173   linConfig.bField = constField;
0174   linConfig.propagator = propagator;
0175   AnalyticalLinearizer linFactory(linConfig);
0176
0177   NumericalLinearizer::Config numLinConfig(constField, propagator);
0178   NumericalLinearizer numLinFactory(numLinConfig);
0179
0180   // Linearizer for 0 field and corresponding state
0181   StraightAnalyticalLinearizer::Config straightLinConfig;
0182   straightLinConfig.propagator = straightPropagator;
0183   StraightAnalyticalLinearizer straightLinFactory(straightLinConfig);
0184
0185   StraightNumericalLinearizer::Config numStraightLinConfig(straightPropagator);
0186   StraightNumericalLinearizer numStraightLinFactory(numStraightLinConfig);
0187
0188   MagneticFieldProvider::Cache fieldCache =
0189       constField->makeCache(magFieldContext);
0190   MagneticFieldProvider::Cache zeroFieldCache =
0191       zeroField->makeCache(magFieldContext);
0192
0193   // Lambda for comparing outputs of the two linearization methods
0194   // We compare the linearization result at the PCA to "linPoint"
0195   auto checkLinearizers = [&fieldCache, &zeroFieldCache](
0196                               auto& lin1, auto& lin2,
0197                               const BoundTrackParameters& track,
0198                               const Vector4& linPoint,
0199                               const auto& geometryContext,
0200                               const auto& fieldContext) {
0201     // In addition to comparing the output of the linearizers, we check that
0202     // they return non-zero quantities
0203     BoundVector vecBoundZero = BoundVector::Zero();
0204     BoundSquareMatrix matBoundZero = BoundSquareMatrix::Zero();
0205     ActsMatrix<eBoundSize, 4> matBound2SPZero =
0206         ActsMatrix<eBoundSize, 4>::Zero();
0207     ActsMatrix<eBoundSize, 3> matBound2MomZero =
0208         ActsMatrix<eBoundSize, 3>::Zero();
0209
0210     // We check that the entries of the output quantities either
0211     // -) have a relative difference of less than "relTol"
0212     // or
0213     // -) are both smaller than "small"
0214     double relTol = 5e-4;
0215     double small = 5e-4;
0216
0217     std::shared_ptr<PerigeeSurface> perigee =
0218         Surface::makeShared<PerigeeSurface>(VectorHelpers::position(linPoint));
0219
0220     const LinearizedTrack linTrack1 =
0221         lin1.linearizeTrack(track, linPoint[3], *perigee, geometryContext,
0222                             fieldContext, fieldCache)
0223             .value();
0224     const LinearizedTrack linTrack2 =
0225         lin2.linearizeTrack(track, linPoint[3], *perigee, geometryContext,
0226                             fieldContext, zeroFieldCache)
0227             .value();
0228
0229     // There should be no problem here because both linearizers compute
0230     // "parametersAtPCA" the same way
0231     CHECK_CLOSE_OR_SMALL(linTrack1.parametersAtPCA, linTrack2.parametersAtPCA,
0232                          relTol, small);
0233     BOOST_CHECK_NE(linTrack1.parametersAtPCA, vecBoundZero);
0234     BOOST_CHECK_NE(linTrack2.parametersAtPCA, vecBoundZero);
0235
0236     // Compare position Jacobians
0237     CHECK_CLOSE_OR_SMALL((linTrack1.positionJacobian),
0238                          (linTrack2.positionJacobian), relTol, small);
0239     BOOST_CHECK_NE(linTrack1.positionJacobian, matBound2SPZero);
0240     BOOST_CHECK_NE(linTrack2.positionJacobian, matBound2SPZero);
0241
0242     // Compare momentum Jacobians
0243     CHECK_CLOSE_OR_SMALL((linTrack1.momentumJacobian),
0244                          (linTrack2.momentumJacobian), relTol, small);
0245     BOOST_CHECK_NE(linTrack1.momentumJacobian, matBound2MomZero);
0246     BOOST_CHECK_NE(linTrack2.momentumJacobian, matBound2MomZero);
0247
0248     // Again, both methods compute "covarianceAtPCA" the same way => this
0249     // check should always work
0250     CHECK_CLOSE_OR_SMALL(linTrack1.covarianceAtPCA, linTrack2.covarianceAtPCA,
0251                          relTol, small);
0252     BOOST_CHECK_NE(linTrack1.covarianceAtPCA, matBoundZero);
0253     BOOST_CHECK_NE(linTrack2.covarianceAtPCA, matBoundZero);
0254
0255     // Check whether "linPoint" is saved correctly in the LinearizerTrack
0256     // objects
0257     BOOST_CHECK_EQUAL(linTrack1.linearizationPoint, linPoint);
0258     BOOST_CHECK_EQUAL(linTrack2.linearizationPoint, linPoint);
0259
0260     CHECK_CLOSE_OR_SMALL(linTrack1.constantTerm, linTrack2.constantTerm, relTol,
0261                          small);
0262     BOOST_CHECK_NE(linTrack1.constantTerm, vecBoundZero);
0263     BOOST_CHECK_NE(linTrack2.constantTerm, vecBoundZero);
0264   };
0265
0266   // Compare linearizers for all tracks
0267   for (const BoundTrackParameters& trk : tracks) {
0268     BOOST_TEST_CONTEXT("Linearization in constant magnetic field") {
0269       checkLinearizers(linFactory, numLinFactory, trk, vtxPos, geoContext,
0270                        magFieldContext);
0271     }
0272     BOOST_TEST_CONTEXT("Linearization without magnetic field") {
0273       checkLinearizers(straightLinFactory, numStraightLinFactory, trk, vtxPos,
0274                        geoContext, magFieldContext);
0275     }
0276   }
0277 }
0278
0279 }  // namespace Acts::Test