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   explicit InputTrackStub(const BoundTrackParameters& params)
0098       : m_parameters(params) {}
0099
0100   const BoundTrackParameters& parameters() const { return m_parameters; }
0101
0102   // store e.g. link to original objects here
0103
0104  private:
0105   BoundTrackParameters m_parameters;
0106 };
0107
0108 ///
0109 /// @brief Unit test for FullBilloirVertexFitter
0110 /// with default input track type (= BoundTrackParameters)
0111 ///
0112 BOOST_AUTO_TEST_CASE(billoir_vertex_fitter_defaulttrack_test) {
0113   // Set up RNG
0114   int seed = 31415;
0115   std::mt19937 gen(seed);
0116   // Set up constant B-Field
0117   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0118
0119   // Set up Eigenstepper
0120   EigenStepper<> stepper(bField);
0121   // Set up propagator with void navigator
0122   auto propagator = std::make_shared<Propagator<EigenStepper<>>>(stepper);
0123
0124   HelicalTrackLinearizer::Config ltConfig;
0125   ltConfig.bField = bField;
0126   ltConfig.propagator = propagator;
0127   HelicalTrackLinearizer linearizer(ltConfig);
0128
0129   // Constraint for vertex fit
0130   Vertex constraint;
0131   Vertex customConstraint;
0132   // Some arbitrary values
0133   SquareMatrix4 covMatVtx = SquareMatrix4::Zero();
0134   double ns2 = Acts::UnitConstants::ns * Acts::UnitConstants::ns;
0135   covMatVtx(0, 0) = 30_mm2;
0136   covMatVtx(1, 1) = 30_mm2;
0137   covMatVtx(2, 2) = 30_mm2;
0138   covMatVtx(3, 3) = 30 * ns2;
0139   constraint.setFullCovariance(covMatVtx);
0140   constraint.setFullPosition(Vector4(0, 0, 0, 0));
0141   customConstraint.setFullCovariance(covMatVtx);
0142   customConstraint.setFullPosition(Vector4(0, 0, 0, 0));
0143
0144   // Set up Billoir vertex fitter with default tracks
0145   using VertexFitter = FullBilloirVertexFitter;
0146   VertexFitter::Config vertexFitterCfg;
0147   vertexFitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0148   vertexFitterCfg.trackLinearizer
0149       .connect<&HelicalTrackLinearizer::linearizeTrack>(&linearizer);
0150   VertexFitter billoirFitter(vertexFitterCfg);
0151   auto fieldCache = bField->makeCache(magFieldContext);
0152   // Vertexing options for default tracks
0153   VertexingOptions vfOptions(geoContext, magFieldContext);
0154   VertexingOptions vfOptionsConstr(geoContext, magFieldContext, constraint);
0155
0156   // Create a custom std::function to extract BoundTrackParameters from
0157   // user-defined InputTrack
0158   auto extractParameters = [](const InputTrack& params) {
0159     return params.as<InputTrackStub>()->parameters();
0160   };
0161
0162   // Set up Billoir vertex fitter with user-defined input tracks
0163   VertexFitter::Config customVertexFitterCfg;
0164   customVertexFitterCfg.extractParameters.connect(extractParameters);
0165   customVertexFitterCfg.trackLinearizer
0166       .connect<&HelicalTrackLinearizer::linearizeTrack>(&linearizer);
0167   VertexFitter customBilloirFitter(customVertexFitterCfg);
0168   // Vertexing options for custom tracks
0169   VertexingOptions customVfOptions(geoContext, magFieldContext);
0170   VertexingOptions customVfOptionsConstr(geoContext, magFieldContext,
0171                                          customConstraint);
0172
0173   BOOST_TEST_CONTEXT(
0174       "Testing FullBilloirVertexFitter when input track vector is empty.") {
0175     const std::vector<const BoundTrackParameters*> emptyVector;
0176     const std::vector<InputTrack> emptyVectorInput;
0177
0178     // Without constraint
0179     Vertex fittedVertex =
0180         billoirFitter.fit(emptyVectorInput, vfOptions, fieldCache).value();
0181
0182     Vector3 origin(0., 0., 0.);
0183     SquareMatrix4 zeroMat = SquareMatrix4::Zero();
0184     BOOST_CHECK_EQUAL(fittedVertex.position(), origin);
0185     BOOST_CHECK_EQUAL(fittedVertex.fullCovariance(), zeroMat);
0186
0187     // With constraint
0188     fittedVertex =
0189         billoirFitter.fit(emptyVectorInput, vfOptionsConstr, fieldCache)
0190             .value();
0191
0192     BOOST_CHECK_EQUAL(fittedVertex.position(), origin);
0193     BOOST_CHECK_EQUAL(fittedVertex.fullCovariance(), zeroMat);
0194   }
0195
0196   // Number of events
0197   const int nEvents = 100;
0198   for (int eventIdx = 0; eventIdx < nEvents; ++eventIdx) {
0199     // Number of tracks
0200     std::uint32_t nTracks = nTracksDist(gen);
0201
0202     // Create position of vertex and perigee surface
0203     double x = vXYDist(gen);
0204     double y = vXYDist(gen);
0205     double z = vZDist(gen);
0206     double t = vTDist(gen);
0207
0208     Vector4 trueVertex(x, y, z, t);
0209     std::shared_ptr<PerigeeSurface> perigeeSurface =
0210         Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0211
0212     // Calculate d0 and z0 corresponding to the vertex position
0213     double d0V = std::hypot(x, y);
0214     double z0V = z;
0215
0216     // Vector to store track objects used for vertex fit
0217     std::vector<BoundTrackParameters> tracks;
0218     std::vector<InputTrackStub> customTracks;
0219
0220     // Calculate random track emerging from vicinity of vertex position
0221     for (std::uint32_t iTrack = 0; iTrack < nTracks; iTrack++) {
0222       // Charge
0223       double q = qDist(gen) < 0 ? -1. : 1.;
0224
0225       // Track parameters
0226       BoundVector paramVec;
0227       paramVec << d0V + d0Dist(gen), z0V + z0Dist(gen), phiDist(gen),
0228           thetaDist(gen), q / pTDist(gen), t + tDist(gen);
0229
0230       // Resolutions
0231       double resD0 = resIPDist(gen);
0232       double resZ0 = resIPDist(gen);
0233       double resPh = resAngDist(gen);
0234       double resTh = resAngDist(gen);
0235       double resQp = resQoPDist(gen);
0236       double resT = resTDist(gen);
0237
0238       // Fill vector of track objects with simple covariance matrix
0239       Covariance covMat;
0240
0241       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0242           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0243           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0.,
0244           resT * resT;
0245       tracks.emplace_back(BoundTrackParameters(perigeeSurface, paramVec, covMat,
0246                                                ParticleHypothesis::pion()));
0247       customTracks.emplace_back(
0248           BoundTrackParameters(perigeeSurface, paramVec, std::move(covMat),
0249                                ParticleHypothesis::pion()));
0250     }
0251
0252     std::vector<InputTrack> inputTracks;
0253     for (const auto& trk : tracks) {
0254       inputTracks.push_back(InputTrack{&trk});
0255     }
0256
0257     std::vector<InputTrack> customInputTracks;
0258     for (const auto& trk : customTracks) {
0259       customInputTracks.push_back(InputTrack{&trk});
0260     }
0261
0262     auto fit = [&trueVertex, &nTracks, &fieldCache](const auto& fitter,
0263                                                     const auto& trksPtr,
0264                                                     const auto& vfOpts) {
0265       auto fittedVertex = fitter.fit(trksPtr, vfOpts, fieldCache).value();
0266       if (!fittedVertex.tracks().empty()) {
0267         CHECK_CLOSE_ABS(fittedVertex.position(), trueVertex.head(3), 1_mm);
0268         auto tracksAtVtx = fittedVertex.tracks();
0269         auto trackAtVtx = tracksAtVtx[0];
0270         CHECK_CLOSE_ABS(fittedVertex.time(), trueVertex[3], 1_ns);
0271       }
0272
0273       std::cout << "\nFitting " << nTracks << " tracks" << std::endl;
0274       std::cout << "True Vertex:\n" << trueVertex << std::endl;
0275       std::cout << "Fitted Vertex:\n"
0276                 << fittedVertex.fullPosition() << std::endl;
0277     };
0278
0279     BOOST_TEST_CONTEXT(
0280         "Testing FullBilloirVertexFitter without vertex constraint.") {
0281       fit(billoirFitter, inputTracks, vfOptions);
0282     }
0283     BOOST_TEST_CONTEXT(
0284         "Testing FullBilloirVertexFitter with vertex constraint.") {
0285       fit(billoirFitter, inputTracks, vfOptionsConstr);
0286     }
0287     BOOST_TEST_CONTEXT(
0288         "Testing FullBilloirVertexFitter with custom tracks (no vertex "
0289         "constraint).") {
0290       fit(customBilloirFitter, customInputTracks, customVfOptions);
0291     }
0292     BOOST_TEST_CONTEXT(
0293         "Testing FullBilloirVertexFitter with custom tracks (with vertex "
0294         "constraint).") {
0295       fit(customBilloirFitter, customInputTracks, customVfOptionsConstr);
0296     }
0297   }
0298 }
0299 }  // namespace Acts::Test