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