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