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