Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:11:37

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