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