Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-23 07:36:30

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/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/BoundTrackParameters.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Geometry/GeometryIdentifier.hpp"
0017 #include "Acts/MagneticField/ConstantBField.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Propagator/EigenStepper.hpp"
0020 #include "Acts/Propagator/Propagator.hpp"
0021 #include "Acts/Surfaces/PerigeeSurface.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Utilities/Result.hpp"
0024 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0025 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0026 #include "Acts/Vertexing/KalmanVertexUpdater.hpp"
0027 #include "Acts/Vertexing/LinearizedTrack.hpp"
0028 #include "Acts/Vertexing/TrackAtVertex.hpp"
0029 #include "Acts/Vertexing/Vertex.hpp"
0030 
0031 #include <cmath>
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <optional>
0036 #include <random>
0037 #include <utility>
0038 #include <vector>
0039 
0040 using namespace Acts;
0041 using namespace Acts::UnitLiterals;
0042 
0043 namespace ActsTests {
0044 
0045 using Covariance = BoundMatrix;
0046 using Propagator = Acts::Propagator<EigenStepper<>>;
0047 using Linearizer = HelicalTrackLinearizer;
0048 
0049 // Create a test context
0050 GeometryContext geoContext = GeometryContext::dangerouslyDefaultConstruct();
0051 MagneticFieldContext magFieldContext = MagneticFieldContext();
0052 
0053 // Vertex x/y position distribution
0054 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0055 // Vertex z position distribution
0056 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0057 // Track d0 distribution
0058 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0059 // Track z0 distribution
0060 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0061 // Track pT distribution
0062 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0063 // Track phi distribution
0064 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0065                                                std::numbers::pi);
0066 // Track theta distribution
0067 std::uniform_real_distribution<double> thetaDist(1., std::numbers::pi - 1.);
0068 // Track charge helper distribution
0069 std::uniform_real_distribution<double> qDist(-1, 1);
0070 // Track IP resolution distribution
0071 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0072 // Track angular distribution
0073 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0074 // Track q/p resolution distribution
0075 std::uniform_real_distribution<double> resQoPDist(-0.01, 0.01);
0076 // Number of vertices per test event distribution
0077 
0078 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0079 ///
0080 /// @brief Unit test for KalmanVertexUpdater
0081 ///
0082 BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) {
0083   bool debug = false;
0084 
0085   // Number of tests
0086   unsigned int nTests = 10;
0087 
0088   // Set up RNG
0089   int mySeed = 31415;
0090   std::mt19937 gen(mySeed);
0091 
0092   // Set up constant B-Field
0093   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0094 
0095   // Set up Eigenstepper
0096   EigenStepper<> stepper(bField);
0097 
0098   // Set up propagator with void navigator
0099   auto propagator = std::make_shared<Propagator>(stepper);
0100 
0101   // Linearizer for BoundTrackParameters type test
0102   Linearizer::Config ltConfig;
0103   ltConfig.bField = bField;
0104   ltConfig.propagator = propagator;
0105   Linearizer linearizer(ltConfig);
0106   auto fieldCache = bField->makeCache(magFieldContext);
0107 
0108   // Create perigee surface at origin
0109   std::shared_ptr<PerigeeSurface> perigeeSurface =
0110       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0111 
0112   // Creates a random tracks around origin and a random vertex.
0113   // VertexUpdater adds track to vertex and updates the position
0114   // which should afterwards be closer to the origin/track
0115   for (unsigned int i = 0; i < nTests; ++i) {
0116     if (debug) {
0117       std::cout << "Test " << i + 1 << std::endl;
0118     }
0119     // Construct positive or negative charge randomly
0120     double q = std::copysign(1., qDist(gen));
0121 
0122     // Construct random track parameters around origin
0123     BoundVector paramVec;
0124 
0125     paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen),
0126         q / pTDist(gen), 0.;
0127 
0128     if (debug) {
0129       std::cout << "Creating track parameters: " << paramVec << std::endl;
0130     }
0131 
0132     // Fill vector of track objects with simple covariance matrix
0133     Covariance covMat;
0134 
0135     // Resolutions
0136     double res_d0 = resIPDist(gen);
0137     double res_z0 = resIPDist(gen);
0138     double res_ph = resAngDist(gen);
0139     double res_th = resAngDist(gen);
0140     double res_qp = resQoPDist(gen);
0141 
0142     covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0., 0.,
0143         0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
0144         res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0., 0.,
0145         0., 0., 0., 1.;
0146     BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat),
0147                                 ParticleHypothesis::pion());
0148 
0149     std::shared_ptr<PerigeeSurface> perigee =
0150         Surface::makeShared<PerigeeSurface>(Vector3::Zero());
0151 
0152     // Linearized state of the track
0153     LinearizedTrack linTrack =
0154         linearizer
0155             .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
0156                             fieldCache)
0157             .value();
0158 
0159     // Create TrackAtVertex
0160     TrackAtVertex trkAtVtx(0., params, InputTrack{&params});
0161 
0162     // Set linearized state of trackAtVertex
0163     trkAtVtx.linearizedState = linTrack;
0164     trkAtVtx.isLinearized = true;
0165 
0166     // Create a vertex
0167     Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen));
0168     Vertex vtx(vtxPos);
0169     vtx.setFullCovariance(SquareMatrix4::Identity() * 0.01);
0170 
0171     // Update trkAtVertex with assumption of originating from vtx
0172     KalmanVertexUpdater::updateVertexWithTrack(vtx, trkAtVtx, 3);
0173 
0174     if (debug) {
0175       std::cout << "Old vertex position: " << vtxPos << std::endl;
0176       std::cout << "New vertex position: " << vtx.position() << std::endl;
0177     }
0178 
0179     double oldDistance = vtxPos.norm();
0180     double newDistance = vtx.position().norm();
0181 
0182     if (debug) {
0183       std::cout << "Old distance: " << oldDistance << std::endl;
0184       std::cout << "New distance: " << newDistance << std::endl;
0185     }
0186 
0187     // After update, vertex should be closer to the track
0188     BOOST_CHECK_LT(newDistance, oldDistance);
0189 
0190     // Note: KalmanVertexUpdater updates the vertex w.r.t. the
0191     // newly given track, but does NOT add the track to the
0192     // TrackAtVertex list. Has to be done manually after calling
0193     // the update method.
0194     BOOST_CHECK(vtx.tracks().empty());
0195 
0196   }  // end for loop
0197 
0198 }  // end test case
0199 
0200 ///
0201 /// @brief Unit test for KalmanVertexTrackUpdater
0202 ///
0203 BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) {
0204   bool debug = true;
0205 
0206   // Number of tests
0207   unsigned int nTests = 10;
0208 
0209   // Set up RNG
0210   int mySeed = 31415;
0211   std::mt19937 gen(mySeed);
0212 
0213   // Set up constant B-Field
0214   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0215 
0216   // Set up Eigenstepper
0217   EigenStepper<> stepper(bField);
0218 
0219   // Set up propagator with void navigator
0220   auto propagator = std::make_shared<Propagator>(stepper);
0221 
0222   // Set up ImpactPointEstimator, used for comparisons later
0223   ImpactPointEstimator::Config ip3dEstConfig(bField, propagator);
0224   ImpactPointEstimator ip3dEst(ip3dEstConfig);
0225   ImpactPointEstimator::State state{bField->makeCache(magFieldContext)};
0226 
0227   // Set up HelicalTrackLinearizer, needed for linearizing the tracks
0228   // Linearizer for BoundTrackParameters type test
0229   Linearizer::Config ltConfig;
0230   ltConfig.bField = bField;
0231   ltConfig.propagator = propagator;
0232   Linearizer linearizer(ltConfig);
0233   auto fieldCache = bField->makeCache(magFieldContext);
0234 
0235   // Create perigee surface at origin
0236   std::shared_ptr<PerigeeSurface> perigeeSurface =
0237       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0238 
0239   // Create random tracks around origin and a random vertex.
0240   // Update tracks with the assumption that they originate from
0241   // the vertex position and check if they are closer to the
0242   // vertex after the update process
0243   for (unsigned int i = 0; i < nTests; ++i) {
0244     // Construct positive or negative charge randomly
0245     double q = std::copysign(1., qDist(gen));
0246 
0247     // Construct random track parameters
0248     BoundVector paramVec;
0249 
0250     paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen),
0251         q / pTDist(gen), 0.;
0252 
0253     if (debug) {
0254       std::cout << "Creating track parameters: " << paramVec << std::endl;
0255     }
0256 
0257     // Fill vector of track objects with simple covariance matrix
0258     Covariance covMat;
0259 
0260     // Resolutions
0261     double res_d0 = resIPDist(gen);
0262     double res_z0 = resIPDist(gen);
0263     double res_ph = resAngDist(gen);
0264     double res_th = resAngDist(gen);
0265     double res_qp = resQoPDist(gen);
0266 
0267     covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0., 0.,
0268         0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
0269         res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0., 0.,
0270         0., 0., 0., 1.;
0271     BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat),
0272                                 ParticleHypothesis::pion());
0273 
0274     std::shared_ptr<PerigeeSurface> perigee =
0275         Surface::makeShared<PerigeeSurface>(Vector3::Zero());
0276 
0277     // Linearized state of the track
0278     LinearizedTrack linTrack =
0279         linearizer
0280             .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
0281                             fieldCache)
0282             .value();
0283 
0284     // Create TrackAtVertex
0285     TrackAtVertex trkAtVtx(0., params, InputTrack{&params});
0286 
0287     // Set linearized state of trackAtVertex
0288     trkAtVtx.linearizedState = linTrack;
0289     trkAtVtx.isLinearized = true;
0290 
0291     // Copy parameters for later comparison of old and new version
0292     auto fittedParamsCopy = trkAtVtx.fittedParams;
0293 
0294     // Create a vertex
0295     Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen));
0296     Vertex vtx(vtxPos);
0297 
0298     // Update trkAtVertex with assumption of originating from vtx
0299     KalmanVertexUpdater::updateTrackWithVertex(trkAtVtx, vtx, 3);
0300 
0301     // The old distance
0302     double oldDistance =
0303         ip3dEst.calculateDistance(geoContext, fittedParamsCopy, vtxPos, state)
0304             .value();
0305 
0306     // The new distance after update
0307     double newDistance =
0308         ip3dEst
0309             .calculateDistance(geoContext, trkAtVtx.fittedParams, vtxPos, state)
0310             .value();
0311     if (debug) {
0312       std::cout << "Old distance: " << oldDistance << std::endl;
0313       std::cout << "New distance: " << newDistance << std::endl;
0314     }
0315 
0316     // Parameters should have changed
0317     BOOST_CHECK_NE(fittedParamsCopy, trkAtVtx.fittedParams);
0318 
0319     // After update, track should be closer to the vertex
0320     BOOST_CHECK_LT(newDistance, oldDistance);
0321 
0322   }  // end for loop
0323 
0324 }  // end test case
0325 
0326 BOOST_AUTO_TEST_SUITE_END()
0327 
0328 }  // namespace ActsTests