File indexing completed on 2026-01-07 09:25:04
0001
0002
0003
0004
0005 #include "IterativeVertexFinder.h"
0006
0007 #include <Acts/Definitions/Units.hpp>
0008 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0009 #include <Acts/Propagator/EigenStepper.hpp>
0010 #include <Acts/Propagator/Propagator.hpp>
0011 #include <Acts/Propagator/VoidNavigator.hpp>
0012 #include <Acts/Utilities/Logger.hpp>
0013 #include <Acts/Utilities/Result.hpp>
0014 #include <Acts/Utilities/detail/ContextType.hpp>
0015 #include <Acts/Vertexing/FullBilloirVertexFitter.hpp>
0016 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0017 #include <Acts/Vertexing/IVertexFinder.hpp>
0018 #include <Acts/Vertexing/ImpactPointEstimator.hpp>
0019 #include <Acts/Vertexing/IterativeVertexFinder.hpp>
0020 #include <Acts/Vertexing/LinearizedTrack.hpp>
0021 #include <Acts/Vertexing/TrackAtVertex.hpp>
0022 #include <Acts/Vertexing/Vertex.hpp>
0023 #include <Acts/Vertexing/VertexingOptions.hpp>
0024 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
0025 #include <ActsExamples/EventData/Track.hpp>
0026 #include <ActsExamples/EventData/Trajectories.hpp>
0027 #include <edm4eic/Cov4f.h>
0028 #include <edm4eic/ReconstructedParticleCollection.h>
0029 #include <edm4eic/Track.h>
0030 #include <edm4eic/TrackParameters.h>
0031 #include <edm4eic/Trajectory.h>
0032 #include <edm4eic/unit_system.h>
0033 #include <edm4hep/Vector2f.h>
0034 #include <edm4hep/Vector4f.h>
0035 #include <fmt/core.h>
0036 #include <fmt/format.h>
0037 #include <podio/RelationRange.h>
0038 #include <Eigen/Core>
0039 #include <cmath>
0040 #include <string>
0041 #include <utility>
0042
0043 #include "extensions/spdlog/SpdlogToActs.h"
0044
0045 void eicrecon::IterativeVertexFinder::init(std::shared_ptr<const ActsGeometryProvider> geo_svc,
0046 std::shared_ptr<spdlog::logger> log) {
0047
0048 m_log = log;
0049
0050 m_geoSvc = geo_svc;
0051
0052 m_BField =
0053 std::dynamic_pointer_cast<const eicrecon::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
0054 m_fieldctx = eicrecon::BField::BFieldVariant(m_BField);
0055 }
0056
0057 std::unique_ptr<edm4eic::VertexCollection> eicrecon::IterativeVertexFinder::produce(
0058 std::vector<const ActsExamples::Trajectories*> trajectories,
0059 const edm4eic::ReconstructedParticleCollection* reconParticles) {
0060
0061 auto outputVertices = std::make_unique<edm4eic::VertexCollection>();
0062
0063 using Propagator = Acts::Propagator<Acts::EigenStepper<>>;
0064 using Linearizer = Acts::HelicalTrackLinearizer;
0065 using VertexFitter = Acts::FullBilloirVertexFitter;
0066 using ImpactPointEstimator = Acts::ImpactPointEstimator;
0067 using VertexSeeder = Acts::ZScanVertexFinder;
0068 using VertexFinder = Acts::IterativeVertexFinder;
0069 using VertexFinderOptions = Acts::VertexingOptions;
0070
0071 ACTS_LOCAL_LOGGER(eicrecon::getSpdlogLogger("IVF", m_log));
0072
0073 Acts::EigenStepper<> stepper(m_BField);
0074
0075
0076 auto propagator = std::make_shared<Propagator>(stepper, Acts::VoidNavigator{},
0077 logger().cloneWithSuffix("Prop"));
0078
0079
0080 Linearizer::Config linearizerCfg;
0081 linearizerCfg.bField = m_BField;
0082 linearizerCfg.propagator = propagator;
0083 Linearizer linearizer(linearizerCfg, logger().cloneWithSuffix("HelLin"));
0084
0085
0086 VertexFitter::Config vertexFitterCfg;
0087 vertexFitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0088 vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0089 VertexFitter vertexFitter(vertexFitterCfg);
0090
0091
0092 ImpactPointEstimator::Config ipEstCfg(m_BField, propagator);
0093 ImpactPointEstimator ipEst(ipEstCfg);
0094 VertexSeeder::Config seederCfg(ipEst);
0095 seederCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0096 auto seeder = std::make_shared<VertexSeeder>(seederCfg);
0097
0098
0099 VertexFinder::Config finderCfg(std::move(vertexFitter), std::move(seeder), std::move(ipEst));
0100 finderCfg.maxVertices = m_cfg.maxVertices;
0101 finderCfg.reassignTracksAfterFirstFit = m_cfg.reassignTracksAfterFirstFit;
0102 finderCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0103 finderCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0104 finderCfg.field = m_BField;
0105 VertexFinder finder(std::move(finderCfg));
0106 Acts::IVertexFinder::State state(std::in_place_type<VertexFinder::State>, *m_BField, m_fieldctx);
0107 VertexFinderOptions finderOpts(m_geoctx, m_fieldctx);
0108
0109 std::vector<Acts::InputTrack> inputTracks;
0110
0111 for (const auto& trajectory : trajectories) {
0112 auto tips = trajectory->tips();
0113 if (tips.empty()) {
0114 continue;
0115 }
0116
0117 for (auto& tip : tips) {
0118 ActsExamples::TrackParameters par = trajectory->trackParameters(tip);
0119
0120 inputTracks.emplace_back(&(trajectory->trackParameters(tip)));
0121 m_log->trace("Track local position at input = {} mm, {} mm",
0122 par.localPosition().x() / Acts::UnitConstants::mm,
0123 par.localPosition().y() / Acts::UnitConstants::mm);
0124 }
0125 }
0126
0127 std::vector<Acts::Vertex> vertices;
0128 auto result = finder.find(inputTracks, finderOpts, state);
0129 if (result.ok()) {
0130 vertices = std::move(result.value());
0131 }
0132
0133 for (const auto& vtx : vertices) {
0134 edm4eic::Cov4f cov(vtx.fullCovariance()(0, 0), vtx.fullCovariance()(1, 1),
0135 vtx.fullCovariance()(2, 2), vtx.fullCovariance()(3, 3),
0136 vtx.fullCovariance()(0, 1), vtx.fullCovariance()(0, 2),
0137 vtx.fullCovariance()(0, 3), vtx.fullCovariance()(1, 2),
0138 vtx.fullCovariance()(1, 3), vtx.fullCovariance()(2, 3));
0139 auto eicvertex = outputVertices->create();
0140 eicvertex.setType(1);
0141 eicvertex.setChi2((float)vtx.fitQuality().first);
0142 eicvertex.setNdf((float)vtx.fitQuality().second);
0143 eicvertex.setPosition({
0144 (float)vtx.position().x(),
0145 (float)vtx.position().y(),
0146 (float)vtx.position().z(),
0147 (float)vtx.time(),
0148 });
0149 eicvertex.setPositionError(cov);
0150
0151 for (const auto& t : vtx.tracks()) {
0152 const auto& par = Acts::InputTrack::extractParameters(t.originalParams);
0153 m_log->trace("Track local position from vertex = {} mm, {} mm",
0154 par.localPosition().x() / Acts::UnitConstants::mm,
0155 par.localPosition().y() / Acts::UnitConstants::mm);
0156 float loc_a = par.localPosition().x();
0157 float loc_b = par.localPosition().y();
0158
0159 for (const auto& part : *reconParticles) {
0160 const auto& tracks = part.getTracks();
0161 for (const auto& trk : tracks) {
0162 const auto& traj = trk.getTrajectory();
0163 const auto& trkPars = traj.getTrackParameters();
0164 for (const auto& par : trkPars) {
0165 const double EPSILON = 1.0e-4;
0166 if (std::abs((par.getLoc().a / edm4eic::unit::mm) - (loc_a / Acts::UnitConstants::mm)) <
0167 EPSILON &&
0168 std::abs((par.getLoc().b / edm4eic::unit::mm) - (loc_b / Acts::UnitConstants::mm)) <
0169 EPSILON) {
0170 m_log->trace(
0171 "From ReconParticles, track local position [Loc a, Loc b] = {} mm, {} mm",
0172 par.getLoc().a / edm4eic::unit::mm, par.getLoc().b / edm4eic::unit::mm);
0173 eicvertex.addToAssociatedParticles(part);
0174 }
0175 }
0176 }
0177 }
0178 }
0179 m_log->debug("One vertex found at (x,y,z) = ({}, {}, {}) mm.",
0180 vtx.position().x() / Acts::UnitConstants::mm,
0181 vtx.position().y() / Acts::UnitConstants::mm,
0182 vtx.position().z() / Acts::UnitConstants::mm);
0183
0184 }
0185
0186 return outputVertices;
0187 }