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