Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-11-01 07:55:19

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