Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-02 07:52:31

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