Core/Vertexing/FullBilloirVertexFitterTests.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/EventData/BoundTrackParameters.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Geometry/GeometryIdentifier.hpp"
0017 #include "Acts/MagneticField/ConstantBField.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Propagator/EigenStepper.hpp"
0020 #include "Acts/Propagator/Propagator.hpp"
0021 #include "Acts/Surfaces/PerigeeSurface.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Utilities/Result.hpp"
0024 #include "Acts/Vertexing/FullBilloirVertexFitter.hpp"
0025 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0026 #include "Acts/Vertexing/TrackAtVertex.hpp"
0027 #include "Acts/Vertexing/Vertex.hpp"
0028 #include "Acts/Vertexing/VertexingOptions.hpp"
0029 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0030
0031 #include <cmath>
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <random>
0036 #include <utility>
0037 #include <vector>
0038
0039 using namespace Acts;
0040 using namespace Acts::UnitLiterals;
0041
0042 namespace ActsTests {
0043
0044 using Covariance = BoundMatrix;
0045
0046 // Create a test context
0047 GeometryContext geoContext = GeometryContext::dangerouslyDefaultConstruct();
0048 MagneticFieldContext magFieldContext = MagneticFieldContext();
0049
0050 // 4D vertex distributions
0051 // x-/y-position
0052 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0053 // z-position
0054 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0055 // time
0056 std::uniform_real_distribution<double> vTDist(-1_ns, 1_ns);
0057
0058 // Track parameter distributions
0059 // d0
0060 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0061 // z0
0062 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0063 // pT
0064 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0065 // phi
0066 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0067                                                std::numbers::pi);
0068 // theta
0069 std::uniform_real_distribution<double> thetaDist(1., std::numbers::pi - 1.);
0070 // charge helper
0071 std::uniform_real_distribution<double> qDist(-1, 1);
0072 // time
0073 std::uniform_real_distribution<double> tDist(-0.002_ns, 0.002_ns);
0074
0075 // Track parameter resolution distributions
0076 // impact parameters
0077 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0078 // angles
0079 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0080 // q/p
0081 std::uniform_real_distribution<double> resQoPDist(-0.1, 0.1);
0082 // Track time resolution distribution
0083 std::uniform_real_distribution<double> resTDist(0.1_ns, 0.2_ns);
0084
0085 // Number of tracks distritbution
0086 std::uniform_int_distribution<std::uint32_t> nTracksDist(3, 10);
0087
0088 // Dummy user-defined InputTrack type
0089 struct InputTrackStub {
0090   explicit InputTrackStub(const BoundTrackParameters& params)
0091       : m_parameters(params) {}
0092
0093   const BoundTrackParameters& parameters() const { return m_parameters; }
0094
0095   // store e.g. link to original objects here
0096
0097  private:
0098   BoundTrackParameters m_parameters;
0099 };
0100
0101 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0102
0103 ///
0104 /// @brief Unit test for FullBilloirVertexFitter
0105 /// with default input track type (= BoundTrackParameters)
0106 ///
0107 BOOST_AUTO_TEST_CASE(billoir_vertex_fitter_defaulttrack_test) {
0108   // Set up RNG
0109   int seed = 31415;
0110   std::mt19937 gen(seed);
0111   // Set up constant B-Field
0112   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0113
0114   // Set up Eigenstepper
0115   EigenStepper<> stepper(bField);
0116   // Set up propagator with void navigator
0117   auto propagator = std::make_shared<Propagator<EigenStepper<>>>(stepper);
0118
0119   HelicalTrackLinearizer::Config ltConfig;
0120   ltConfig.bField = bField;
0121   ltConfig.propagator = propagator;
0122   HelicalTrackLinearizer linearizer(ltConfig);
0123
0124   // Constraint for vertex fit
0125   Vertex constraint;
0126   Vertex customConstraint;
0127   // Some arbitrary values
0128   SquareMatrix4 covMatVtx = SquareMatrix4::Zero();
0129   double ns2 = Acts::UnitConstants::ns * Acts::UnitConstants::ns;
0130   covMatVtx(0, 0) = 30_mm2;
0131   covMatVtx(1, 1) = 30_mm2;
0132   covMatVtx(2, 2) = 30_mm2;
0133   covMatVtx(3, 3) = 30 * ns2;
0134   constraint.setFullCovariance(covMatVtx);
0135   constraint.setFullPosition(Vector4(0, 0, 0, 0));
0136   customConstraint.setFullCovariance(covMatVtx);
0137   customConstraint.setFullPosition(Vector4(0, 0, 0, 0));
0138
0139   // Set up Billoir vertex fitter with default tracks
0140   using VertexFitter = FullBilloirVertexFitter;
0141   VertexFitter::Config vertexFitterCfg;
0142   vertexFitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0143   vertexFitterCfg.trackLinearizer
0144       .connect<&HelicalTrackLinearizer::linearizeTrack>(&linearizer);
0145   VertexFitter billoirFitter(vertexFitterCfg);
0146   auto fieldCache = bField->makeCache(magFieldContext);
0147   // Vertexing options for default tracks
0148   VertexingOptions vfOptions(geoContext, magFieldContext);
0149   VertexingOptions vfOptionsConstr(geoContext, magFieldContext, constraint);
0150
0151   // Create a custom std::function to extract BoundTrackParameters from
0152   // user-defined InputTrack
0153   auto extractParameters = [](const InputTrack& params) {
0154     return params.as<InputTrackStub>()->parameters();
0155   };
0156
0157   // Set up Billoir vertex fitter with user-defined input tracks
0158   VertexFitter::Config customVertexFitterCfg;
0159   customVertexFitterCfg.extractParameters.connect(extractParameters);
0160   customVertexFitterCfg.trackLinearizer
0161       .connect<&HelicalTrackLinearizer::linearizeTrack>(&linearizer);
0162   VertexFitter customBilloirFitter(customVertexFitterCfg);
0163   // Vertexing options for custom tracks
0164   VertexingOptions customVfOptions(geoContext, magFieldContext);
0165   VertexingOptions customVfOptionsConstr(geoContext, magFieldContext,
0166                                          customConstraint);
0167
0168   BOOST_TEST_CONTEXT(
0169       "Testing FullBilloirVertexFitter when input track vector is empty.") {
0170     const std::vector<const BoundTrackParameters*> emptyVector;
0171     const std::vector<InputTrack> emptyVectorInput;
0172
0173     // Without constraint
0174     Vertex fittedVertex =
0175         billoirFitter.fit(emptyVectorInput, vfOptions, fieldCache).value();
0176
0177     Vector3 origin(0., 0., 0.);
0178     SquareMatrix4 zeroMat = SquareMatrix4::Zero();
0179     BOOST_CHECK_EQUAL(fittedVertex.position(), origin);
0180     BOOST_CHECK_EQUAL(fittedVertex.fullCovariance(), zeroMat);
0181
0182     // With constraint
0183     fittedVertex =
0184         billoirFitter.fit(emptyVectorInput, vfOptionsConstr, fieldCache)
0185             .value();
0186
0187     BOOST_CHECK_EQUAL(fittedVertex.position(), origin);
0188     BOOST_CHECK_EQUAL(fittedVertex.fullCovariance(), zeroMat);
0189   }
0190
0191   // Number of events
0192   const int nEvents = 100;
0193   for (int eventIdx = 0; eventIdx < nEvents; ++eventIdx) {
0194     // Number of tracks
0195     std::uint32_t nTracks = nTracksDist(gen);
0196
0197     // Create position of vertex and perigee surface
0198     double x = vXYDist(gen);
0199     double y = vXYDist(gen);
0200     double z = vZDist(gen);
0201     double t = vTDist(gen);
0202
0203     Vector4 trueVertex(x, y, z, t);
0204     std::shared_ptr<PerigeeSurface> perigeeSurface =
0205         Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0206
0207     // Calculate d0 and z0 corresponding to the vertex position
0208     double d0V = std::hypot(x, y);
0209     double z0V = z;
0210
0211     // Vector to store track objects used for vertex fit
0212     std::vector<BoundTrackParameters> tracks;
0213     std::vector<InputTrackStub> customTracks;
0214
0215     // Calculate random track emerging from vicinity of vertex position
0216     for (std::uint32_t iTrack = 0; iTrack < nTracks; iTrack++) {
0217       // Charge
0218       double q = std::copysign(1., qDist(gen));
0219
0220       // Track parameters
0221       BoundVector paramVec;
0222       paramVec << d0V + d0Dist(gen), z0V + z0Dist(gen), phiDist(gen),
0223           thetaDist(gen), q / pTDist(gen), t + tDist(gen);
0224
0225       // Resolutions
0226       double resD0 = resIPDist(gen);
0227       double resZ0 = resIPDist(gen);
0228       double resPh = resAngDist(gen);
0229       double resTh = resAngDist(gen);
0230       double resQp = resQoPDist(gen);
0231       double resT = resTDist(gen);
0232
0233       // Fill vector of track objects with simple covariance matrix
0234       Covariance covMat;
0235
0236       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0237           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0238           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
0239           resT * resT;
0240       tracks.emplace_back(BoundTrackParameters(perigeeSurface, paramVec, covMat,
0241                                                ParticleHypothesis::pion()));
0242       customTracks.emplace_back(
0243           BoundTrackParameters(perigeeSurface, paramVec, std::move(covMat),
0244                                ParticleHypothesis::pion()));
0245     }
0246
0247     std::vector<InputTrack> inputTracks;
0248     for (const auto& trk : tracks) {
0249       inputTracks.push_back(InputTrack{&trk});
0250     }
0251
0252     std::vector<InputTrack> customInputTracks;
0253     for (const auto& trk : customTracks) {
0254       customInputTracks.push_back(InputTrack{&trk});
0255     }
0256
0257     auto fit = [&trueVertex, &nTracks, &fieldCache](const auto& fitter,
0258                                                     const auto& trksPtr,
0259                                                     const auto& vfOpts) {
0260       auto fittedVertex = fitter.fit(trksPtr, vfOpts, fieldCache).value();
0261       if (!fittedVertex.tracks().empty()) {
0262         CHECK_CLOSE_ABS(fittedVertex.position(), trueVertex.head(3), 1_mm);
0263         auto tracksAtVtx = fittedVertex.tracks();
0264         auto trackAtVtx = tracksAtVtx[0];
0265         CHECK_CLOSE_ABS(fittedVertex.time(), trueVertex[3], 1_ns);
0266       }
0267
0268       std::cout << "\nFitting " << nTracks << " tracks" << std::endl;
0269       std::cout << "True Vertex:\n" << trueVertex << std::endl;
0270       std::cout << "Fitted Vertex:\n"
0271                 << fittedVertex.fullPosition() << std::endl;
0272     };
0273
0274     BOOST_TEST_CONTEXT(
0275         "Testing FullBilloirVertexFitter without vertex constraint.") {
0276       fit(billoirFitter, inputTracks, vfOptions);
0277     }
0278     BOOST_TEST_CONTEXT(
0279         "Testing FullBilloirVertexFitter with vertex constraint.") {
0280       fit(billoirFitter, inputTracks, vfOptionsConstr);
0281     }
0282     BOOST_TEST_CONTEXT(
0283         "Testing FullBilloirVertexFitter with custom tracks (no vertex "
0284         "constraint).") {
0285       fit(customBilloirFitter, customInputTracks, customVfOptions);
0286     }
0287     BOOST_TEST_CONTEXT(
0288         "Testing FullBilloirVertexFitter with custom tracks (with vertex "
0289         "constraint).") {
0290       fit(customBilloirFitter, customInputTracks, customVfOptionsConstr);
0291     }
0292   }
0293 }
0294
0295 BOOST_AUTO_TEST_SUITE_END()
0296
0297 }  // namespace ActsTests