Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:13:02

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{{{densityCfg}}};
0073   auto state = finder.makeState(magFieldContext);
0074 
0075   // Start creating some track parameters
0076   Covariance covMat = Covariance::Identity();
0077   std::shared_ptr<PerigeeSurface> perigeeSurface =
0078       Surface::makeShared<PerigeeSurface>(pos0);
0079 
0080   // Test finder for some fixed track parameter values
0081   auto params1a =
0082       BoundTrackParameters::create(
0083           perigeeSurface, geoContext, makeVector4(pos1a, 0), mom1a.normalized(),
0084           1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
0085           .value();
0086   auto params1b =
0087       BoundTrackParameters::create(
0088           perigeeSurface, geoContext, makeVector4(pos1b, 0), mom1b.normalized(),
0089           -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
0090           .value();
0091   auto params1c =
0092       BoundTrackParameters::create(
0093           perigeeSurface, geoContext, makeVector4(pos1c, 0), mom1c.normalized(),
0094           1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
0095           .value();
0096 
0097   // Vectors of track parameters in different orders
0098   std::vector<InputTrack> vec1 = {InputTrack{&params1a}, InputTrack{&params1b},
0099                                   InputTrack{&params1c}};
0100   std::vector<InputTrack> vec2 = {InputTrack{&params1c}, InputTrack{&params1a},
0101                                   InputTrack{&params1b}};
0102 
0103   auto res1 = finder.find(vec1, vertexingOptions, state);
0104   auto res2 = finder.find(vec2, vertexingOptions, state);
0105 
0106   if (!res1.ok()) {
0107     std::cout << res1.error().message() << std::endl;
0108   }
0109 
0110   if (!res2.ok()) {
0111     std::cout << res2.error().message() << std::endl;
0112   }
0113 
0114   if (res1.ok() && res2.ok()) {
0115     BOOST_CHECK(!(*res1).empty());
0116     BOOST_CHECK(!(*res2).empty());
0117     Vector3 result1 = (*res1).back().position();
0118     Vector3 result2 = (*res2).back().position();
0119     BOOST_CHECK_EQUAL(result1, result2);
0120   }
0121 }
0122 
0123 ///
0124 /// @brief Unit test for TrackDensityVertexFinder using same configuration
0125 /// and values as VertexSeedFinderTestAlg in Athena implementation
0126 ///
0127 BOOST_AUTO_TEST_CASE(track_density_finder_constr_test) {
0128   // Define some track parameter properties
0129   Vector3 pos0{0, 0, 0};
0130   Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
0131   Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
0132   Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
0133   Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
0134   Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
0135   Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
0136 
0137   // From Athena VertexSeedFinderTestAlg
0138   double const expectedZResult = -13.013;
0139 
0140   // Create constraint for seed finding
0141   Vector3 constraintPos{1.7_mm, 1.3_mm, -6_mm};
0142   SquareMatrix3 constrCov = ActsSquareMatrix<3>::Identity();
0143 
0144   Vertex constraint(constraintPos);
0145   constraint.setCovariance(constrCov);
0146 
0147   // Finder options
0148   VertexingOptions vertexingOptions(geoContext, magFieldContext, constraint);
0149   GaussianTrackDensity::Config densityCfg;
0150   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0151   TrackDensityVertexFinder finder{{{densityCfg}}};
0152   auto state = finder.makeState(magFieldContext);
0153 
0154   // Start creating some track parameters
0155   Covariance covMat = Covariance::Identity();
0156   std::shared_ptr<PerigeeSurface> perigeeSurface =
0157       Surface::makeShared<PerigeeSurface>(pos0);
0158 
0159   // Test finder for some fixed track parameter values
0160   auto params1a =
0161       BoundTrackParameters::create(
0162           perigeeSurface, geoContext, makeVector4(pos1a, 0), mom1a.normalized(),
0163           1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
0164           .value();
0165   auto params1b =
0166       BoundTrackParameters::create(
0167           perigeeSurface, geoContext, makeVector4(pos1b, 0), mom1b.normalized(),
0168           -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
0169           .value();
0170   auto params1c =
0171       BoundTrackParameters::create(
0172           perigeeSurface, geoContext, makeVector4(pos1c, 0), mom1c.normalized(),
0173           -1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
0174           .value();
0175 
0176   // Vector of track parameters
0177   std::vector<InputTrack> vec1 = {InputTrack{&params1a}, InputTrack{&params1b},
0178                                   InputTrack{&params1c}};
0179 
0180   auto res = finder.find(vec1, vertexingOptions, state);
0181 
0182   if (!res.ok()) {
0183     std::cout << res.error().message() << std::endl;
0184   }
0185 
0186   if (res.ok()) {
0187     BOOST_CHECK(!(*res).empty());
0188     Vector3 result = (*res).back().position();
0189 
0190     BOOST_CHECK_EQUAL(result[eX], constraintPos[eX]);
0191     BOOST_CHECK_EQUAL(result[eY], constraintPos[eY]);
0192     CHECK_CLOSE_ABS(result[eZ], expectedZResult, 0.001_mm);
0193   }
0194 }
0195 
0196 const double zVertexPos = 12.;
0197 // x position
0198 std::normal_distribution<double> xdist(1_mm, 0.1_mm);
0199 // y position
0200 std::normal_distribution<double> ydist(-0.7_mm, 0.1_mm);
0201 // z1 position
0202 std::normal_distribution<double> z1dist(zVertexPos * 1_mm, 1_mm);
0203 // z2 position
0204 std::normal_distribution<double> z2dist(-3_mm, 0.5_mm);
0205 // Track pT distribution
0206 std::uniform_real_distribution<double> pTDist(0.1_GeV, 100_GeV);
0207 // Track phi distribution
0208 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0209                                                std::numbers::pi);
0210 // Track eta distribution
0211 std::uniform_real_distribution<double> etaDist(-4., 4.);
0212 
0213 ///
0214 /// @brief Unit test for TrackDensityVertexFinder using same configuration
0215 /// and values as VertexSeedFinderTestAlg in Athena implementation
0216 ///
0217 BOOST_AUTO_TEST_CASE(track_density_finder_random_test) {
0218   Covariance covMat = Covariance::Identity();
0219 
0220   // Perigee surface for track parameters
0221   Vector3 pos0{0, 0, 0};
0222   std::shared_ptr<PerigeeSurface> perigeeSurface =
0223       Surface::makeShared<PerigeeSurface>(pos0);
0224 
0225   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0226   GaussianTrackDensity::Config densityCfg;
0227   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0228   TrackDensityVertexFinder finder{{{densityCfg}}};
0229   auto state = finder.makeState(magFieldContext);
0230 
0231   int mySeed = 31415;
0232   std::mt19937 gen(mySeed);
0233   unsigned int nTracks = 200;
0234 
0235   std::vector<BoundTrackParameters> trackVec;
0236   trackVec.reserve(nTracks);
0237 
0238   // Create nTracks tracks for test case
0239   for (unsigned int i = 0; i < nTracks; i++) {
0240     // The position of the particle
0241     Vector3 pos(xdist(gen), ydist(gen), 0);
0242 
0243     // Create momentum and charge of track
0244     double pt = pTDist(gen);
0245     double phi = phiDist(gen);
0246     double eta = etaDist(gen);
0247     double charge = etaDist(gen) > 0 ? 1 : -1;
0248 
0249     // project the position on the surface
0250     Vector3 direction = makeDirectionFromPhiEta(phi, eta);
0251     auto intersection =
0252         perigeeSurface->intersect(geoContext, pos, direction).closest();
0253     pos = intersection.position();
0254 
0255     // Produce most of the tracks at near z1 position,
0256     // some near z2. Highest track density then expected at z1
0257     pos[eZ] = ((i % 4) == 0) ? z2dist(gen) : z1dist(gen);
0258 
0259     trackVec.push_back(BoundTrackParameters::create(
0260                            perigeeSurface, geoContext, makeVector4(pos, 0),
0261                            direction, charge / pt, covMat,
0262                            ParticleHypothesis::pion())
0263                            .value());
0264   }
0265 
0266   std::vector<InputTrack> inputTracks;
0267   for (const auto& trk : trackVec) {
0268     inputTracks.emplace_back(&trk);
0269   }
0270 
0271   auto res3 = finder.find(inputTracks, vertexingOptions, state);
0272   if (!res3.ok()) {
0273     std::cout << res3.error().message() << std::endl;
0274   }
0275 
0276   if (res3.ok()) {
0277     BOOST_CHECK(!(*res3).empty());
0278     Vector3 result = (*res3).back().position();
0279     CHECK_CLOSE_ABS(result[eZ], zVertexPos, 1_mm);
0280   }
0281 }
0282 
0283 // Dummy user-defined InputTrackStub type
0284 struct InputTrackStub {
0285   InputTrackStub(const BoundTrackParameters& params) : 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{{{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                               perigeeSurface, geoContext, makeVector4(pos1a, 0),
0339                               mom1a, 1_e / mom1a.norm(), covMat,
0340                               ParticleHypothesis::pion())
0341                               .value());
0342   InputTrackStub params1b(BoundTrackParameters::create(
0343                               perigeeSurface, geoContext, makeVector4(pos1b, 0),
0344                               mom1b, -1_e / mom1b.norm(), covMat,
0345                               ParticleHypothesis::pion())
0346                               .value());
0347   InputTrackStub params1c(BoundTrackParameters::create(
0348                               perigeeSurface, geoContext, 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 }  // namespace Acts::Test