File indexing completed on 2026-05-23 07:36:30
0001
0002
0003
0004
0005
0006
0007
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
0050 GeometryContext geoContext = GeometryContext::dangerouslyDefaultConstruct();
0051 MagneticFieldContext magFieldContext = MagneticFieldContext();
0052
0053
0054 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0055
0056 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0057
0058 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0059
0060 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0061
0062 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0063
0064 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0065 std::numbers::pi);
0066
0067 std::uniform_real_distribution<double> thetaDist(1., std::numbers::pi - 1.);
0068
0069 std::uniform_real_distribution<double> qDist(-1, 1);
0070
0071 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0072
0073 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0074
0075 std::uniform_real_distribution<double> resQoPDist(-0.01, 0.01);
0076
0077
0078 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0079
0080
0081
0082 BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) {
0083 bool debug = false;
0084
0085
0086 unsigned int nTests = 10;
0087
0088
0089 int mySeed = 31415;
0090 std::mt19937 gen(mySeed);
0091
0092
0093 auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0094
0095
0096 EigenStepper<> stepper(bField);
0097
0098
0099 auto propagator = std::make_shared<Propagator>(stepper);
0100
0101
0102 Linearizer::Config ltConfig;
0103 ltConfig.bField = bField;
0104 ltConfig.propagator = propagator;
0105 Linearizer linearizer(ltConfig);
0106 auto fieldCache = bField->makeCache(magFieldContext);
0107
0108
0109 std::shared_ptr<PerigeeSurface> perigeeSurface =
0110 Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0111
0112
0113
0114
0115 for (unsigned int i = 0; i < nTests; ++i) {
0116 if (debug) {
0117 std::cout << "Test " << i + 1 << std::endl;
0118 }
0119
0120 double q = std::copysign(1., qDist(gen));
0121
0122
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
0133 Covariance covMat;
0134
0135
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
0153 LinearizedTrack linTrack =
0154 linearizer
0155 .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
0156 fieldCache)
0157 .value();
0158
0159
0160 TrackAtVertex trkAtVtx(0., params, InputTrack{¶ms});
0161
0162
0163 trkAtVtx.linearizedState = linTrack;
0164 trkAtVtx.isLinearized = true;
0165
0166
0167 Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen));
0168 Vertex vtx(vtxPos);
0169 vtx.setFullCovariance(SquareMatrix4::Identity() * 0.01);
0170
0171
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
0188 BOOST_CHECK_LT(newDistance, oldDistance);
0189
0190
0191
0192
0193
0194 BOOST_CHECK(vtx.tracks().empty());
0195
0196 }
0197
0198 }
0199
0200
0201
0202
0203 BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) {
0204 bool debug = true;
0205
0206
0207 unsigned int nTests = 10;
0208
0209
0210 int mySeed = 31415;
0211 std::mt19937 gen(mySeed);
0212
0213
0214 auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0215
0216
0217 EigenStepper<> stepper(bField);
0218
0219
0220 auto propagator = std::make_shared<Propagator>(stepper);
0221
0222
0223 ImpactPointEstimator::Config ip3dEstConfig(bField, propagator);
0224 ImpactPointEstimator ip3dEst(ip3dEstConfig);
0225 ImpactPointEstimator::State state{bField->makeCache(magFieldContext)};
0226
0227
0228
0229 Linearizer::Config ltConfig;
0230 ltConfig.bField = bField;
0231 ltConfig.propagator = propagator;
0232 Linearizer linearizer(ltConfig);
0233 auto fieldCache = bField->makeCache(magFieldContext);
0234
0235
0236 std::shared_ptr<PerigeeSurface> perigeeSurface =
0237 Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0238
0239
0240
0241
0242
0243 for (unsigned int i = 0; i < nTests; ++i) {
0244
0245 double q = std::copysign(1., qDist(gen));
0246
0247
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
0258 Covariance covMat;
0259
0260
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
0278 LinearizedTrack linTrack =
0279 linearizer
0280 .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
0281 fieldCache)
0282 .value();
0283
0284
0285 TrackAtVertex trkAtVtx(0., params, InputTrack{¶ms});
0286
0287
0288 trkAtVtx.linearizedState = linTrack;
0289 trkAtVtx.isLinearized = true;
0290
0291
0292 auto fittedParamsCopy = trkAtVtx.fittedParams;
0293
0294
0295 Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen));
0296 Vertex vtx(vtxPos);
0297
0298
0299 KalmanVertexUpdater::updateTrackWithVertex(trkAtVtx, vtx, 3);
0300
0301
0302 double oldDistance =
0303 ip3dEst.calculateDistance(geoContext, fittedParamsCopy, vtxPos, state)
0304 .value();
0305
0306
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
0317 BOOST_CHECK_NE(fittedParamsCopy, trkAtVtx.fittedParams);
0318
0319
0320 BOOST_CHECK_LT(newDistance, oldDistance);
0321
0322 }
0323
0324 }
0325
0326 BOOST_AUTO_TEST_SUITE_END()
0327
0328 }