Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-16 07:37:22

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/Common.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/BoundTrackParameters.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0018 #include "Acts/Surfaces/PerigeeSurface.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/Intersection.hpp"
0021 #include "Acts/Utilities/Result.hpp"
0022 #include "Acts/Utilities/UnitVectors.hpp"
0023 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
0024 #include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
0025 #include "Acts/Vertexing/Vertex.hpp"
0026 #include "Acts/Vertexing/VertexingOptions.hpp"
0027 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0028 
0029 #include <iostream>
0030 #include <memory>
0031 #include <numbers>
0032 #include <random>
0033 #include <system_error>
0034 #include <vector>
0035 
0036 using namespace Acts;
0037 using namespace Acts::UnitLiterals;
0038 using Acts::VectorHelpers::makeVector4;
0039 
0040 namespace ActsTests {
0041 
0042 using Covariance = BoundMatrix;
0043 
0044 // Create a test context
0045 GeometryContext geoContext = GeometryContext::dangerouslyDefaultConstruct();
0046 MagneticFieldContext magFieldContext = MagneticFieldContext();
0047 
0048 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0049 ///
0050 /// @brief Unit test for TrackDensityVertexFinder using same configuration
0051 /// and values as VertexSeedFinderTestAlg in Athena implementation, i.e.
0052 /// tests if reordering tracks returns the same result
0053 ///
0054 BOOST_AUTO_TEST_CASE(track_density_finder_test) {
0055   // Define some track parameter properties
0056   Vector3 pos0{0, 0, 0};
0057   Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
0058   Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
0059   Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
0060   Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
0061   Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
0062   Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
0063 
0064   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0065   GaussianTrackDensity::Config densityCfg;
0066   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0067   TrackDensityVertexFinder finder{
0068       TrackDensityVertexFinder::Config{Acts::GaussianTrackDensity(densityCfg)}};
0069   auto state = finder.makeState(magFieldContext);
0070 
0071   // Start creating some track parameters
0072   Covariance covMat = Covariance::Identity();
0073   std::shared_ptr<PerigeeSurface> perigeeSurface =
0074       Surface::makeShared<PerigeeSurface>(pos0);
0075 
0076   // Test finder for some fixed track parameter values
0077   auto params1a =
0078       BoundTrackParameters::create(
0079           geoContext, perigeeSurface, makeVector4(pos1a, 0), mom1a.normalized(),
0080           1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
0081           .value();
0082   auto params1b =
0083       BoundTrackParameters::create(
0084           geoContext, perigeeSurface, makeVector4(pos1b, 0), mom1b.normalized(),
0085           -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
0086           .value();
0087   auto params1c =
0088       BoundTrackParameters::create(
0089           geoContext, perigeeSurface, makeVector4(pos1c, 0), mom1c.normalized(),
0090           1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
0091           .value();
0092 
0093   // Vectors of track parameters in different orders
0094   std::vector<InputTrack> vec1 = {InputTrack{&params1a}, InputTrack{&params1b},
0095                                   InputTrack{&params1c}};
0096   std::vector<InputTrack> vec2 = {InputTrack{&params1c}, InputTrack{&params1a},
0097                                   InputTrack{&params1b}};
0098 
0099   auto res1 = finder.find(vec1, vertexingOptions, state);
0100   auto res2 = finder.find(vec2, vertexingOptions, state);
0101 
0102   if (!res1.ok()) {
0103     std::cout << res1.error().message() << std::endl;
0104   }
0105 
0106   if (!res2.ok()) {
0107     std::cout << res2.error().message() << std::endl;
0108   }
0109 
0110   if (res1.ok() && res2.ok()) {
0111     BOOST_CHECK(!(*res1).empty());
0112     BOOST_CHECK(!(*res2).empty());
0113     Vector3 result1 = (*res1).back().position();
0114     Vector3 result2 = (*res2).back().position();
0115     BOOST_CHECK_EQUAL(result1, result2);
0116   }
0117 }
0118 
0119 ///
0120 /// @brief Unit test for TrackDensityVertexFinder using same configuration
0121 /// and values as VertexSeedFinderTestAlg in Athena implementation
0122 ///
0123 BOOST_AUTO_TEST_CASE(track_density_finder_constr_test) {
0124   // Define some track parameter properties
0125   Vector3 pos0{0, 0, 0};
0126   Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
0127   Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
0128   Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
0129   Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
0130   Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
0131   Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
0132 
0133   // From Athena VertexSeedFinderTestAlg
0134   double const expectedZResult = -13.013;
0135 
0136   // Create constraint for seed finding
0137   Vector3 constraintPos{1.7_mm, 1.3_mm, -6_mm};
0138   SquareMatrix3 constrCov = SquareMatrix<3>::Identity();
0139 
0140   Vertex constraint(constraintPos);
0141   constraint.setCovariance(constrCov);
0142 
0143   // Finder options
0144   VertexingOptions vertexingOptions(geoContext, magFieldContext, constraint);
0145   GaussianTrackDensity::Config densityCfg;
0146   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0147   TrackDensityVertexFinder finder{
0148       TrackDensityVertexFinder::Config{Acts::GaussianTrackDensity(densityCfg)}};
0149   auto state = finder.makeState(magFieldContext);
0150 
0151   // Start creating some track parameters
0152   Covariance covMat = Covariance::Identity();
0153   std::shared_ptr<PerigeeSurface> perigeeSurface =
0154       Surface::makeShared<PerigeeSurface>(pos0);
0155 
0156   // Test finder for some fixed track parameter values
0157   auto params1a =
0158       BoundTrackParameters::create(
0159           geoContext, perigeeSurface, makeVector4(pos1a, 0), mom1a.normalized(),
0160           1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
0161           .value();
0162   auto params1b =
0163       BoundTrackParameters::create(
0164           geoContext, perigeeSurface, makeVector4(pos1b, 0), mom1b.normalized(),
0165           -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
0166           .value();
0167   auto params1c =
0168       BoundTrackParameters::create(
0169           geoContext, perigeeSurface, makeVector4(pos1c, 0), mom1c.normalized(),
0170           -1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
0171           .value();
0172 
0173   // Vector of track parameters
0174   std::vector<InputTrack> vec1 = {InputTrack{&params1a}, InputTrack{&params1b},
0175                                   InputTrack{&params1c}};
0176 
0177   auto res = finder.find(vec1, vertexingOptions, state);
0178 
0179   if (!res.ok()) {
0180     std::cout << res.error().message() << std::endl;
0181   }
0182 
0183   if (res.ok()) {
0184     BOOST_CHECK(!(*res).empty());
0185     Vector3 result = (*res).back().position();
0186 
0187     BOOST_CHECK_EQUAL(result[eX], constraintPos[eX]);
0188     BOOST_CHECK_EQUAL(result[eY], constraintPos[eY]);
0189     CHECK_CLOSE_ABS(result[eZ], expectedZResult, 0.001_mm);
0190   }
0191 }
0192 
0193 const double zVertexPos = 12.;
0194 // x position
0195 std::normal_distribution<double> xdist(1_mm, 0.1_mm);
0196 // y position
0197 std::normal_distribution<double> ydist(-0.7_mm, 0.1_mm);
0198 // z1 position
0199 std::normal_distribution<double> z1dist(zVertexPos * 1_mm, 1_mm);
0200 // z2 position
0201 std::normal_distribution<double> z2dist(-3_mm, 0.5_mm);
0202 // Track pT distribution
0203 std::uniform_real_distribution<double> pTDist(0.1_GeV, 100_GeV);
0204 // Track phi distribution
0205 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0206                                                std::numbers::pi);
0207 // Track eta distribution
0208 std::uniform_real_distribution<double> etaDist(-4., 4.);
0209 
0210 ///
0211 /// @brief Unit test for TrackDensityVertexFinder using same configuration
0212 /// and values as VertexSeedFinderTestAlg in Athena implementation
0213 ///
0214 BOOST_AUTO_TEST_CASE(track_density_finder_random_test) {
0215   Covariance covMat = Covariance::Identity();
0216 
0217   // Perigee surface for track parameters
0218   Vector3 pos0{0, 0, 0};
0219   std::shared_ptr<PerigeeSurface> perigeeSurface =
0220       Surface::makeShared<PerigeeSurface>(pos0);
0221 
0222   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0223   GaussianTrackDensity::Config densityCfg;
0224   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0225   TrackDensityVertexFinder finder{
0226       TrackDensityVertexFinder::Config{Acts::GaussianTrackDensity(densityCfg)}};
0227   auto state = finder.makeState(magFieldContext);
0228 
0229   int mySeed = 31415;
0230   std::mt19937 gen(mySeed);
0231   unsigned int nTracks = 200;
0232 
0233   std::vector<BoundTrackParameters> trackVec;
0234   trackVec.reserve(nTracks);
0235 
0236   // Create nTracks tracks for test case
0237   for (unsigned int i = 0; i < nTracks; i++) {
0238     // The position of the particle
0239     Vector3 pos(xdist(gen), ydist(gen), 0);
0240 
0241     // Create momentum and charge of track
0242     double pt = pTDist(gen);
0243     double phi = phiDist(gen);
0244     double eta = etaDist(gen);
0245     double charge = std::copysign(1., etaDist(gen));
0246 
0247     // project the position on the surface
0248     Vector3 direction = makeDirectionFromPhiEta(phi, eta);
0249     Intersection3D intersection =
0250         perigeeSurface->intersect(geoContext, pos, direction).closest();
0251     pos = intersection.position();
0252 
0253     // Produce most of the tracks at near z1 position,
0254     // some near z2. Highest track density then expected at z1
0255     pos[eZ] = ((i % 4) == 0) ? z2dist(gen) : z1dist(gen);
0256 
0257     trackVec.push_back(BoundTrackParameters::create(
0258                            geoContext, perigeeSurface, makeVector4(pos, 0),
0259                            direction, charge / pt, covMat,
0260                            ParticleHypothesis::pion())
0261                            .value());
0262   }
0263 
0264   std::vector<InputTrack> inputTracks;
0265   for (const auto& trk : trackVec) {
0266     inputTracks.emplace_back(&trk);
0267   }
0268 
0269   auto res3 = finder.find(inputTracks, vertexingOptions, state);
0270   if (!res3.ok()) {
0271     std::cout << res3.error().message() << std::endl;
0272   }
0273 
0274   if (res3.ok()) {
0275     BOOST_CHECK(!(*res3).empty());
0276     Vector3 result = (*res3).back().position();
0277     CHECK_CLOSE_ABS(result[eZ], zVertexPos, 1_mm);
0278   }
0279 }
0280 
0281 // Dummy user-defined InputTrackStub type
0282 struct InputTrackStub {
0283   explicit InputTrackStub(const BoundTrackParameters& params)
0284       : m_parameters(params) {}
0285 
0286   const BoundTrackParameters& parameters() const { return m_parameters; }
0287 
0288   // store e.g. link to original objects here
0289 
0290  private:
0291   BoundTrackParameters m_parameters;
0292 };
0293 
0294 ///
0295 /// @brief Unit test for TrackDensityVertexFinder with user-defined input track
0296 /// type with same values as in other tests
0297 ///
0298 BOOST_AUTO_TEST_CASE(track_density_finder_usertrack_test) {
0299   // Define some track parameter properties
0300   Vector3 pos0{0, 0, 0};
0301   Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
0302   Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
0303   Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
0304   Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
0305   Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
0306   Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
0307 
0308   // From Athena VertexSeedFinderTestAlg
0309   double const expectedZResult = -13.013;
0310 
0311   // Create constraint for seed finding
0312   Vector3 constraintPos{1.7_mm, 1.3_mm, -6_mm};
0313   SquareMatrix3 constrCov = SquareMatrix3::Identity();
0314 
0315   Vertex constraint(constraintPos);
0316   constraint.setCovariance(constrCov);
0317 
0318   // Finder options
0319   VertexingOptions vertexingOptions(geoContext, magFieldContext, constraint);
0320 
0321   auto extractParameters = [](const InputTrack& params) {
0322     return params.as<InputTrackStub>()->parameters();
0323   };
0324 
0325   GaussianTrackDensity::Config densityCfg;
0326   densityCfg.extractParameters.connect(extractParameters);
0327   TrackDensityVertexFinder finder{
0328       TrackDensityVertexFinder::Config{Acts::GaussianTrackDensity(densityCfg)}};
0329   auto state = finder.makeState(magFieldContext);
0330 
0331   // Start creating some track parameters
0332   Covariance covMat = Covariance::Identity();
0333   std::shared_ptr<PerigeeSurface> perigeeSurface =
0334       Surface::makeShared<PerigeeSurface>(pos0);
0335 
0336   // Test finder for some fixed track parameter values
0337   InputTrackStub params1a(BoundTrackParameters::create(
0338                               geoContext, perigeeSurface, makeVector4(pos1a, 0),
0339                               mom1a, 1_e / mom1a.norm(), covMat,
0340                               ParticleHypothesis::pion())
0341                               .value());
0342   InputTrackStub params1b(BoundTrackParameters::create(
0343                               geoContext, perigeeSurface, makeVector4(pos1b, 0),
0344                               mom1b, -1_e / mom1b.norm(), covMat,
0345                               ParticleHypothesis::pion())
0346                               .value());
0347   InputTrackStub params1c(BoundTrackParameters::create(
0348                               geoContext, perigeeSurface, makeVector4(pos1c, 0),
0349                               mom1c, -1_e / mom1c.norm(), covMat,
0350                               ParticleHypothesis::pion())
0351                               .value());
0352 
0353   // Vector of track parameters
0354   std::vector<InputTrack> vec1 = {InputTrack{&params1a}, InputTrack{&params1b},
0355                                   InputTrack{&params1c}};
0356 
0357   auto res = finder.find(vec1, vertexingOptions, state);
0358 
0359   if (!res.ok()) {
0360     std::cout << res.error().message() << std::endl;
0361   }
0362 
0363   if (res.ok()) {
0364     BOOST_CHECK(!(*res).empty());
0365     Vector3 result = (*res).back().position();
0366 
0367     BOOST_CHECK_EQUAL(result[eX], constraintPos[eX]);
0368     BOOST_CHECK_EQUAL(result[eY], constraintPos[eY]);
0369     CHECK_CLOSE_ABS(result[eZ], expectedZResult, 0.001_mm);
0370   }
0371 }
0372 
0373 BOOST_AUTO_TEST_SUITE_END()
0374 
0375 }  // namespace ActsTests