algorithms/tracking/IterativeVertexFinder.cc

0001 // Created by Joe Osborn
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004
0005 #include "IterativeVertexFinder.h"
0006
0007 #include <Acts/Definitions/TrackParametrization.hpp>
0008 #include <Acts/Definitions/Units.hpp>
0009 #if Acts_VERSION_MAJOR >= 46
0010 #include <Acts/EventData/BoundTrackParameters.hpp>
0011 #else
0012 #include <Acts/EventData/TrackParameters.hpp>
0013 #endif
0014 #include <Acts/EventData/TrackProxy.hpp>
0015 #include <Acts/Propagator/EigenStepper.hpp>
0016 #include <Acts/Propagator/Propagator.hpp>
0017 #include <Acts/Propagator/VoidNavigator.hpp>
0018 #include <Acts/Surfaces/Surface.hpp>
0019 #include <Acts/Utilities/Logger.hpp>
0020 #include <Acts/Utilities/Result.hpp>
0021 #include <Acts/Vertexing/FullBilloirVertexFitter.hpp>
0022 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0023 #include <Acts/Vertexing/IVertexFinder.hpp>
0024 #include <Acts/Vertexing/ImpactPointEstimator.hpp>
0025 #include <Acts/Vertexing/IterativeVertexFinder.hpp>
0026 #include <Acts/Vertexing/LinearizedTrack.hpp>
0027 #include <Acts/Vertexing/TrackAtVertex.hpp>
0028 #include <Acts/Vertexing/Vertex.hpp>
0029 #include <Acts/Vertexing/VertexingOptions.hpp>
0030 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
0031 #include <ActsExamples/EventData/Track.hpp>
0032 #include <edm4eic/Cov4f.h>
0033 #include <edm4eic/ReconstructedParticleCollection.h>
0034 #include <edm4eic/Track.h>
0035 #include <edm4eic/TrackParameters.h>
0036 #include <edm4eic/Trajectory.h>
0037 #include <edm4eic/unit_system.h>
0038 #include <edm4hep/Vector2f.h>
0039 #include <edm4hep/Vector4f.h>
0040 #include <podio/RelationRange.h>
0041 #include <spdlog/common.h>
0042 #include <cmath>
0043 #include <string>
0044 #include <tuple>
0045 #include <utility>
0046 #include <vector>
0047
0048 #include "ActsGeometryProvider.h"
0049 #include "algorithms/tracking/IterativeVertexFinderConfig.h"
0050 #include "extensions/spdlog/SpdlogToActs.h"
0051
0052 void eicrecon::IterativeVertexFinder::process(const Input& input, const Output& output) const {
0053   const auto [trackStates, tracks, reconParticles] = input;
0054   auto [outputVertices]                            = output;
0055
0056   // Construct ConstTrackContainer from underlying containers
0057   auto trackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(*trackStates);
0058   auto trackContainer      = std::make_shared<Acts::ConstVectorTrackContainer>(*tracks);
0059   ActsExamples::ConstTrackContainer constTracks(trackContainer, trackStateContainer);
0060
0061   using Propagator           = Acts::Propagator<Acts::EigenStepper<>>;
0062   using Linearizer           = Acts::HelicalTrackLinearizer;
0063   using VertexFitter         = Acts::FullBilloirVertexFitter;
0064   using ImpactPointEstimator = Acts::ImpactPointEstimator;
0065   using VertexSeeder         = Acts::ZScanVertexFinder;
0066   using VertexFinder         = Acts::IterativeVertexFinder;
0067   using VertexFinderOptions  = Acts::VertexingOptions;
0068
0069   // Convert algorithm log level to Acts log level
0070   const auto spdlog_level = static_cast<spdlog::level::level_enum>(this->level());
0071   const auto acts_level   = eicrecon::SpdlogToActsLevel(spdlog_level);
0072   ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("IVF", acts_level));
0073
0074   Acts::EigenStepper<> stepper(m_BField);
0075
0076   // Set up propagator with void navigator
0077   auto propagator = std::make_shared<Propagator>(stepper, Acts::VoidNavigator{},
0078                                                  logger().cloneWithSuffix("Prop"));
0079
0080   // Setup the track linearizer
0081   Linearizer::Config linearizerCfg;
0082   linearizerCfg.bField     = m_BField;
0083   linearizerCfg.propagator = propagator;
0084   Linearizer linearizer(linearizerCfg, logger().cloneWithSuffix("HelLin"));
0085
0086   // Setup the vertex fitter
0087   VertexFitter::Config vertexFitterCfg;
0088   vertexFitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0089   vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0090   VertexFitter vertexFitter(vertexFitterCfg);
0091
0092   // Setup the seed finder
0093   ImpactPointEstimator::Config ipEstCfg(m_BField, propagator);
0094   ImpactPointEstimator ipEst(ipEstCfg);
0095   VertexSeeder::Config seederCfg(ipEst);
0096   seederCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0097   auto seeder = std::make_shared<VertexSeeder>(seederCfg);
0098
0099   // Set up the actual vertex finder
0100   VertexFinder::Config finderCfg(std::move(vertexFitter), std::move(seeder), std::move(ipEst));
0101   finderCfg.maxVertices                 = m_cfg.maxVertices;
0102   finderCfg.reassignTracksAfterFirstFit = m_cfg.reassignTracksAfterFirstFit;
0103   finderCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0104   finderCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0105   finderCfg.field = m_BField;
0106   VertexFinder finder(std::move(finderCfg));
0107
0108   // Get run-scoped contexts from service
0109   const auto& gctx = m_geoSvc->getActsGeometryContext();
0110   const auto& mctx = m_geoSvc->getActsMagneticFieldContext();
0111
0112   Acts::IVertexFinder::State state(std::in_place_type<VertexFinder::State>, *m_BField, mctx);
0113
0114   VertexFinderOptions finderOpts(gctx, mctx);
0115
0116   std::vector<Acts::InputTrack> inputTracks;
0117   std::vector<Acts::BoundTrackParameters> trackParameters;
0118   trackParameters.reserve(constTracks.size());
0119
0120   for (const auto& track : constTracks) {
0121     // Filter tracks based on minimum number of measurements (hits)
0122     if (track.nMeasurements() < m_cfg.minTrackHits) {
0123       trace("Track rejected: {} measurements < {} minimum required measurements",
0124             track.nMeasurements(), m_cfg.minTrackHits);
0125       continue;
0126     }
0127
0128     // Create BoundTrackParameters and store it
0129     trackParameters.emplace_back(track.referenceSurface().getSharedPtr(), track.parameters(),
0130                                  track.covariance(), track.particleHypothesis());
0131     // Create InputTrack from stored parameters
0132     inputTracks.emplace_back(&trackParameters.back());
0133     trace("Track local position at input = {} mm, {} mm",
0134           track.parameters()[Acts::eBoundLoc0] / Acts::UnitConstants::mm,
0135           track.parameters()[Acts::eBoundLoc1] / Acts::UnitConstants::mm);
0136   }
0137
0138   std::vector<Acts::Vertex> vertices;
0139   auto result = finder.find(inputTracks, finderOpts, state);
0140   if (result.ok()) {
0141     vertices = std::move(result.value());
0142   }
0143
0144   for (const auto& vtx : vertices) {
0145     static constexpr auto mm_to_mm =
0146         static_cast<float>(edm4eic::unit::mm / Acts::UnitConstants::mm);
0147     static constexpr auto ns_to_ns =
0148         static_cast<float>(edm4eic::unit::ns / Acts::UnitConstants::ns);
0149     static constexpr float mm2_to_mm2     = mm_to_mm * mm_to_mm;
0150     static constexpr float ns2_to_ns2     = ns_to_ns * ns_to_ns;
0151     static constexpr float mm_ns_to_mm_ns = mm_to_mm * ns_to_ns;
0152     edm4eic::Cov4f cov(
0153         vtx.fullCovariance()(0, 0) * mm2_to_mm2, vtx.fullCovariance()(1, 1) * mm2_to_mm2,
0154         vtx.fullCovariance()(2, 2) * mm2_to_mm2, vtx.fullCovariance()(3, 3) * ns2_to_ns2,
0155         vtx.fullCovariance()(0, 1) * mm2_to_mm2, vtx.fullCovariance()(0, 2) * mm2_to_mm2,
0156         vtx.fullCovariance()(0, 3) * mm_ns_to_mm_ns, vtx.fullCovariance()(1, 2) * mm2_to_mm2,
0157         vtx.fullCovariance()(1, 3) * mm_ns_to_mm_ns, vtx.fullCovariance()(2, 3) * mm_ns_to_mm_ns);
0158     auto eicvertex = outputVertices->create();
0159     eicvertex.setType(1); // boolean flag if vertex is primary vertex of event
0160     eicvertex.setChi2((float)vtx.fitQuality().first); // chi2
0161     eicvertex.setNdf((float)vtx.fitQuality().second); // ndf
0162     eicvertex.setPosition({
0163         static_cast<float>(vtx.position().x() * mm_to_mm),
0164         static_cast<float>(vtx.position().y() * mm_to_mm),
0165         static_cast<float>(vtx.position().z() * mm_to_mm),
0166         static_cast<float>(vtx.time() * ns_to_ns),
0167     });                              // vtxposition
0168     eicvertex.setPositionError(cov); // covariance
0169
0170     for (const auto& t : vtx.tracks()) {
0171       const auto& par = Acts::InputTrack::extractParameters(t.originalParams);
0172       trace("Track local position from vertex = {} mm, {} mm",
0173             par.localPosition().x() / Acts::UnitConstants::mm,
0174             par.localPosition().y() / Acts::UnitConstants::mm);
0175       float loc_a = par.localPosition().x();
0176       float loc_b = par.localPosition().y();
0177
0178       for (const auto& part : *reconParticles) {
0179         const auto& tracks = part.getTracks();
0180         for (const auto& trk : tracks) {
0181           const auto& traj    = trk.getTrajectory();
0182           const auto& trkPars = traj.getTrackParameters();
0183           for (const auto& par : trkPars) {
0184             const double EPSILON = 1.0e-4; // mm
0185             if (std::abs((par.getLoc().a / edm4eic::unit::mm) - (loc_a / Acts::UnitConstants::mm)) <
0186                     EPSILON &&
0187                 std::abs((par.getLoc().b / edm4eic::unit::mm) - (loc_b / Acts::UnitConstants::mm)) <
0188                     EPSILON) {
0189               trace("From ReconParticles, track local position [Loc a, Loc b] = {} mm, {} mm",
0190                     par.getLoc().a / edm4eic::unit::mm, par.getLoc().b / edm4eic::unit::mm);
0191               eicvertex.addToAssociatedParticles(part);
0192             } // endif
0193           } // end for par
0194         } // end for trk
0195       } // end for part
0196     } // end for t
0197     debug("One vertex found at (x,y,z) = ({}, {}, {}) mm.",
0198           vtx.position().x() / Acts::UnitConstants::mm,
0199           vtx.position().y() / Acts::UnitConstants::mm,
0200           vtx.position().z() / Acts::UnitConstants::mm);
0201
0202   } // end for vtx
0203 }