algorithms/tracking/SecondaryVertexFinder.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Dongwi H. Dongwi (Bishoy)
0003
0004 #include "SecondaryVertexFinder.h"
0005
0006 #include <Acts/Definitions/TrackParametrization.hpp>
0007 #include <Acts/Definitions/Units.hpp>
0008 #include <Acts/EventData/GenericBoundTrackParameters.hpp>
0009 #include <Acts/EventData/TrackParameters.hpp>
0010 #include <Acts/EventData/TrackProxy.hpp>
0011 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0012 #include <Acts/EventData/VectorTrackContainer.hpp>
0013 #include <Acts/Propagator/EigenStepper.hpp>
0014 #include <Acts/Propagator/Propagator.hpp>
0015 #include <Acts/Propagator/VoidNavigator.hpp>
0016 #include <Acts/Surfaces/Surface.hpp>
0017 #include <Acts/Utilities/AnnealingUtility.hpp>
0018 #include <Acts/Utilities/Logger.hpp>
0019 #include <Acts/Utilities/Result.hpp>
0020 #include <Acts/Vertexing/AdaptiveGridTrackDensity.hpp>
0021 #include <Acts/Vertexing/IVertexFinder.hpp>
0022 #include <Acts/Vertexing/LinearizedTrack.hpp>
0023 #include <Acts/Vertexing/TrackAtVertex.hpp>
0024 #include <ActsExamples/EventData/Track.hpp>
0025 #include <edm4eic/Cov4f.h>
0026 #include <edm4eic/Track.h>
0027 #include <edm4eic/TrackParameters.h>
0028 #include <edm4eic/Trajectory.h>
0029 #include <edm4eic/unit_system.h>
0030 #include <edm4hep/Vector2f.h>
0031 #include <edm4hep/Vector4f.h>
0032 #include <podio/RelationRange.h>
0033 #include <spdlog/common.h>
0034 #include <cmath>
0035 #include <limits>
0036 #include <memory>
0037 #include <string>
0038 #include <tuple>
0039 #include <utility>
0040
0041 #include "ActsGeometryProvider.h"
0042 #include "SecondaryVertexFinderConfig.h"
0043 #include "extensions/spdlog/SpdlogToActs.h"
0044
0045 namespace eicrecon {
0046
0047 void SecondaryVertexFinder::process(const SecondaryVertexFinder::Input& input,
0048                                     const SecondaryVertexFinder::Output& output) const {
0049   auto [recotracks, trackStates, tracks] = input;
0050   auto [primaryVertices, outputVertices] = output;
0051
0052   Acts::EigenStepper<> stepperSec(m_BField);
0053
0054   //Need to make sure that the track container is not actually empty
0055   if (tracks->size_impl() > 0) {
0056     // Calculate primary vertex using AMVF
0057     calculatePrimaryVertex(*recotracks, trackStates, tracks, stepperSec, *primaryVertices);
0058     // Primary vertex collection container to be used in Sec. Vertex fitting
0059     calculateSecondaryVertex(*recotracks, trackStates, tracks, stepperSec, *outputVertices);
0060   }
0061 }
0062
0063 //Quickly calculate the PV using the Adaptive Multi-vertex Finder
0064 void SecondaryVertexFinder::calculatePrimaryVertex(
0065     const edm4eic::ReconstructedParticleCollection& reconParticles,
0066     const Acts::ConstVectorMultiTrajectory* trackStates,
0067     const Acts::ConstVectorTrackContainer* tracks, Acts::EigenStepper<> stepperSec,
0068     edm4eic::VertexCollection& prmVertices) const {
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("AMVF", acts_level));
0073
0074   // Geometry and field contexts
0075   const auto& gctx = m_geoSvc->getActsGeometryContext();
0076   const auto& mctx = m_geoSvc->getActsMagneticFieldContext();
0077
0078   // Set-up the propagator
0079   using PropagatorSec = Acts::Propagator<Acts::EigenStepper<>>;
0080
0081   // Set up propagator with void navigator
0082   auto propagatorSec = std::make_shared<PropagatorSec>(stepperSec, Acts::VoidNavigator{},
0083                                                        logger().cloneWithSuffix("Prop"));
0084
0085   // Set up track density used during vertex seeding
0086   Acts::AdaptiveGridTrackDensity::Config trkDensityConfig;
0087   // Bin extent in z-direction
0088   trkDensityConfig.spatialBinExtent = m_cfg.spatialBinExtent;
0089   // Bin extent in t-direction
0090   trkDensityConfig.temporalBinExtent = m_cfg.temporalBinExtent;
0091   trkDensityConfig.useTime           = m_cfg.useTime;
0092   Acts::AdaptiveGridTrackDensity trkDensity(trkDensityConfig);
0093
0094   // Setup the track linearizer
0095   LinearizerSec::Config linearizerConfigSec(m_BField, propagatorSec);
0096   std::unique_ptr<const Acts::Logger> linearizer_log = logger().cloneWithSuffix("Linearizer");
0097   LinearizerSec linearizerSec(linearizerConfigSec, std::move(linearizer_log));
0098
0099   //Staring multivertex fitter
0100   // Set up deterministic annealing with user-defined temperatures
0101   Acts::AnnealingUtility::Config annealingConfig;
0102   annealingConfig.setOfTemperatures = {9., 1.0};
0103   Acts::AnnealingUtility annealingUtility(annealingConfig);
0104   // Setup the vertex fitter
0105   ImpactPointEstimator::Config ipEstConfig(m_BField, propagatorSec);
0106   ImpactPointEstimator ipEst(ipEstConfig);
0107   VertexFitterSec::Config vertexFitterConfigSec(ipEst);
0108
0109   vertexFitterConfigSec.annealingTool     = annealingUtility;
0110   vertexFitterConfigSec.minWeight         = m_cfg.minWeight;
0111   vertexFitterConfigSec.maxDistToLinPoint = m_cfg.maxDistToLinPoint;
0112   vertexFitterConfigSec.doSmoothing       = m_cfg.doSmoothing;
0113   vertexFitterConfigSec.useTime           = m_cfg.useTime;
0114   vertexFitterConfigSec.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0115   vertexFitterConfigSec.trackLinearizer.connect<&LinearizerSec::linearizeTrack>(&linearizerSec);
0116   VertexFitterSec vertexFitterSec(std::move(vertexFitterConfigSec));
0117
0118   // Set up vertex seeder and finder
0119   using seedFinder = Acts::AdaptiveGridDensityVertexFinder;
0120   std::shared_ptr<const Acts::IVertexFinder> seeder;
0121   seedFinder::Config seederConfig(trkDensity);
0122   seederConfig.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0123   seeder = std::make_shared<seedFinder>(seedFinder::Config{seederConfig});
0124
0125   VertexFinderSec::Config vertexfinderConfigSec(std::move(vertexFitterSec), std::move(seeder),
0126                                                 std::move(ipEst), m_BField);
0127
0128   // The vertex finder state
0129   // Set the initial variance of the 4D vertex position. Since time is on a
0130   // numerical scale, we have to provide a greater value in the corresponding
0131   // dimension.
0132   vertexfinderConfigSec.initialVariances = m_cfg.initialVariances;
0133   //Use time for Sec. Vertex
0134   vertexfinderConfigSec.useTime                    = m_cfg.useTime;
0135   vertexfinderConfigSec.tracksMaxZinterval         = m_cfg.tracksMaxZinterval;
0136   vertexfinderConfigSec.maxIterations              = m_cfg.maxIterations;
0137   vertexfinderConfigSec.doFullSplitting            = m_cfg.doFullSplitting;
0138   vertexfinderConfigSec.tracksMaxSignificance      = m_cfg.tracksMaxSignificance;
0139   vertexfinderConfigSec.maxMergeVertexSignificance = m_cfg.maxMergeVertexSignificance;
0140
0141   if (m_cfg.useTime) {
0142     // When using time, we have an extra contribution to the chi2 by the time
0143     // coordinate.
0144     vertexfinderConfigSec.tracksMaxSignificance      = m_cfg.tracksMaxSignificance;
0145     vertexfinderConfigSec.maxMergeVertexSignificance = m_cfg.maxMergeVertexSignificance;
0146   }
0147
0148   vertexfinderConfigSec.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0149
0150   vertexfinderConfigSec.bField = m_BField;
0151   VertexFinderSec finder(std::move(vertexfinderConfigSec));
0152   // Instantiate the finder
0153   auto stateSec = finder.makeState(mctx);
0154
0155   VertexFinderOptionsSec vfOptions(gctx, mctx);
0156
0157   // Construct ConstTrackContainer from underlying containers
0158   auto trackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(*trackStates);
0159   auto trackContainer      = std::make_shared<Acts::ConstVectorTrackContainer>(*tracks);
0160   ActsExamples::ConstTrackContainer constTracks(trackContainer, trackStateContainer);
0161
0162   std::vector<Acts::InputTrack> inputTracks;
0163   std::vector<Acts::BoundTrackParameters> trackParameters;
0164   trackParameters.reserve(constTracks.size());
0165
0166   for (const auto& track : constTracks) {
0167     // Create BoundTrackParameters and store it
0168     trackParameters.emplace_back(track.referenceSurface().getSharedPtr(), track.parameters(),
0169                                  track.covariance(), track.particleHypothesis());
0170     // Create InputTrack from stored parameters
0171     inputTracks.emplace_back(&trackParameters.back());
0172     trace("Track local position at input = {} mm, {} mm",
0173           track.parameters()[Acts::eBoundLoc0] / Acts::UnitConstants::mm,
0174           track.parameters()[Acts::eBoundLoc1] / Acts::UnitConstants::mm);
0175   }
0176
0177   std::vector<Acts::Vertex> vertices;
0178   auto result = finder.find(inputTracks, vfOptions, stateSec);
0179   if (result.ok()) {
0180     vertices = std::move(result.value());
0181   }
0182
0183   for (const auto& vtx : vertices) {
0184     edm4eic::Cov4f cov(vtx.fullCovariance()(0, 0), vtx.fullCovariance()(1, 1),
0185                        vtx.fullCovariance()(2, 2), vtx.fullCovariance()(3, 3),
0186                        vtx.fullCovariance()(0, 1), vtx.fullCovariance()(0, 2),
0187                        vtx.fullCovariance()(0, 3), vtx.fullCovariance()(1, 2),
0188                        vtx.fullCovariance()(1, 3), vtx.fullCovariance()(2, 3));
0189     auto eicvertex = prmVertices.create();
0190     eicvertex.setType(1); // boolean flag if vertex is primary vertex of event
0191     eicvertex.setChi2(static_cast<float>(vtx.fitQuality().first)); // chi2
0192     eicvertex.setNdf(static_cast<float>(vtx.fitQuality().second)); // ndf
0193     eicvertex.setPosition({
0194         static_cast<float>(vtx.position().x()),
0195         static_cast<float>(vtx.position().y()),
0196         static_cast<float>(vtx.position().z()),
0197         static_cast<float>(vtx.time()),
0198     });                              // vtxposition
0199     eicvertex.setPositionError(cov); // covariance
0200
0201     for (const auto& t : vtx.tracks()) {
0202       const auto& par = vertexfinderConfigSec.extractParameters(t.originalParams);
0203       trace("Track local position from vertex = {} mm, {} mm",
0204             par.localPosition().x() / Acts::UnitConstants::mm,
0205             par.localPosition().y() / Acts::UnitConstants::mm);
0206       float loc_a = par.localPosition().x();
0207       float loc_b = par.localPosition().y();
0208
0209       for (const auto& part : reconParticles) {
0210         const auto& tracks = part.getTracks();
0211         for (const auto& trk : tracks) {
0212           const auto& traj    = trk.getTrajectory();
0213           const auto& trkPars = traj.getTrackParameters();
0214           for (const auto& par : trkPars) {
0215             double EPSILON = std::numeric_limits<double>::epsilon(); // mm
0216             if (std::abs((par.getLoc().a / edm4eic::unit::mm) - (loc_a / Acts::UnitConstants::mm)) <
0217                     EPSILON &&
0218                 std::abs((par.getLoc().b / edm4eic::unit::mm) - (loc_b / Acts::UnitConstants::mm)) <
0219                     EPSILON) {
0220               trace("From ReconParticles, track local position [Loc a, Loc b] = {} mm, {} mm",
0221                     par.getLoc().a / edm4eic::unit::mm, par.getLoc().b / edm4eic::unit::mm);
0222               eicvertex.addToAssociatedParticles(part);
0223             } // endif
0224           } // end for par
0225         } // end for trk
0226       } // end for part
0227     } // end for t
0228     debug("One AMVF vertex found at (x,y,z) = ({}, {}, {}) mm.",
0229           vtx.position().x() / Acts::UnitConstants::mm,
0230           vtx.position().y() / Acts::UnitConstants::mm,
0231           vtx.position().z() / Acts::UnitConstants::mm);
0232   } // end for vtx
0233 }
0234
0235 void SecondaryVertexFinder::calculateSecondaryVertex(
0236     const edm4eic::ReconstructedParticleCollection& reconParticles,
0237     const Acts::ConstVectorMultiTrajectory* trackStates,
0238     const Acts::ConstVectorTrackContainer* tracks, Acts::EigenStepper<> stepperSec,
0239     edm4eic::VertexCollection& secVertices) const {
0240   // Convert algorithm log level to Acts log level
0241   const auto spdlog_level = static_cast<spdlog::level::level_enum>(this->level());
0242   const auto acts_level   = eicrecon::SpdlogToActsLevel(spdlog_level);
0243   ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("AMVF", acts_level));
0244
0245   // Geometry and field contexts
0246   const auto& gctx = m_geoSvc->getActsGeometryContext();
0247   const auto& mctx = m_geoSvc->getActsMagneticFieldContext();
0248
0249   // Set-up the propagator
0250   using PropagatorSec = Acts::Propagator<Acts::EigenStepper<>>;
0251
0252   // Set up propagator with void navigator
0253   auto propagatorSec = std::make_shared<PropagatorSec>(stepperSec, Acts::VoidNavigator{},
0254                                                        logger().cloneWithSuffix("Prop"));
0255
0256   // Set up track density used during vertex seeding
0257   Acts::AdaptiveGridTrackDensity::Config trkDensityConfig;
0258   // Bin extent in z-direction
0259   trkDensityConfig.spatialBinExtent = m_cfg.spatialBinExtent;
0260   // Bin extent in t-direction
0261   trkDensityConfig.temporalBinExtent = m_cfg.temporalBinExtent;
0262   trkDensityConfig.useTime           = m_cfg.useTime;
0263   Acts::AdaptiveGridTrackDensity trkDensity(trkDensityConfig);
0264
0265   // Setup the track linearizer
0266   LinearizerSec::Config linearizerConfigSec(m_BField, propagatorSec);
0267   std::unique_ptr<const Acts::Logger> linearizer_log = logger().cloneWithSuffix("Linearizer");
0268   LinearizerSec linearizerSec(linearizerConfigSec, std::move(linearizer_log));
0269
0270   //Staring multivertex fitter
0271   // Set up deterministic annealing with user-defined temperatures
0272   Acts::AnnealingUtility::Config annealingConfig;
0273   annealingConfig.setOfTemperatures = {9., 1.0};
0274   Acts::AnnealingUtility annealingUtility(annealingConfig);
0275   // Setup the vertex fitter
0276   ImpactPointEstimator::Config ipEstConfig(m_BField, propagatorSec);
0277   ImpactPointEstimator ipEst(ipEstConfig);
0278   VertexFitterSec::Config vertexFitterConfigSec(ipEst);
0279
0280   vertexFitterConfigSec.annealingTool     = annealingUtility;
0281   vertexFitterConfigSec.minWeight         = m_cfg.minWeight;
0282   vertexFitterConfigSec.maxDistToLinPoint = m_cfg.maxDistToLinPoint;
0283   vertexFitterConfigSec.doSmoothing       = m_cfg.doSmoothing;
0284   vertexFitterConfigSec.useTime           = m_cfg.useTime;
0285   vertexFitterConfigSec.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0286   vertexFitterConfigSec.trackLinearizer.connect<&LinearizerSec::linearizeTrack>(&linearizerSec);
0287   VertexFitterSec vertexFitterSec(std::move(vertexFitterConfigSec));
0288
0289   // Set up vertex seeder and finder
0290   using seedFinder = Acts::AdaptiveGridDensityVertexFinder;
0291   std::shared_ptr<const Acts::IVertexFinder> seeder;
0292   seedFinder::Config seederConfig(trkDensity);
0293   seederConfig.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0294   seeder = std::make_shared<seedFinder>(seedFinder::Config{seederConfig});
0295
0296   VertexFinderSec::Config vertexfinderConfigSec(std::move(vertexFitterSec), std::move(seeder),
0297                                                 std::move(ipEst), m_BField);
0298
0299   // The vertex finder state
0300   // Set the initial variance of the 4D vertex position. Since time is on a
0301   // numerical scale, we have to provide a greater value in the corresponding
0302   // dimension.
0303   vertexfinderConfigSec.initialVariances = m_cfg.initialVariances;
0304   //Use time for Sec. Vertex
0305   vertexfinderConfigSec.useTime                    = m_cfg.useTime;
0306   vertexfinderConfigSec.useSeedConstraint          = m_cfg.useSeedConstraint;
0307   vertexfinderConfigSec.tracksMaxZinterval         = m_cfg.tracksMaxZintervalSec;
0308   vertexfinderConfigSec.maxIterations              = m_cfg.maxSecIterations;
0309   vertexfinderConfigSec.doFullSplitting            = m_cfg.doFullSplitting;
0310   vertexfinderConfigSec.tracksMaxSignificance      = m_cfg.tracksMaxSignificance;
0311   vertexfinderConfigSec.maxMergeVertexSignificance = m_cfg.maxMergeVertexSignificance;
0312
0313   if (m_cfg.useTime) {
0314     // When using time, we have an extra contribution to the chi2 by the time
0315     // coordinate.
0316     vertexfinderConfigSec.tracksMaxSignificance      = m_cfg.tracksMaxSignificance;
0317     vertexfinderConfigSec.maxMergeVertexSignificance = m_cfg.maxMergeVertexSignificance;
0318   }
0319
0320   vertexfinderConfigSec.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0321
0322   vertexfinderConfigSec.bField = m_BField;
0323   VertexFinderSec finder(std::move(vertexfinderConfigSec));
0324   // Instantiate the finder
0325   auto stateSec = finder.makeState(mctx);
0326
0327   VertexFinderOptionsSec vfOptions(gctx, mctx);
0328
0329   //--->Add Prm Vertex container here
0330   // Construct ConstTrackContainer from underlying containers
0331   auto trackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(*trackStates);
0332   auto trackContainer      = std::make_shared<Acts::ConstVectorTrackContainer>(*tracks);
0333   ActsExamples::ConstTrackContainer constTracks(trackContainer, trackStateContainer);
0334
0335   // Build BoundTrackParameters for all tracks upfront
0336   std::vector<Acts::BoundTrackParameters> allTrackParameters;
0337   allTrackParameters.reserve(constTracks.size());
0338   for (const auto& track : constTracks) {
0339     allTrackParameters.emplace_back(track.referenceSurface().getSharedPtr(), track.parameters(),
0340                                     track.covariance(), track.particleHypothesis());
0341   }
0342
0343   std::vector<Acts::InputTrack> inputTracks;
0344   for (unsigned int i = 0; i < allTrackParameters.size() - 1; i++) {
0345     for (unsigned int j = i + 1; j < allTrackParameters.size(); j++) {
0346       inputTracks.emplace_back(&allTrackParameters[i]);
0347       inputTracks.emplace_back(&allTrackParameters[j]);
0348       // run the vertex finder for both tracks
0349       std::vector<Acts::Vertex> verticesSec;
0350       auto resultSecondary = finder.find(inputTracks, vfOptions, stateSec);
0351       if (resultSecondary.ok()) {
0352         verticesSec = std::move(resultSecondary.value());
0353       }
0354
0355       for (const auto& secvertex : verticesSec) {
0356         edm4eic::Cov4f cov(secvertex.fullCovariance()(0, 0), secvertex.fullCovariance()(1, 1),
0357                            secvertex.fullCovariance()(2, 2), secvertex.fullCovariance()(3, 3),
0358                            secvertex.fullCovariance()(0, 1), secvertex.fullCovariance()(0, 2),
0359                            secvertex.fullCovariance()(0, 3), secvertex.fullCovariance()(1, 2),
0360                            secvertex.fullCovariance()(1, 3), secvertex.fullCovariance()(2, 3));
0361         auto eicvertex = secVertices.create();
0362         eicvertex.setType(0); // boolean flag if vertex is primary vertex of event
0363         eicvertex.setChi2(static_cast<float>(secvertex.fitQuality().first)); // chi2
0364         eicvertex.setNdf(static_cast<float>(secvertex.fitQuality().second)); // ndf
0365         eicvertex.setPosition({
0366             static_cast<float>(secvertex.position().x()),
0367             static_cast<float>(secvertex.position().y()),
0368             static_cast<float>(secvertex.position().z()),
0369             static_cast<float>(secvertex.time()),
0370         });                              // vtxposition
0371         eicvertex.setPositionError(cov); // covariance
0372
0373         for (const auto& t : secvertex.tracks()) {
0374           const auto& par = vertexfinderConfigSec.extractParameters(t.originalParams);
0375           trace("Track local position from vertex = {} mm, {} mm",
0376                 par.localPosition().x() / Acts::UnitConstants::mm,
0377                 par.localPosition().y() / Acts::UnitConstants::mm);
0378           float loc_a = par.localPosition().x();
0379           float loc_b = par.localPosition().y();
0380
0381           for (const auto& part : reconParticles) {
0382             const auto& tracks = part.getTracks();
0383             for (const auto& trk : tracks) {
0384               const auto& traj    = trk.getTrajectory();
0385               const auto& trkPars = traj.getTrackParameters();
0386               for (const auto& par : trkPars) {
0387                 double EPSILON = std::numeric_limits<double>::epsilon(); // mm
0388                 if (std::abs((par.getLoc().a / edm4eic::unit::mm) -
0389                              (loc_a / Acts::UnitConstants::mm)) < EPSILON &&
0390                     std::abs((par.getLoc().b / edm4eic::unit::mm) -
0391                              (loc_b / Acts::UnitConstants::mm)) < EPSILON) {
0392                   trace("From ReconParticles, track local position [Loc a, Loc b] = {} mm, {} mm",
0393                         par.getLoc().a / edm4eic::unit::mm, par.getLoc().b / edm4eic::unit::mm);
0394                   eicvertex.addToAssociatedParticles(part);
0395                 } // endif
0396               } // end for par
0397             } // end for trk
0398           } // end for part
0399         } // end for t
0400       } // end for secvertex
0401
0402       // empty the vector for the next set of tracks
0403       inputTracks.clear();
0404     } //end of int j=i+1
0405   } // end of int i=0; i<constTracks->size()
0406 }
0407
0408 } // namespace eicrecon