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