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