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