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