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/EventData/GenericBoundTrackParameters.hpp>
0011 #include <Acts/EventData/GenericParticleHypothesis.hpp>
0012 #include <Acts/EventData/ParticleHypothesis.hpp>
0013 #include <Acts/EventData/TrackParameters.hpp>
0014 #include <Acts/Propagator/EigenStepper.hpp>
0015 #include <Acts/Propagator/Propagator.hpp>
0016 #include <Acts/Propagator/detail/VoidPropagatorComponents.hpp>
0017 #include <Acts/Utilities/Logger.hpp>
0018 #include <Acts/Utilities/Result.hpp>
0019 #include <Acts/Utilities/VectorHelpers.hpp>
0020 #include <Acts/Vertexing/FullBilloirVertexFitter.hpp>
0021 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0022 #include <Acts/Vertexing/ImpactPointEstimator.hpp>
0023 #include <Acts/Vertexing/IterativeVertexFinder.hpp>
0024 #include <Acts/Vertexing/Vertex.hpp>
0025 #include <Acts/Vertexing/VertexingOptions.hpp>
0026 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
0027 #include <ActsExamples/EventData/Trajectories.hpp>
0028 #include <boost/container/vector.hpp>
0029 #include <edm4eic/Cov4f.h>
0030 #include <math.h>
0031 #include <Eigen/Core>
0032 #include <Eigen/Geometry>
0033 #include <Eigen/LU>
0034 #include <algorithm>
0035 #include <optional>
0036 #include <tuple>
0037 #include <utility>
0038
0039 #include "extensions/spdlog/SpdlogToActs.h"
0040
0041 void eicrecon::IterativeVertexFinder::init(std::shared_ptr<const ActsGeometryProvider> geo_svc,
0042                                            std::shared_ptr<spdlog::logger> log) {
0043
0044   m_log = log;
0045
0046   m_geoSvc = geo_svc;
0047
0048   m_BField =
0049       std::dynamic_pointer_cast<const eicrecon::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
0050   m_fieldctx = eicrecon::BField::BFieldVariant(m_BField);
0051 }
0052
0053 std::unique_ptr<edm4eic::VertexCollection> eicrecon::IterativeVertexFinder::produce(
0054     std::vector<const ActsExamples::Trajectories*> trajectories) {
0055
0056   auto outputVertices = std::make_unique<edm4eic::VertexCollection>();
0057
0058   using Propagator        = Acts::Propagator<Acts::EigenStepper<>>;
0059   using PropagatorOptions = Acts::PropagatorOptions<>;
0060   using Linearizer        = Acts::HelicalTrackLinearizer<Propagator>;
0061   using VertexFitter      = Acts::FullBilloirVertexFitter<Acts::BoundTrackParameters, Linearizer>;
0062   using ImpactPointEstimator = Acts::ImpactPointEstimator<Acts::BoundTrackParameters, Propagator>;
0063   using VertexSeeder         = Acts::ZScanVertexFinder<VertexFitter>;
0064   using VertexFinder         = Acts::IterativeVertexFinder<VertexFitter, VertexSeeder>;
0065   using VertexFinderOptions  = Acts::VertexingOptions<Acts::BoundTrackParameters>;
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>(
0073     stepper, Acts::detail::VoidNavigator{}, logger().cloneWithSuffix("Prop"));
0074   Acts::PropagatorOptions opts(m_geoctx, m_fieldctx);
0075
0076   // Setup the vertex fitter
0077   VertexFitter::Config vertexFitterCfg;
0078   VertexFitter vertexFitter(vertexFitterCfg);
0079   // Setup the track linearizer
0080   Linearizer::Config linearizerCfg(m_BField, propagator);
0081   Linearizer linearizer(linearizerCfg, logger().cloneWithSuffix("HelLin"));
0082   // Setup the seed finder
0083   ImpactPointEstimator::Config ipEstCfg(m_BField, propagator);
0084   ImpactPointEstimator ipEst(ipEstCfg);
0085   VertexSeeder::Config seederCfg(ipEst);
0086   VertexSeeder seeder(seederCfg);
0087   // Set up the actual vertex finder
0088   VertexFinder::Config finderCfg(std::move(vertexFitter), std::move(linearizer),
0089                                  std::move(seeder), std::move(ipEst));
0090   finderCfg.maxVertices                 = m_cfg.maxVertices;
0091   finderCfg.reassignTracksAfterFirstFit = m_cfg.reassignTracksAfterFirstFit;
0092   #if Acts_VERSION_MAJOR >= 31
0093   VertexFinder finder(std::move(finderCfg));
0094   #else
0095   VertexFinder finder(finderCfg);
0096   #endif
0097   VertexFinder::State state(*m_BField, m_fieldctx);
0098   VertexFinderOptions finderOpts(m_geoctx, m_fieldctx);
0099
0100   std::vector<const Acts::BoundTrackParameters*> inputTrackPointers;
0101
0102   for (const auto& trajectory : trajectories) {
0103     auto tips = trajectory->tips();
0104     if (tips.empty()) {
0105       continue;
0106     }
0107     /// CKF can provide multiple track trajectories for a single input seed
0108     for (auto& tip : tips) {
0109       inputTrackPointers.push_back(&(trajectory->trackParameters(tip)));
0110     }
0111   }
0112
0113   std::vector<Acts::Vertex<Acts::BoundTrackParameters>> vertices;
0114   auto result = finder.find(inputTrackPointers, finderOpts, state);
0115   if (result.ok()) {
0116     vertices = std::move(result.value());
0117   }
0118
0119   for (const auto& vtx : vertices) {
0120     edm4eic::Cov4f cov(vtx.fullCovariance()(0,0), vtx.fullCovariance()(1,1), vtx.fullCovariance()(2,2), vtx.fullCovariance()(3,3),
0121                        vtx.fullCovariance()(0,1), vtx.fullCovariance()(0,2), vtx.fullCovariance()(0,3),
0122                        vtx.fullCovariance()(1,2), vtx.fullCovariance()(1,3),
0123                        vtx.fullCovariance()(2,3));
0124     auto eicvertex = outputVertices->create();
0125     eicvertex.setType(1);                                  // boolean flag if vertex is primary vertex of event
0126     eicvertex.setChi2((float)vtx.fitQuality().first);      // chi2
0127     eicvertex.setNdf((float)vtx.fitQuality().second);      // ndf
0128     eicvertex.setPosition({
0129          (float)vtx.position().x(),
0130          (float)vtx.position().y(),
0131          (float)vtx.position().z(),
0132          (float)vtx.time(),
0133     }); // vtxposition
0134     eicvertex.setPositionError(cov);                          // covariance
0135   }
0136
0137   return std::move(outputVertices);
0138 }