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/Direction.hpp"
0016 #include "Acts/Definitions/TrackParametrization.hpp"
0017 #include "Acts/Definitions/Units.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/Geometry/GeometryContext.hpp"
0020 #include "Acts/Geometry/GeometryIdentifier.hpp"
0021 #include "Acts/MagneticField/ConstantBField.hpp"
0022 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0023 #include "Acts/Propagator/EigenStepper.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Surfaces/PerigeeSurface.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0028 #include "Acts/Utilities/Result.hpp"
0029 #include "Acts/Vertexing/FullBilloirVertexFitter.hpp"
0030 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0031 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0032 #include "Acts/Vertexing/VertexingOptions.hpp"
0033 #include "Acts/Vertexing/ZScanVertexFinder.hpp"
0034 
0035 #include <algorithm>
0036 #include <array>
0037 #include <cmath>
0038 #include <functional>
0039 #include <iostream>
0040 #include <memory>
0041 #include <random>
0042 #include <string>
0043 #include <system_error>
0044 #include <tuple>
0045 #include <type_traits>
0046 #include <utility>
0047 #include <vector>
0048 
0049 namespace bdata = boost::unit_test::data;
0050 using namespace Acts::UnitLiterals;
0051 
0052 namespace Acts::Test {
0053 
0054 using Covariance = BoundSquareMatrix;
0055 using Propagator = Acts::Propagator<EigenStepper<>>;
0056 using Linearizer_t = HelicalTrackLinearizer;
0057 
0058 // Create a test context
0059 GeometryContext geoContext = GeometryContext();
0060 MagneticFieldContext magFieldContext = MagneticFieldContext();
0061 
0062 // Vertex x/y position distribution
0063 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0064 // Vertex z position distribution
0065 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0066 // Track d0 distribution
0067 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0068 // Track z0 distribution
0069 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0070 // Track pT distribution
0071 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0072 // Track phi distribution
0073 std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
0074 // Track theta distribution
0075 std::uniform_real_distribution<double> thetaDist(1.0, M_PI - 1.0);
0076 // Track charge helper distribution
0077 std::uniform_real_distribution<double> qDist(-1, 1);
0078 // Track IP resolution distribution
0079 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0080 // Track angular distribution
0081 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0082 // Track q/p resolution distribution
0083 std::uniform_real_distribution<double> resQoPDist(-0.01, 0.01);
0084 
0085 ///
0086 /// @brief Unit test for ZScanVertexFinder
0087 ///
0088 BOOST_AUTO_TEST_CASE(zscan_finder_test) {
0089   unsigned int nTests = 50;
0090 
0091   for (unsigned int iTest = 0; iTest < nTests; ++iTest) {
0092     // Number of tracks
0093     unsigned int nTracks = 30;
0094 
0095     // Set up RNG
0096     int mySeed = 31415;
0097     std::mt19937 gen(mySeed);
0098 
0099     // Set up constant B-Field
0100     auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0101 
0102     // Set up Eigenstepper
0103     EigenStepper<> stepper(bField);
0104 
0105     // Set up propagator with void navigator
0106     auto propagator = std::make_shared<Propagator>(stepper);
0107 
0108     using BilloirFitter = FullBilloirVertexFitter;
0109 
0110     // Create perigee surface
0111     std::shared_ptr<PerigeeSurface> perigeeSurface =
0112         Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0113 
0114     // Create position of vertex and perigee surface
0115     double x = vXYDist(gen);
0116     double y = vXYDist(gen);
0117     double z = vZDist(gen);
0118 
0119     // Calculate d0 and z0 corresponding to vertex position
0120     double d0_v = std::hypot(x, y);
0121     double z0_v = z;
0122 
0123     // Start constructing nTracks tracks in the following
0124     std::vector<BoundTrackParameters> tracks;
0125 
0126     // Construct random track emerging from vicinity of vertex position
0127     // Vector to store track objects used for vertex fit
0128     for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
0129       // Construct positive or negative charge randomly
0130       double q = qDist(gen) < 0 ? -1. : 1.;
0131 
0132       // Construct random track parameters
0133       BoundVector paramVec = BoundVector::Zero();
0134       paramVec[eBoundLoc0] = d0_v + d0Dist(gen);
0135       paramVec[eBoundLoc1] = z0_v + z0Dist(gen);
0136       paramVec[eBoundPhi] = phiDist(gen);
0137       paramVec[eBoundTheta] = thetaDist(gen);
0138       paramVec[eBoundQOverP] = q / pTDist(gen);
0139 
0140       // Resolutions
0141       double resD0 = resIPDist(gen);
0142       double resZ0 = resIPDist(gen);
0143       double resPh = resAngDist(gen);
0144       double resTh = resAngDist(gen);
0145       double resQp = resQoPDist(gen);
0146 
0147       // Fill vector of track objects with simple covariance matrix
0148       Covariance covMat;
0149 
0150       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0151           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0152           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
0153 
0154       tracks.emplace_back(perigeeSurface, paramVec, std::move(covMat),
0155                           ParticleHypothesis::pion());
0156     }
0157 
0158     std::vector<InputTrack> inputTracks;
0159     for (const auto& trk : tracks) {
0160       inputTracks.emplace_back(&trk);
0161     }
0162 
0163     using VertexFinder = ZScanVertexFinder;
0164 
0165     ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0166     ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0167 
0168     VertexFinder::Config cfg(ipEstimator);
0169     cfg.extractParameters.connect<&InputTrack::extractParameters>();
0170 
0171     VertexFinder finder(cfg);
0172 
0173     VertexingOptions vertexingOptions(geoContext, magFieldContext);
0174 
0175     auto state = finder.makeState(magFieldContext);
0176     auto res = finder.find(inputTracks, vertexingOptions, state);
0177 
0178     BOOST_CHECK(res.ok());
0179 
0180     if (!res.ok()) {
0181       std::cout << res.error().message() << std::endl;
0182     }
0183 
0184     if (res.ok()) {
0185       BOOST_CHECK(!(*res).empty());
0186       Vector3 result = (*res).back().position();
0187       CHECK_CLOSE_ABS(result[eZ], z, 1_mm);
0188     }
0189   }
0190 }
0191 
0192 // Dummy user-defined InputTrackStub type
0193 struct InputTrackStub {
0194   InputTrackStub(const BoundTrackParameters& params) : m_parameters(params) {}
0195 
0196   const BoundTrackParameters& parameters() const { return m_parameters; }
0197 
0198   // store e.g. link to original objects here
0199 
0200  private:
0201   BoundTrackParameters m_parameters;
0202 };
0203 
0204 ///
0205 /// @brief Unit test for ZScanVertexFinder with user-defined input track type
0206 ///
0207 BOOST_AUTO_TEST_CASE(zscan_finder_usertrack_test) {
0208   unsigned int nTests = 50;
0209 
0210   for (unsigned int iTest = 0; iTest < nTests; ++iTest) {
0211     // Number of tracks
0212     unsigned int nTracks = 30;
0213 
0214     // Set up RNG
0215     int mySeed = 31415;
0216     std::mt19937 gen(mySeed);
0217 
0218     // Set up constant B-Field
0219     auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0220 
0221     // Set up Eigenstepper
0222     EigenStepper<> stepper(bField);
0223 
0224     // Set up propagator with void navigator
0225     auto propagator = std::make_shared<Propagator>(stepper);
0226 
0227     using BilloirFitter = FullBilloirVertexFitter;
0228 
0229     // Create perigee surface
0230     std::shared_ptr<PerigeeSurface> perigeeSurface =
0231         Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0232 
0233     // Create position of vertex and perigee surface
0234     double x = vXYDist(gen);
0235     double y = vXYDist(gen);
0236     double z = vZDist(gen);
0237 
0238     // Calculate d0 and z0 corresponding to vertex position
0239     double d0_v = std::hypot(x, y);
0240     double z0_v = z;
0241 
0242     // Start constructing nTracks tracks in the following
0243     std::vector<InputTrackStub> tracks;
0244 
0245     // Construct random track emerging from vicinity of vertex position
0246     // Vector to store track objects used for vertex fit
0247     for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
0248       // Construct positive or negative charge randomly
0249       double q = qDist(gen) < 0 ? -1. : 1.;
0250 
0251       // Construct random track parameters
0252       BoundVector paramVec;
0253       double z0track = z0_v + z0Dist(gen);
0254       paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
0255           q / pTDist(gen), 0.;
0256 
0257       // Resolutions
0258       double resD0 = resIPDist(gen);
0259       double resZ0 = resIPDist(gen);
0260       double resPh = resAngDist(gen);
0261       double resTh = resAngDist(gen);
0262       double resQp = resQoPDist(gen);
0263 
0264       // Fill vector of track objects with simple covariance matrix
0265       Covariance covMat;
0266 
0267       covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0.,
0268           0., 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh,
0269           0., 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
0270 
0271       tracks.emplace_back(BoundTrackParameters(perigeeSurface, paramVec,
0272                                                std::move(covMat),
0273                                                ParticleHypothesis::pion()));
0274     }
0275 
0276     std::vector<InputTrack> inputTracks;
0277     for (const auto& trk : tracks) {
0278       inputTracks.emplace_back(&trk);
0279     }
0280 
0281     using VertexFinder = ZScanVertexFinder;
0282 
0283     ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0284     ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0285 
0286     VertexFinder::Config cfg(ipEstimator);
0287 
0288     // Create a custom std::function to extract BoundTrackParameters from
0289     // user-defined InputTrackStub
0290     auto extractParameters = [](const InputTrack& params) {
0291       return params.as<InputTrackStub>()->parameters();
0292     };
0293 
0294     cfg.extractParameters.connect(extractParameters);
0295     VertexFinder finder(cfg);
0296     auto state = finder.makeState(magFieldContext);
0297 
0298     VertexingOptions vertexingOptions(geoContext, magFieldContext);
0299 
0300     auto res = finder.find(inputTracks, vertexingOptions, state);
0301 
0302     BOOST_CHECK(res.ok());
0303 
0304     if (!res.ok()) {
0305       std::cout << res.error().message() << std::endl;
0306     }
0307 
0308     if (res.ok()) {
0309       BOOST_CHECK(!(*res).empty());
0310       Vector3 result = (*res).back().position();
0311       CHECK_CLOSE_ABS(result[eZ], z, 1_mm);
0312     }
0313   }
0314 }
0315 
0316 }  // namespace Acts::Test