diff --git a/CI/physmon/phys_perf_mon.sh b/CI/physmon/phys_perf_mon.sh index 82308f1a906..81288e902a9 100755 --- a/CI/physmon/phys_perf_mon.sh +++ b/CI/physmon/phys_perf_mon.sh @@ -4,14 +4,13 @@ set -e set -x -mode=${1:all} +mode=${1:-all} if ! [[ $mode = @(all|kalman|gsf|fullchains|vertexing|simulation) ]]; then echo "Usage: $0 (outdir)" exit 1 fi -outdir=${2:physmon} -[ -z "$outdir" ] && outdir=physmon +outdir=${2:-physmon} mkdir -p $outdir refdir=CI/physmon/reference diff --git a/CI/physmon/reference/performance_ivf_truth_estimated_hist.root b/CI/physmon/reference/performance_ivf_truth_estimated_hist.root index 1d2048cf09e..379d5bddde4 100644 Binary files a/CI/physmon/reference/performance_ivf_truth_estimated_hist.root and b/CI/physmon/reference/performance_ivf_truth_estimated_hist.root differ diff --git a/CMakeLists.txt b/CMakeLists.txt index 26e3183f8bb..80b8987ccae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,8 @@ option(ACTS_BUILD_PLUGIN_DD4HEP "Build DD4hep plugin" OFF) option(ACTS_BUILD_PLUGIN_PODIO "Build Podio plugin" OFF) option(ACTS_BUILD_PLUGIN_EDM4HEP "Build EDM4hep plugin" OFF) option(ACTS_BUILD_PLUGIN_FPEMON "Build FPE monitoring plugin" OFF) +option(ACTS_BUILD_PLUGIN_GEOMODEL "Build GeoModel plugin" OFF) +option(ACTS_USE_SYSTEM_GEOMODEL "Use a system-provided GeoModel installation" ${ACTS_USE_SYSTEM_LIBS}) option(ACTS_BUILD_PLUGIN_GEANT4 "Build Geant4 plugin" OFF) option(ACTS_BUILD_PLUGIN_EXATRKX "Build the Exa.TrkX plugin" OFF) option(ACTS_EXATRKX_ENABLE_ONNX "Build the Onnx backend for the exatrkx plugin" OFF) @@ -186,6 +188,7 @@ set(_acts_autodiff_version 0.6) set(_acts_boost_version 1.71.0) set(_acts_dd4hep_version 1.21) set(_acts_edm4hep_version 0.7) +set(_acts_geomodel_version 4.6.0) set(_acts_podio_version 0.16) set(_acts_doxygen_version 1.9.4) set(_acts_eigen3_version 3.3.7) @@ -287,6 +290,15 @@ if(ACTS_BUILD_PLUGIN_DD4HEP) find_package(Python 3.8 REQUIRED COMPONENTS Interpreter) find_package(DD4hep ${_acts_dd4hep_version} REQUIRED CONFIG COMPONENTS DDCore DDDetectors) endif() +if(ACTS_BUILD_PLUGIN_GEOMODEL) + if(ACTS_USE_SYSTEM_GEOMODEL) + message(STATUS "Using system installation of GeoModel") + find_package(GeoModel ${_acts_geomodel_version} REQUIRED CONFIG) + else() + add_subdirectory(thirdparty/GeoModel) + endif() +endif() + if(ACTS_BUILD_PLUGIN_JSON) if(ACTS_USE_SYSTEM_NLOHMANN_JSON) message(STATUS "Using system installation of nlohmann::json") diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index 2ac69d01991..17c194a25ea 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -405,7 +405,6 @@ class Gx2Fitter { // Update: // - Waiting for a current surface auto surface = navigator.currentSurface(state.navigation); - // std::string direction = state.stepping.navDir.toString(); if (surface != nullptr && surface->associatedDetectorElement() != nullptr) { ++result.surfaceCount; diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp index db9933cf87a..83153b627ae 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp @@ -31,18 +31,16 @@ namespace Acts { /// Unlike the GridDensityVertexFinder, this seeder implements an adaptive /// version where the density grid grows bigger with added tracks. class AdaptiveGridDensityVertexFinder final : public IVertexFinder { - using GridDensity = AdaptiveGridTrackDensity; - public: - using DensityMap = typename GridDensity::DensityMap; + using DensityMap = AdaptiveGridTrackDensity::DensityMap; /// @brief The Config struct struct Config { ///@param gDensity The grid density - Config(const GridDensity& gDensity) : gridDensity(gDensity) {} + Config(const AdaptiveGridTrackDensity& gDensity) : gridDensity(gDensity) {} // The grid density object - GridDensity gridDensity; + AdaptiveGridTrackDensity gridDensity; // Cache the main grid and the density contributions (trackGrid and z-bin) // for every single track. @@ -129,5 +127,3 @@ class AdaptiveGridDensityVertexFinder final : public IVertexFinder { }; } // namespace Acts - -#include "AdaptiveGridDensityVertexFinder.ipp" diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp index da84374a69b..7dbc7d7783a 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp @@ -15,6 +15,7 @@ #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/AMVFInfo.hpp" +#include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp" #include "Acts/Vertexing/IVertexFinder.hpp" #include "Acts/Vertexing/ImpactPointEstimator.hpp" #include "Acts/Vertexing/TrackLinearizer.hpp" @@ -25,19 +26,9 @@ namespace Acts { /// @brief Implements an iterative vertex finder -/// -//////////////////////////////////////////////////////////// -/// -/// Brief description of the algorithm implemented: -/// TODO -/// -//////////////////////////////////////////////////////////// -/// -/// @tparam vfitter_t Vertex fitter type -/// @tparam sfinder_t Seed finder type -template class AdaptiveMultiVertexFinder final : public IVertexFinder { - using FitterState_t = typename vfitter_t::State; + using VertexFitter = AdaptiveMultiVertexFitter; + using VertexFitterState = VertexFitter::State; public: /// Configuration struct @@ -48,7 +39,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @param sfinder The seed finder /// @param ipEst ImpactPointEstimator /// @param bIn Input magnetic field - Config(vfitter_t fitter, std::shared_ptr sfinder, + Config(VertexFitter fitter, std::shared_ptr sfinder, ImpactPointEstimator ipEst, std::shared_ptr bIn) : vertexFitter(std::move(fitter)), @@ -57,7 +48,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { bField{std::move(bIn)} {} // Vertex fitter - vfitter_t vertexFitter; + VertexFitter vertexFitter; // Vertex seed finder std::shared_ptr seedFinder; @@ -230,7 +221,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// from seed track collection in last iteration /// /// @return The seed vertex - Result doSeeding( + Result> doSeeding( const std::vector& trackVector, Vertex& currentConstraint, const VertexingOptions& vertexingOptions, IVertexFinder::State& seedFinderState, @@ -264,7 +255,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @param vertexingOptions Vertexing options Result addCompatibleTracksToVertex( const std::vector& tracks, Vertex& vtx, - FitterState_t& fitterState, + VertexFitterState& fitterState, const VertexingOptions& vertexingOptions) const; /// @brief Method that tries to recover from cases where no tracks @@ -282,7 +273,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { Result canRecoverFromNoCompatibleTracks( const std::vector& allTracks, const std::vector& seedTracks, Vertex& vtx, - const Vertex& currentConstraint, FitterState_t& fitterState, + const Vertex& currentConstraint, VertexFitterState& fitterState, const VertexingOptions& vertexingOptions) const; /// @brief Method that tries to prepare the vertex for the fit @@ -299,7 +290,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { Result canPrepareVertexForFit( const std::vector& allTracks, const std::vector& seedTracks, Vertex& vtx, - const Vertex& currentConstraint, FitterState_t& fitterState, + const Vertex& currentConstraint, VertexFitterState& fitterState, const VertexingOptions& vertexingOptions) const; /// @brief Method that checks if vertex is a good vertex and if @@ -313,7 +304,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @return pair(nCompatibleTracks, isGoodVertex) std::pair checkVertexAndCompatibleTracks( Vertex& vtx, const std::vector& seedTracks, - FitterState_t& fitterState, bool useVertexConstraintInFit) const; + VertexFitterState& fitterState, bool useVertexConstraintInFit) const; /// @brief Method that removes all tracks that are compatible with /// current vertex from seedTracks @@ -325,7 +316,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// removed void removeCompatibleTracksFromSeedTracks( Vertex& vtx, std::vector& seedTracks, - FitterState_t& fitterState, + VertexFitterState& fitterState, std::vector& removedSeedTracks) const; /// @brief Method that tries to remove an incompatible track @@ -341,7 +332,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @return Incompatible track was removed bool removeTrackIfIncompatible(Vertex& vtx, std::vector& seedTracks, - FitterState_t& fitterState, + VertexFitterState& fitterState, std::vector& removedSeedTracks, const GeometryContext& geoCtx) const; @@ -354,7 +345,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// /// @return Keep new vertex bool keepNewVertex(Vertex& vtx, const std::vector& allVertices, - FitterState_t& fitterState) const; + VertexFitterState& fitterState) const; /// @brief Method that evaluates if the new vertex candidate is /// merged with one of the previously found vertices @@ -376,7 +367,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @param vertexingOptions Vertexing options Result deleteLastVertex( Vertex& vtx, std::vector>& allVertices, - std::vector& allVerticesPtr, FitterState_t& fitterState, + std::vector& allVerticesPtr, VertexFitterState& fitterState, const VertexingOptions& vertexingOptions) const; /// @brief Prepares the output vector of vertices @@ -387,9 +378,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder { /// @return The output vertex collection Result> getVertexOutputList( const std::vector& allVerticesPtr, - FitterState_t& fitterState) const; + VertexFitterState& fitterState) const; }; } // namespace Acts - -#include "Acts/Vertexing/AdaptiveMultiVertexFinder.ipp" diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp index ce803dd68e2..ecfcd10d602 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp @@ -271,5 +271,3 @@ class AdaptiveMultiVertexFitter { }; } // namespace Acts - -#include "Acts/Vertexing/AdaptiveMultiVertexFitter.ipp" diff --git a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.hpp b/Core/include/Acts/Vertexing/FullBilloirVertexFitter.hpp index e409748c497..b4989895dc2 100644 --- a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.hpp +++ b/Core/include/Acts/Vertexing/FullBilloirVertexFitter.hpp @@ -105,5 +105,3 @@ class FullBilloirVertexFitter { }; } // namespace Acts - -#include "FullBilloirVertexFitter.ipp" diff --git a/Core/include/Acts/Vertexing/GaussianTrackDensity.hpp b/Core/include/Acts/Vertexing/GaussianTrackDensity.hpp index b35f459bf5e..932219c6b50 100644 --- a/Core/include/Acts/Vertexing/GaussianTrackDensity.hpp +++ b/Core/include/Acts/Vertexing/GaussianTrackDensity.hpp @@ -22,7 +22,6 @@ namespace Acts { /// their d0 and z0 perigee parameters (mean value) and covariance /// matrices (determining the width of the function) class GaussianTrackDensity { - // @TODO: Remove template public: /// @brief Struct to store information for a single track struct TrackEntry { @@ -202,5 +201,3 @@ class GaussianTrackDensity { }; } // namespace Acts - -#include "Acts/Vertexing/GaussianTrackDensity.ipp" diff --git a/Core/include/Acts/Vertexing/GridDensityVertexFinder.ipp b/Core/include/Acts/Vertexing/GridDensityVertexFinder.ipp index 9c3585f5ba6..50af5deee3b 100644 --- a/Core/include/Acts/Vertexing/GridDensityVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/GridDensityVertexFinder.ipp @@ -29,10 +29,9 @@ auto Acts::GridDensityVertexFinder::find( } if (!couldRemoveTracks) { // No tracks were removed anymore - // Return empty seed, i.e. vertex at constraint position + // Return empty seed // (Note: Upstream finder should check for this break condition) - std::vector seedVec{vertexingOptions.constraint}; - return seedVec; + return std::vector{}; } } else { state.mainGrid = MainGridVector::Zero(); @@ -94,9 +93,7 @@ auto Acts::GridDensityVertexFinder::find( returnVertex.setFullCovariance(seedCov); - std::vector seedVec{returnVertex}; - - return seedVec; + return std::vector{returnVertex}; } template diff --git a/Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp b/Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp index 8b855bcfa4e..4f6015db549 100644 --- a/Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp +++ b/Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp @@ -98,5 +98,3 @@ class HelicalTrackLinearizer { }; } // namespace Acts - -#include "HelicalTrackLinearizer.ipp" diff --git a/Core/include/Acts/Vertexing/IterativeVertexFinder.hpp b/Core/include/Acts/Vertexing/IterativeVertexFinder.hpp index 97ac15f2ce6..f3ff4aeacb3 100644 --- a/Core/include/Acts/Vertexing/IterativeVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/IterativeVertexFinder.hpp @@ -14,7 +14,6 @@ #include "Acts/MagneticField/MagneticFieldProvider.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" -#include "Acts/Vertexing/FsmwMode1dFinder.hpp" #include "Acts/Vertexing/FullBilloirVertexFitter.hpp" #include "Acts/Vertexing/HelicalTrackLinearizer.hpp" #include "Acts/Vertexing/IVertexFinder.hpp" @@ -22,7 +21,6 @@ #include "Acts/Vertexing/TrackLinearizer.hpp" #include "Acts/Vertexing/Vertex.hpp" #include "Acts/Vertexing/VertexingOptions.hpp" -#include "Acts/Vertexing/ZScanVertexFinder.hpp" #include @@ -56,12 +54,9 @@ namespace Acts { /// from tracksAtVertex if not compatible. /// 5. Add vertex to vertexCollection /// 6. Repeat until no seedTracks are left or max. number of vertices found -/// -//////////////////////////////////////////////////////////// -/// -/// @tparam vfitter_t Vertex fitter type -template class IterativeVertexFinder final : public IVertexFinder { + using VertexFitter = FullBilloirVertexFitter; + public: /// Configuration struct struct Config { @@ -71,14 +66,14 @@ class IterativeVertexFinder final : public IVertexFinder { /// @param lin Track linearizer /// @param sfinder The seed finder /// @param est ImpactPointEstimator - Config(vfitter_t fitter, std::shared_ptr sfinder, + Config(VertexFitter fitter, std::shared_ptr sfinder, ImpactPointEstimator est) : vertexFitter(std::move(fitter)), seedFinder(std::move(sfinder)), ipEst(std::move(est)) {} /// Vertex fitter - vfitter_t vertexFitter; + VertexFitter vertexFitter; /// Track linearizer TrackLinearizer trackLinearizer; @@ -143,33 +138,7 @@ class IterativeVertexFinder final : public IVertexFinder { /// @param logger The logging instance IterativeVertexFinder(Config cfg, std::unique_ptr logger = getDefaultLogger( - "IterativeVertexFinder", Logging::INFO)) - : m_cfg(std::move(cfg)), m_logger(std::move(logger)) { - if (!m_cfg.extractParameters.connected()) { - throw std::invalid_argument( - "IterativeVertexFinder: " - "No function to extract parameters " - "provided."); - } - - if (!m_cfg.trackLinearizer.connected()) { - throw std::invalid_argument( - "IterativeVertexFinder: " - "No track linearizer provided."); - } - - if (!m_cfg.seedFinder) { - throw std::invalid_argument( - "IterativeVertexFinder: " - "No seed finder provided."); - } - - if (!m_cfg.field) { - throw std::invalid_argument( - "IterativeVertexFinder: " - "No magnetic field provider provided."); - } - } + "IterativeVertexFinder", Logging::INFO)); /// @brief Finds vertices corresponding to input trackVector /// @@ -206,11 +175,14 @@ class IterativeVertexFinder final : public IVertexFinder { /// @brief Method that calls seed finder to retrieve a vertex seed /// + /// @param state The state object /// @param seedTracks Seeding tracks /// @param vertexingOptions Vertexing options - Result getVertexSeed(State& state, - const std::vector& seedTracks, - const VertexingOptions& vertexingOptions) const; + /// + /// @return Vertex seed + Result> getVertexSeed( + State& state, const std::vector& seedTracks, + const VertexingOptions& vertexingOptions) const; /// @brief Removes all tracks in tracksToRemove from seedTracks /// @@ -288,5 +260,3 @@ class IterativeVertexFinder final : public IVertexFinder { }; } // namespace Acts - -#include "IterativeVertexFinder.ipp" diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp deleted file mode 100644 index 9ce58e47238..00000000000 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ /dev/null @@ -1,60 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/EventData/TrackParameters.hpp" -#include "Acts/Utilities/Logger.hpp" -#include "Acts/Utilities/Result.hpp" -#include "Acts/Vertexing/KalmanVertexUpdater.hpp" -#include "Acts/Vertexing/LinearizedTrack.hpp" -#include "Acts/Vertexing/TrackAtVertex.hpp" -#include "Acts/Vertexing/Vertex.hpp" - -namespace Acts { -namespace KalmanVertexTrackUpdater { - -/// KalmanVertexTrackUpdater -/// Based on -/// Ref. (1): -/// R. Frühwirth et al. -/// Vertex reconstruction and track bundling at the lep collider using robust -/// algorithms -/// Computer Physics Comm.: 96 (1996) 189 -/// Chapter 2.1 - -/// @brief Refits a single track with the knowledge of -/// the vertex it has originated from -/// -/// @param track Track to update -/// @param vtxPosFull full vertex position -/// @param vtxCovFull full vertex covariance matrix -/// @param nDimVertex number of dimensions of the vertex. Can be 3 (if we only -/// fit its spatial coordinates) or 4 (if we also fit time). -void update(TrackAtVertexRef track, const Vector4& vtxPosFull, - const SquareMatrix4& vtxCovFull, unsigned int nDimVertex); - -namespace detail { - -/// @brief Calculates a covariance matrix for the refitted track parameters -/// -/// @param wMat W_k matrix from Ref. (1) -/// @param crossCovVP Cross-covariance matrix between vertex position and track -/// momentum -/// @param vtxCov Vertex covariance matrix -/// @param newTrkParams Refitted track parameters -template -inline BoundMatrix calculateTrackCovariance( - const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, - const ActsSquareMatrix& vtxCov, - const BoundVector& newTrkParams); - -} // Namespace detail - -} // Namespace KalmanVertexTrackUpdater -} // Namespace Acts diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp index 3baecfbf7fc..77d5a8bf37e 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019-2023 CERN for the benefit of the Acts project +// Copyright (C) 2019-2024 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -8,12 +8,11 @@ #pragma once -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Utilities/Result.hpp" -#include "Acts/Vertexing/TrackAtVertex.hpp" -#include "Acts/Vertexing/Vertex.hpp" - namespace Acts { + +class Vertex; +struct TrackAtVertex; + namespace KalmanVertexUpdater { /// KalmanVertexUpdater @@ -35,91 +34,36 @@ namespace KalmanVertexUpdater { /// CERN-THESIS-2010-027 /// Section 5.3.5 -/// Cache object, the comments indicate the names of the variables in Ref. (1) -/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only -/// fit its spatial coordinates) or 4 (if we also fit time). -template -struct Cache { - using VertexVector = ActsVector; - using VertexMatrix = ActsSquareMatrix; - // \tilde{x_k} - VertexVector newVertexPos = VertexVector::Zero(); - // C_k - VertexMatrix newVertexCov = VertexMatrix::Zero(); - // C_k^-1 - VertexMatrix newVertexWeight = VertexMatrix::Zero(); - // C_{k-1}^-1 - VertexMatrix oldVertexWeight = VertexMatrix::Zero(); - // W_k - SquareMatrix3 wMat = SquareMatrix3::Zero(); -}; - /// @brief Updates vertex with knowledge of new track /// @note KalmanVertexUpdater updates the vertex when trk is added to the fit. /// However, it does not add the track to the TrackAtVertex list. This to be /// done manually after calling the method. /// -/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only -/// fit its spatial coordinates) or 4 (if we also fit time). /// /// @param vtx Vertex to be updated /// @param trk Track to be used for updating the vertex -template -void updateVertexWithTrack(Vertex& vtx, TrackAtVertex& trk); - -namespace detail { -void updateVertexWithTrack(Vector4& vtxPos, SquareMatrix4& vtxCov, - std::pair& fitQuality, - TrackAtVertexRef trk, int sign, +/// @param nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +void updateVertexWithTrack(Vertex& vtx, TrackAtVertex& trk, unsigned int nDimVertex); -// These two functions only exist so we can compile calculateUpdate in a -// compilation unit -void calculateUpdate3(const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, - Cache<3>& cache); - -void calculateUpdate4(const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, - Cache<4>& cache); -} // namespace detail +/// Based on +/// Ref. (1): +/// R. Frühwirth et al. +/// Vertex reconstruction and track bundling at the lep collider using robust +/// algorithms +/// Computer Physics Comm.: 96 (1996) 189 +/// Chapter 2.1 -/// @brief Calculates updated vertex position and covariance as well as the -/// updated track momentum when adding/removing linTrack. Saves the result in -/// cache. +/// @brief Refits a single track with the knowledge of +/// the vertex it has originated from /// -/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// @param track Track to update +/// @param vtx Vertex to use for updating the track +/// @param nDimVertex number of dimensions of the vertex. Can be 3 (if we only /// fit its spatial coordinates) or 4 (if we also fit time). -/// -/// @param vtxPos Vertex position -/// @param vtxCov Vertex covariance matrix -/// @param linTrack Linearized track to be added or removed -/// @param trackWeight Track weight -/// @param sign +1 (add track) or -1 (remove track) -/// @note Tracks are removed during the smoothing procedure to compute -/// the chi2 of the track wrt the updated vertex position -/// @param[in,out] cache A cache to store the results of this function -template -void calculateUpdate(const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, - Cache& cache) { - static_assert(nDimVertex == 3 || nDimVertex == 4, - "The vertex dimension must either be 3 (when fitting the " - "spatial coordinates) or 4 (when fitting the spatial " - "coordinates + time)."); - if constexpr (nDimVertex == 3) { - detail::calculateUpdate3(vtxPos, vtxCov, linTrack, trackWeight, sign, - cache); - } else if constexpr (nDimVertex == 4) { - detail::calculateUpdate4(vtxPos, vtxCov, linTrack, trackWeight, sign, - cache); - } -} - -} // Namespace KalmanVertexUpdater -} // Namespace Acts +void updateTrackWithVertex(TrackAtVertex& track, const Vertex& vtx, + unsigned int nDimVertex); -#include "KalmanVertexUpdater.ipp" +} // namespace KalmanVertexUpdater +} // namespace Acts diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp deleted file mode 100644 index 18854240bcb..00000000000 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ /dev/null @@ -1,18 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019-2023 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include - -template -void Acts::KalmanVertexUpdater::updateVertexWithTrack(Vertex& vtx, - TrackAtVertex& trk) { - std::pair fitQuality = vtx.fitQuality(); - detail::updateVertexWithTrack(vtx.fullPosition(), vtx.fullCovariance(), - fitQuality, trk, 1, nDimVertex); - vtx.setFitQuality(fitQuality); -} diff --git a/Core/include/Acts/Vertexing/NumericalTrackLinearizer.hpp b/Core/include/Acts/Vertexing/NumericalTrackLinearizer.hpp index da7e48a237a..1eddcd3908e 100644 --- a/Core/include/Acts/Vertexing/NumericalTrackLinearizer.hpp +++ b/Core/include/Acts/Vertexing/NumericalTrackLinearizer.hpp @@ -122,5 +122,3 @@ class NumericalTrackLinearizer { }; } // namespace Acts - -#include "NumericalTrackLinearizer.ipp" diff --git a/Core/include/Acts/Vertexing/TrackAtVertex.hpp b/Core/include/Acts/Vertexing/TrackAtVertex.hpp index 67eb9f2fa25..7072ab422d0 100644 --- a/Core/include/Acts/Vertexing/TrackAtVertex.hpp +++ b/Core/include/Acts/Vertexing/TrackAtVertex.hpp @@ -124,25 +124,6 @@ struct TrackAtVertex { bool isLinearized = false; }; -struct TrackAtVertexRef { - BoundTrackParameters& fittedParams; - double& chi2Track; - double& ndf; - double& vertexCompatibility; - double& trackWeight; - LinearizedTrack& linearizedState; - bool isLinearized; - - TrackAtVertexRef(TrackAtVertex& track) - : fittedParams(track.fittedParams), - chi2Track(track.chi2Track), - ndf(track.ndf), - vertexCompatibility(track.vertexCompatibility), - trackWeight(track.trackWeight), - linearizedState(track.linearizedState), - isLinearized(track.isLinearized) {} -}; - } // namespace Acts template <> diff --git a/Core/include/Acts/Vertexing/TrackDensityVertexFinder.hpp b/Core/include/Acts/Vertexing/TrackDensityVertexFinder.hpp index d81db1535d6..b5125015168 100644 --- a/Core/include/Acts/Vertexing/TrackDensityVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/TrackDensityVertexFinder.hpp @@ -29,15 +29,12 @@ namespace Acts { /// of the maximum of all summed track density functions. /// /// Ref. (1): https://cds.cern.ch/record/2670380 -/// -/// @tparam track_density_t The track density type -template class TrackDensityVertexFinder final : public IVertexFinder { public: /// @brief The Config struct struct Config { // The track density estimator - track_density_t trackDensityEstimator; + GaussianTrackDensity trackDensityEstimator; }; /// State struct for fulfilling interface @@ -78,5 +75,3 @@ class TrackDensityVertexFinder final : public IVertexFinder { }; } // namespace Acts - -#include "TrackDensityVertexFinder.ipp" diff --git a/Core/include/Acts/Vertexing/Vertex.hpp b/Core/include/Acts/Vertexing/Vertex.hpp index 34554edc815..da825157a0a 100644 --- a/Core/include/Acts/Vertexing/Vertex.hpp +++ b/Core/include/Acts/Vertexing/Vertex.hpp @@ -114,5 +114,3 @@ class Vertex { }; } // namespace Acts - -#include "Vertex.ipp" diff --git a/Core/include/Acts/Vertexing/Vertex.ipp b/Core/include/Acts/Vertexing/Vertex.ipp deleted file mode 100644 index 7dfe99e7ef7..00000000000 --- a/Core/include/Acts/Vertexing/Vertex.ipp +++ /dev/null @@ -1,107 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -inline Acts::Vertex::Vertex(const Vector3& position) { - m_position[ePos0] = position[ePos0]; - m_position[ePos1] = position[ePos1]; - m_position[ePos2] = position[ePos2]; -} - -inline Acts::Vertex::Vertex(const Vector4& position) : m_position(position) {} - -inline Acts::Vertex::Vertex(const Vector3& position, - const SquareMatrix3& covariance, - std::vector tracks) - : m_tracksAtVertex(std::move(tracks)) { - m_position[ePos0] = position[ePos0]; - m_position[ePos1] = position[ePos1]; - m_position[ePos2] = position[ePos2]; - m_covariance.block<3, 3>(ePos0, ePos0) = covariance; -} - -inline Acts::Vertex::Vertex(const Vector4& position, - const SquareMatrix4& covariance, - std::vector tracks) - : m_position(position), - m_covariance(covariance), - m_tracksAtVertex(std::move(tracks)) {} - -inline Acts::Vector3 Acts::Vertex::position() const { - return VectorHelpers::position(m_position); -} - -inline Acts::ActsScalar Acts::Vertex::time() const { - return m_position[eTime]; -} - -inline const Acts::Vector4& Acts::Vertex::fullPosition() const { - return m_position; -} - -inline Acts::Vector4& Acts::Vertex::fullPosition() { - return m_position; -} - -inline Acts::SquareMatrix3 Acts::Vertex::covariance() const { - return m_covariance.block<3, 3>(ePos0, ePos0); -} - -inline const Acts::SquareMatrix4& Acts::Vertex::fullCovariance() const { - return m_covariance; -} - -inline Acts::SquareMatrix4& Acts::Vertex::fullCovariance() { - return m_covariance; -} - -inline const std::vector& Acts::Vertex::tracks() const { - return m_tracksAtVertex; -} - -inline std::pair Acts::Vertex::fitQuality() const { - return std::pair(m_chiSquared, m_numberDoF); -} - -inline void Acts::Vertex::setPosition(const Vector3& position, - ActsScalar time) { - m_position[ePos0] = position[ePos0]; - m_position[ePos1] = position[ePos1]; - m_position[ePos2] = position[ePos2]; - m_position[eTime] = time; -} - -inline void Acts::Vertex::setFullPosition(const Vector4& fullPosition) { - m_position = fullPosition; -} - -inline void Acts::Vertex::setTime(ActsScalar time) { - m_position[eTime] = time; -} - -inline void Acts::Vertex::setCovariance(const SquareMatrix3& covariance) { - m_covariance.setZero(); - m_covariance.block<3, 3>(ePos0, ePos0) = covariance; -} - -inline void Acts::Vertex::setFullCovariance(const SquareMatrix4& covariance) { - m_covariance = covariance; -} - -inline void Acts::Vertex::setTracksAtVertex(std::vector tracks) { - m_tracksAtVertex = std::move(tracks); -} - -inline void Acts::Vertex::setFitQuality(double chiSquared, double numberDoF) { - m_chiSquared = chiSquared; - m_numberDoF = numberDoF; -} - -inline void Acts::Vertex::setFitQuality(std::pair fitQuality) { - m_chiSquared = fitQuality.first; - m_numberDoF = fitQuality.second; -} diff --git a/Core/include/Acts/Vertexing/ZScanVertexFinder.hpp b/Core/include/Acts/Vertexing/ZScanVertexFinder.hpp index 71417f42483..432036c1007 100644 --- a/Core/include/Acts/Vertexing/ZScanVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/ZScanVertexFinder.hpp @@ -74,14 +74,7 @@ class ZScanVertexFinder final : public IVertexFinder { /// @param logger Logging instance ZScanVertexFinder(const Config& cfg, std::unique_ptr logger = - getDefaultLogger("ZScanVertexFinder", Logging::INFO)) - : m_cfg(cfg), m_logger(std::move(logger)) { - if (!m_cfg.extractParameters.connected()) { - throw std::invalid_argument( - "ZScanVertexFinder: " - "No track parameter extractor provided."); - } - } + getDefaultLogger("ZScanVertexFinder", Logging::INFO)); /// @brief Function that determines single vertex, /// based on z0 values of input tracks, @@ -119,5 +112,3 @@ class ZScanVertexFinder final : public IVertexFinder { }; } // namespace Acts - -#include "ZScanVertexFinder.ipp" diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp b/Core/src/Vertexing/AdaptiveGridDensityVertexFinder.cpp similarity index 90% rename from Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp rename to Core/src/Vertexing/AdaptiveGridDensityVertexFinder.cpp index 14242ae8f22..f60ca2f0e09 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp +++ b/Core/src/Vertexing/AdaptiveGridDensityVertexFinder.cpp @@ -6,10 +6,13 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -inline auto Acts::AdaptiveGridDensityVertexFinder::find( +#include "Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp" + +Acts::Result> +Acts::AdaptiveGridDensityVertexFinder::find( const std::vector& trackVector, const VertexingOptions& vertexingOptions, - IVertexFinder::State& anyState) const -> Result> { + IVertexFinder::State& anyState) const { auto& state = anyState.as(); // Remove density contributions from tracks removed from track collection if (m_cfg.cacheGridStateForTrackRemoval && state.isInitialized && @@ -43,10 +46,9 @@ inline auto Acts::AdaptiveGridDensityVertexFinder::find( if (state.mainDensityMap.empty()) { // No tracks passed selection - // Return empty seed, i.e. vertex at constraint position + // Return empty seed // (Note: Upstream finder should check for this break condition) - std::vector seedVec{vertexingOptions.constraint}; - return seedVec; + return std::vector{}; } double z = 0; @@ -90,13 +92,11 @@ inline auto Acts::AdaptiveGridDensityVertexFinder::find( returnVertex.setFullCovariance(seedCov); - std::vector seedVec{returnVertex}; - - return seedVec; + return std::vector{returnVertex}; } -inline auto Acts::AdaptiveGridDensityVertexFinder::doesPassTrackSelection( - const BoundTrackParameters& trk) const -> bool { +bool Acts::AdaptiveGridDensityVertexFinder::doesPassTrackSelection( + const BoundTrackParameters& trk) const { // Get required track parameters const double d0 = trk.parameters()[BoundIndices::eBoundLoc0]; const double z0 = trk.parameters()[BoundIndices::eBoundLoc1]; diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp b/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp similarity index 84% rename from Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp rename to Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp index 5123b2f06b2..9e6ac04f36b 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp +++ b/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp @@ -1,19 +1,22 @@ // This file is part of the Acts project. // -// Copyright (C) 2020 CERN for the benefit of the Acts project +// Copyright (C) 2020-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp" + #include "Acts/Utilities/AlgebraHelpers.hpp" #include "Acts/Vertexing/VertexingError.hpp" -template -auto Acts::AdaptiveMultiVertexFinder::find( +namespace Acts { + +Acts::Result> AdaptiveMultiVertexFinder::find( const std::vector& allTracks, const VertexingOptions& vertexingOptions, - IVertexFinder::State& anyState) const -> Result> { + IVertexFinder::State& anyState) const { if (allTracks.empty()) { ACTS_ERROR("Empty track collection handed to find method"); return VertexingError::EmptyInput; @@ -27,7 +30,8 @@ auto Acts::AdaptiveMultiVertexFinder::find( // Seed tracks std::vector seedTracks = allTracks; - FitterState_t fitterState(*m_cfg.bField, vertexingOptions.magFieldContext); + VertexFitterState fitterState(*m_cfg.bField, + vertexingOptions.magFieldContext); auto seedFinderState = m_cfg.seedFinder->makeState(state.magContext); std::vector> allVertices; @@ -38,41 +42,43 @@ auto Acts::AdaptiveMultiVertexFinder::find( std::vector removedSeedTracks; while (!seedTracks.empty() && iteration < m_cfg.maxIterations && (m_cfg.addSingleTrackVertices || seedTracks.size() >= 2)) { - // Tracks that are used for searching compatible tracks - // near a vertex candidate - std::vector searchTracks; - if (m_cfg.doRealMultiVertex) { - searchTracks = origTracks; - } else { - searchTracks = seedTracks; - } Vertex currentConstraint = vertexingOptions.constraint; + // Retrieve seed vertex from all remaining seedTracks auto seedResult = doSeeding(seedTracks, currentConstraint, vertexingOptions, seedFinderState, removedSeedTracks); if (!seedResult.ok()) { return seedResult.error(); } - allVertices.push_back(std::make_unique(*seedResult)); + const auto& seedOptional = *seedResult; - Vertex& vtxCandidate = *allVertices.back(); - allVerticesPtr.push_back(&vtxCandidate); - - ACTS_DEBUG("Position of vertex candidate after seeding: " - << vtxCandidate.fullPosition().transpose()); - if (vtxCandidate.position().z() == - vertexingOptions.constraint.position().z()) { + if (!seedOptional.has_value()) { ACTS_DEBUG( "No seed found anymore. Break and stop primary vertex finding."); - allVertices.pop_back(); - allVerticesPtr.pop_back(); break; } + const auto& seedVertex = seedOptional.value(); + + ACTS_DEBUG("Position of vertex candidate after seeding: " + << seedVertex.fullPosition().transpose()); + + allVertices.push_back(std::make_unique(seedVertex)); + Vertex& vtxCandidate = *allVertices.back(); + allVerticesPtr.push_back(&vtxCandidate); // Clear the seed track collection that has been removed in last iteration // now after seed finding is done removedSeedTracks.clear(); + // Tracks that are used for searching compatible tracks near a vertex + // candidate + std::vector searchTracks; + if (m_cfg.doRealMultiVertex) { + searchTracks = origTracks; + } else { + searchTracks = seedTracks; + } + auto prepResult = canPrepareVertexForFit(searchTracks, seedTracks, vtxCandidate, currentConstraint, fitterState, vertexingOptions); @@ -140,12 +146,12 @@ auto Acts::AdaptiveMultiVertexFinder::find( return getVertexOutputList(allVerticesPtr, fitterState); } -template -auto Acts::AdaptiveMultiVertexFinder::doSeeding( +auto AdaptiveMultiVertexFinder::doSeeding( const std::vector& trackVector, Vertex& currentConstraint, const VertexingOptions& vertexingOptions, IVertexFinder::State& seedFinderState, - const std::vector& removedSeedTracks) const -> Result { + const std::vector& removedSeedTracks) const + -> Result> { VertexingOptions seedOptions = vertexingOptions; seedOptions.constraint = currentConstraint; @@ -158,8 +164,15 @@ auto Acts::AdaptiveMultiVertexFinder::doSeeding( if (!seedResult.ok()) { return seedResult.error(); } + const auto& seedVector = *seedResult; + + ACTS_DEBUG("Found " << seedVector.size() << " seeds"); + + if (seedVector.empty()) { + return std::nullopt; + } + Vertex seedVertex = seedVector.back(); - Vertex seedVertex = (*seedResult).back(); // Update constraints according to seed vertex setConstraintAfterSeeding(currentConstraint, seedOptions.useConstraintInFit, seedVertex); @@ -167,10 +180,9 @@ auto Acts::AdaptiveMultiVertexFinder::doSeeding( return seedVertex; } -template -auto Acts::AdaptiveMultiVertexFinder::setConstraintAfterSeeding( +void AdaptiveMultiVertexFinder::setConstraintAfterSeeding( Vertex& currentConstraint, bool useVertexConstraintInFit, - Vertex& seedVertex) const -> void { + Vertex& seedVertex) const { if (useVertexConstraintInFit) { if (!m_cfg.useSeedConstraint) { // Set seed vertex constraint to old constraint before seeding @@ -187,10 +199,9 @@ auto Acts::AdaptiveMultiVertexFinder::setConstraintAfterSeeding( } } -template -auto Acts::AdaptiveMultiVertexFinder::getIPSignificance( +Acts::Result AdaptiveMultiVertexFinder::getIPSignificance( const InputTrack& track, const Vertex& vtx, - const VertexingOptions& vertexingOptions) const -> Result { + const VertexingOptions& vertexingOptions) const { // TODO: In original implementation the covariance of the given vertex is set // to zero. I did the same here now, but consider removing this and just // passing the vtx object to the estimator without changing its covariance. @@ -227,11 +238,10 @@ auto Acts::AdaptiveMultiVertexFinder::getIPSignificance( return significance; } -template -auto Acts::AdaptiveMultiVertexFinder::addCompatibleTracksToVertex( +Acts::Result AdaptiveMultiVertexFinder::addCompatibleTracksToVertex( const std::vector& tracks, Vertex& vtx, - FitterState_t& fitterState, const VertexingOptions& vertexingOptions) const - -> Result { + VertexFitterState& fitterState, + const VertexingOptions& vertexingOptions) const { for (const auto& trk : tracks) { auto params = m_cfg.extractParameters(trk); auto pos = params.position(vertexingOptions.geoContext); @@ -257,13 +267,11 @@ auto Acts::AdaptiveMultiVertexFinder::addCompatibleTracksToVertex( return {}; } -template -auto Acts::AdaptiveMultiVertexFinder:: - canRecoverFromNoCompatibleTracks( - const std::vector& allTracks, - const std::vector& seedTracks, Vertex& vtx, - const Vertex& currentConstraint, FitterState_t& fitterState, - const VertexingOptions& vertexingOptions) const -> Result { +Acts::Result AdaptiveMultiVertexFinder::canRecoverFromNoCompatibleTracks( + const std::vector& allTracks, + const std::vector& seedTracks, Vertex& vtx, + const Vertex& currentConstraint, VertexFitterState& fitterState, + const VertexingOptions& vertexingOptions) const { // Recover from cases where no compatible tracks to vertex // candidate were found // TODO: This is for now how it's done in athena... this look a bit @@ -313,12 +321,11 @@ auto Acts::AdaptiveMultiVertexFinder:: return Result::success(true); } -template -auto Acts::AdaptiveMultiVertexFinder::canPrepareVertexForFit( +Acts::Result AdaptiveMultiVertexFinder::canPrepareVertexForFit( const std::vector& allTracks, const std::vector& seedTracks, Vertex& vtx, - const Vertex& currentConstraint, FitterState_t& fitterState, - const VertexingOptions& vertexingOptions) const -> Result { + const Vertex& currentConstraint, VertexFitterState& fitterState, + const VertexingOptions& vertexingOptions) const { // Add vertex info to fitter state fitterState.vtxInfoMap[&vtx] = VertexInfo(currentConstraint, vtx.fullPosition()); @@ -341,11 +348,9 @@ auto Acts::AdaptiveMultiVertexFinder::canPrepareVertexForFit( return Result::success(*resRec); } -template -auto Acts::AdaptiveMultiVertexFinder::checkVertexAndCompatibleTracks( +std::pair AdaptiveMultiVertexFinder::checkVertexAndCompatibleTracks( Vertex& vtx, const std::vector& seedTracks, - FitterState_t& fitterState, bool useVertexConstraintInFit) const - -> std::pair { + VertexFitterState& fitterState, bool useVertexConstraintInFit) const { bool isGoodVertex = false; int nCompatibleTracks = 0; for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) { @@ -380,12 +385,10 @@ auto Acts::AdaptiveMultiVertexFinder::checkVertexAndCompatibleTracks( return {nCompatibleTracks, isGoodVertex}; } -template -auto Acts::AdaptiveMultiVertexFinder:: - removeCompatibleTracksFromSeedTracks( - Vertex& vtx, std::vector& seedTracks, - FitterState_t& fitterState, - std::vector& removedSeedTracks) const -> void { +auto AdaptiveMultiVertexFinder::removeCompatibleTracksFromSeedTracks( + Vertex& vtx, std::vector& seedTracks, + VertexFitterState& fitterState, + std::vector& removedSeedTracks) const -> void { for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) { const auto& trkAtVtx = fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx)); @@ -406,11 +409,10 @@ auto Acts::AdaptiveMultiVertexFinder:: } } -template -auto Acts::AdaptiveMultiVertexFinder::removeTrackIfIncompatible( +bool AdaptiveMultiVertexFinder::removeTrackIfIncompatible( Vertex& vtx, std::vector& seedTracks, - FitterState_t& fitterState, std::vector& removedSeedTracks, - const GeometryContext& geoCtx) const -> bool { + VertexFitterState& fitterState, std::vector& removedSeedTracks, + const GeometryContext& geoCtx) const { // Try to find the track with highest compatibility double maxCompatibility = 0; @@ -462,10 +464,9 @@ auto Acts::AdaptiveMultiVertexFinder::removeTrackIfIncompatible( return true; } -template -auto Acts::AdaptiveMultiVertexFinder::keepNewVertex( +bool AdaptiveMultiVertexFinder::keepNewVertex( Vertex& vtx, const std::vector& allVertices, - FitterState_t& fitterState) const -> bool { + VertexFitterState& fitterState) const { double contamination = 0.; double contaminationNum = 0; double contaminationDeNom = 0; @@ -490,9 +491,8 @@ auto Acts::AdaptiveMultiVertexFinder::keepNewVertex( return true; } -template -auto Acts::AdaptiveMultiVertexFinder::isMergedVertex( - const Vertex& vtx, const std::vector& allVertices) const -> bool { +bool AdaptiveMultiVertexFinder::isMergedVertex( + const Vertex& vtx, const std::vector& allVertices) const { const Vector4& candidatePos = vtx.fullPosition(); const SquareMatrix4& candidateCov = vtx.fullCovariance(); @@ -544,11 +544,10 @@ auto Acts::AdaptiveMultiVertexFinder::isMergedVertex( return false; } -template -auto Acts::AdaptiveMultiVertexFinder::deleteLastVertex( +Acts::Result AdaptiveMultiVertexFinder::deleteLastVertex( Vertex& vtx, std::vector>& allVertices, - std::vector& allVerticesPtr, FitterState_t& fitterState, - const VertexingOptions& vertexingOptions) const -> Result { + std::vector& allVerticesPtr, VertexFitterState& fitterState, + const VertexingOptions& vertexingOptions) const { allVertices.pop_back(); allVerticesPtr.pop_back(); @@ -584,10 +583,10 @@ auto Acts::AdaptiveMultiVertexFinder::deleteLastVertex( return {}; } -template -auto Acts::AdaptiveMultiVertexFinder::getVertexOutputList( +Acts::Result> +AdaptiveMultiVertexFinder::getVertexOutputList( const std::vector& allVerticesPtr, - FitterState_t& fitterState) const -> Acts::Result> { + VertexFitterState& fitterState) const { std::vector outputVec; for (auto vtx : allVerticesPtr) { auto& outVtx = *vtx; @@ -601,3 +600,4 @@ auto Acts::AdaptiveMultiVertexFinder::getVertexOutputList( } return Result>(outputVec); } +} // namespace Acts diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/src/Vertexing/AdaptiveMultiVertexFitter.cpp similarity index 92% rename from Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp rename to Core/src/Vertexing/AdaptiveMultiVertexFitter.cpp index 5b6aa524115..d644a1d4be3 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/src/Vertexing/AdaptiveMultiVertexFitter.cpp @@ -6,11 +6,12 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include "Acts/Vertexing/KalmanVertexTrackUpdater.hpp" +#include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp" + #include "Acts/Vertexing/KalmanVertexUpdater.hpp" #include "Acts/Vertexing/VertexingError.hpp" -inline Acts::Result Acts::AdaptiveMultiVertexFitter::fit( +Acts::Result Acts::AdaptiveMultiVertexFitter::fit( State& state, const VertexingOptions& vertexingOptions) const { // Reset annealing tool state.annealingState = AnnealingUtility::State(); @@ -112,7 +113,7 @@ inline Acts::Result Acts::AdaptiveMultiVertexFitter::fit( return {}; } -inline Acts::Result Acts::AdaptiveMultiVertexFitter::addVtxToFit( +Acts::Result Acts::AdaptiveMultiVertexFitter::addVtxToFit( State& state, Vertex& newVertex, const VertexingOptions& vertexingOptions) const { if (state.vtxInfoMap[&newVertex].trackLinks.empty()) { @@ -185,12 +186,12 @@ inline Acts::Result Acts::AdaptiveMultiVertexFitter::addVtxToFit( return {}; } -inline bool Acts::AdaptiveMultiVertexFitter::isAlreadyInList( +bool Acts::AdaptiveMultiVertexFitter::isAlreadyInList( Vertex* vtx, const std::vector& vertices) const { return std::find(vertices.begin(), vertices.end(), vtx) != vertices.end(); } -inline Acts::Result Acts::AdaptiveMultiVertexFitter::prepareVertexForFit( +Acts::Result Acts::AdaptiveMultiVertexFitter::prepareVertexForFit( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { // Vertex info object auto& vtxInfo = state.vtxInfoMap[vtx]; @@ -211,8 +212,7 @@ inline Acts::Result Acts::AdaptiveMultiVertexFitter::prepareVertexForFit( return {}; } -inline Acts::Result -Acts::AdaptiveMultiVertexFitter::setAllVertexCompatibilities( +Acts::Result Acts::AdaptiveMultiVertexFitter::setAllVertexCompatibilities( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; @@ -252,7 +252,7 @@ Acts::AdaptiveMultiVertexFitter::setAllVertexCompatibilities( return {}; } -inline Acts::Result Acts::AdaptiveMultiVertexFitter::setWeightsAndUpdate( +Acts::Result Acts::AdaptiveMultiVertexFitter::setWeightsAndUpdate( State& state, const VertexingOptions& vertexingOptions) const { for (auto vtx : state.vertexCollection) { VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; @@ -291,11 +291,8 @@ inline Acts::Result Acts::AdaptiveMultiVertexFitter::setWeightsAndUpdate( // Update the vertex with the new track. The second template argument // corresponds to the number of fitted vertex dimensions (i.e., 3 if we // only fit spatial coordinates and 4 if we also fit time). - if (m_cfg.useTime) { - KalmanVertexUpdater::updateVertexWithTrack<4>(*vtx, trkAtVtx); - } else { - KalmanVertexUpdater::updateVertexWithTrack<3>(*vtx, trkAtVtx); - } + KalmanVertexUpdater::updateVertexWithTrack(*vtx, trkAtVtx, + m_cfg.useTime ? 4 : 3); } else { ACTS_VERBOSE("Track weight too low. Skip track."); } @@ -306,7 +303,7 @@ inline Acts::Result Acts::AdaptiveMultiVertexFitter::setWeightsAndUpdate( return {}; } -inline std::vector +std::vector Acts::AdaptiveMultiVertexFitter::collectTrackToVertexCompatibilities( State& state, const InputTrack& trk) const { // Compatibilities of trk wrt all of its associated vertices @@ -330,8 +327,7 @@ Acts::AdaptiveMultiVertexFitter::collectTrackToVertexCompatibilities( return trkToVtxCompatibilities; } -inline bool Acts::AdaptiveMultiVertexFitter::checkSmallShift( - State& state) const { +bool Acts::AdaptiveMultiVertexFitter::checkSmallShift(State& state) const { for (auto* vtx : state.vertexCollection) { Vector3 diff = state.vtxInfoMap[vtx].oldPosition.template head<3>() - vtx->position(); @@ -344,8 +340,7 @@ inline bool Acts::AdaptiveMultiVertexFitter::checkSmallShift( return true; } -inline void Acts::AdaptiveMultiVertexFitter::doVertexSmoothing( - State& state) const { +void Acts::AdaptiveMultiVertexFitter::doVertexSmoothing(State& state) const { for (const auto vtx : state.vertexCollection) { for (const auto& trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); @@ -354,15 +349,14 @@ inline void Acts::AdaptiveMultiVertexFitter::doVertexSmoothing( // vertex. The second template argument corresponds to the number of // fitted vertex dimensions (i.e., 3 if we only fit spatial coordinates // and 4 if we also fit time). - KalmanVertexTrackUpdater::update(trkAtVtx, vtx->fullPosition(), - vtx->fullCovariance(), - m_cfg.useTime ? 4 : 3); + KalmanVertexUpdater::updateTrackWithVertex(trkAtVtx, *vtx, + m_cfg.useTime ? 4 : 3); } } } } -inline void Acts::AdaptiveMultiVertexFitter::logDebugData( +void Acts::AdaptiveMultiVertexFitter::logDebugData( const State& state, const Acts::GeometryContext& geoContext) const { ACTS_DEBUG("Encountered an error when fitting the following " << state.vertexCollection.size() << " vertices:"); diff --git a/Core/src/Vertexing/CMakeLists.txt b/Core/src/Vertexing/CMakeLists.txt index ec53ddd7fb8..09ea0ae12a0 100644 --- a/Core/src/Vertexing/CMakeLists.txt +++ b/Core/src/Vertexing/CMakeLists.txt @@ -2,8 +2,18 @@ target_sources( ActsCore PRIVATE AdaptiveGridTrackDensity.cpp - KalmanVertexTrackUpdater.cpp KalmanVertexUpdater.cpp FsmwMode1dFinder.cpp VertexingError.cpp + IterativeVertexFinder.cpp + AdaptiveMultiVertexFitter.cpp + AdaptiveGridDensityVertexFinder.cpp + ZScanVertexFinder.cpp + HelicalTrackLinearizer.cpp + FullBilloirVertexFitter.cpp + AdaptiveMultiVertexFinder.cpp + Vertex.cpp + NumericalTrackLinearizer.cpp + TrackDensityVertexFinder.cpp + GaussianTrackDensity.cpp ) diff --git a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp b/Core/src/Vertexing/FullBilloirVertexFitter.cpp similarity index 96% rename from Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp rename to Core/src/Vertexing/FullBilloirVertexFitter.cpp index c640c723a61..44e9c3267c7 100644 --- a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp +++ b/Core/src/Vertexing/FullBilloirVertexFitter.cpp @@ -1,11 +1,13 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Vertexing/FullBilloirVertexFitter.hpp" + #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/MagneticField/ConstantBField.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" @@ -13,17 +15,17 @@ #include "Acts/Vertexing/TrackAtVertex.hpp" #include "Acts/Vertexing/VertexingError.hpp" -namespace Acts::detail { +namespace { /// @struct BilloirTrack /// /// @brief Struct to cache track-specific matrix operations in Billoir fitter struct BilloirTrack { - BilloirTrack(const InputTrack& params) : originalTrack(params) {} + BilloirTrack(const Acts::InputTrack& params) : originalTrack(params) {} BilloirTrack(const BilloirTrack& arg) = default; - InputTrack originalTrack; + Acts::InputTrack originalTrack; double chi2 = 0; // We drop the summation index i from Ref. (1) for better readability @@ -52,9 +54,9 @@ struct BilloirVertex { Acts::Vector4 sumBCinvU = Acts::Vector4::Zero(); }; -} // namespace Acts::detail +} // namespace -inline Acts::Result Acts::FullBilloirVertexFitter::fit( +Acts::Result Acts::FullBilloirVertexFitter::fit( const std::vector& paramVector, const VertexingOptions& vertexingOptions, MagneticFieldProvider::Cache& fieldCache) const { @@ -86,7 +88,7 @@ inline Acts::Result Acts::FullBilloirVertexFitter::fit( ndf += 3; } - std::vector billoirTracks; + std::vector billoirTracks; std::vector trackMomenta; // Initial guess of the 4D vertex position Vector4 linPoint = vertexingOptions.constraint.fullPosition(); @@ -95,7 +97,7 @@ inline Acts::Result Acts::FullBilloirVertexFitter::fit( for (int nIter = 0; nIter < m_cfg.maxIterations; ++nIter) { billoirTracks.clear(); double newChi2 = 0; - detail::BilloirVertex billoirVertex; + BilloirVertex billoirVertex; Vector3 linPointPos = VectorHelpers::position(linPoint); // Make Perigee surface at linPointPos, transverse plane of Perigee @@ -139,7 +141,7 @@ inline Acts::Result Acts::FullBilloirVertexFitter::fit( double fTheta = trackMomenta[iTrack][1]; double fQOvP = trackMomenta[iTrack][2]; double fTime = linPoint[FreeIndices::eFreeTime]; - detail::BilloirTrack billoirTrack(trackContainer); + BilloirTrack billoirTrack(trackContainer); billoirTrack.deltaQ << d0, z0, phi - fPhi, theta - fTheta, qOverP - fQOvP, t0 - fTime; diff --git a/Core/include/Acts/Vertexing/GaussianTrackDensity.ipp b/Core/src/Vertexing/GaussianTrackDensity.cpp similarity index 86% rename from Core/include/Acts/Vertexing/GaussianTrackDensity.ipp rename to Core/src/Vertexing/GaussianTrackDensity.cpp index c978e53891d..63093300e3e 100644 --- a/Core/include/Acts/Vertexing/GaussianTrackDensity.ipp +++ b/Core/src/Vertexing/GaussianTrackDensity.cpp @@ -1,19 +1,20 @@ // This file is part of the Acts project. // -// Copyright (C) 2020 CERN for the benefit of the Acts project +// Copyright (C) 2020-2024 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Vertexing/GaussianTrackDensity.hpp" + #include "Acts/Vertexing/VertexingError.hpp" #include namespace Acts { -inline std::pair -Acts::GaussianTrackDensity::globalMaximumWithWidth( +std::pair Acts::GaussianTrackDensity::globalMaximumWithWidth( State& state, const std::vector& trackList) const { auto result = addTracks(state, trackList); if (!result.ok()) { @@ -63,12 +64,12 @@ Acts::GaussianTrackDensity::globalMaximumWithWidth( std::sqrt(-(maxDensity / maxSecondDerivative))); } -inline double Acts::GaussianTrackDensity::globalMaximum( +double Acts::GaussianTrackDensity::globalMaximum( State& state, const std::vector& trackList) const { return globalMaximumWithWidth(state, trackList).first; } -inline Result Acts::GaussianTrackDensity::addTracks( +Result Acts::GaussianTrackDensity::addTracks( State& state, const std::vector& trackList) const { for (auto trk : trackList) { const BoundTrackParameters& boundParams = m_cfg.extractParameters(trk); @@ -121,7 +122,7 @@ inline Result Acts::GaussianTrackDensity::addTracks( return Result::success(); } -inline std::tuple +std::tuple Acts::GaussianTrackDensity::trackDensityAndDerivatives(State& state, double z) const { GaussianTrackDensityStore densityResult(z); @@ -131,11 +132,9 @@ Acts::GaussianTrackDensity::trackDensityAndDerivatives(State& state, return densityResult.densityAndDerivatives(); } -inline std::tuple -Acts::GaussianTrackDensity::updateMaximum(double newZ, double newValue, - double newSecondDerivative, - double maxZ, double maxValue, - double maxSecondDerivative) const { +std::tuple Acts::GaussianTrackDensity::updateMaximum( + double newZ, double newValue, double newSecondDerivative, double maxZ, + double maxValue, double maxSecondDerivative) const { if (newValue > maxValue) { maxZ = newZ; maxValue = newValue; @@ -144,13 +143,12 @@ Acts::GaussianTrackDensity::updateMaximum(double newZ, double newValue, return {maxZ, maxValue, maxSecondDerivative}; } -inline double Acts::GaussianTrackDensity::stepSize(double y, double dy, - double ddy) const { +double Acts::GaussianTrackDensity::stepSize(double y, double dy, + double ddy) const { return (m_cfg.isGaussianShaped ? (y * dy) / (dy * dy - y * ddy) : -dy / ddy); } -inline void -Acts::GaussianTrackDensity::GaussianTrackDensityStore::addTrackToDensity( +void Acts::GaussianTrackDensity::GaussianTrackDensityStore::addTrackToDensity( const TrackEntry& entry) { // Take track only if it's within bounds if (entry.lowerBound < m_z && m_z < entry.upperBound) { diff --git a/Core/include/Acts/Vertexing/HelicalTrackLinearizer.ipp b/Core/src/Vertexing/HelicalTrackLinearizer.cpp similarity index 98% rename from Core/include/Acts/Vertexing/HelicalTrackLinearizer.ipp rename to Core/src/Vertexing/HelicalTrackLinearizer.cpp index f15e41356e8..e3d4fd27606 100644 --- a/Core/include/Acts/Vertexing/HelicalTrackLinearizer.ipp +++ b/Core/src/Vertexing/HelicalTrackLinearizer.cpp @@ -6,10 +6,12 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Vertexing/HelicalTrackLinearizer.hpp" + #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Vertexing/LinearizerTrackParameters.hpp" -inline Acts::Result +Acts::Result Acts::HelicalTrackLinearizer::linearizeTrack( const BoundTrackParameters& params, double linPointTime, const Surface& perigeeSurface, const Acts::GeometryContext& gctx, diff --git a/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp b/Core/src/Vertexing/IterativeVertexFinder.cpp similarity index 90% rename from Core/include/Acts/Vertexing/IterativeVertexFinder.ipp rename to Core/src/Vertexing/IterativeVertexFinder.cpp index be7ec15d3ce..300b58fb9dd 100644 --- a/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp +++ b/Core/src/Vertexing/IterativeVertexFinder.cpp @@ -1,13 +1,46 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -template -auto Acts::IterativeVertexFinder::find( +#include "Acts/Vertexing/IterativeVertexFinder.hpp" + +#include "Acts/Surfaces/PerigeeSurface.hpp" +#include "Acts/Vertexing/VertexingError.hpp" + +Acts::IterativeVertexFinder::IterativeVertexFinder( + Config cfg, std::unique_ptr logger) + : m_cfg(std::move(cfg)), m_logger(std::move(logger)) { + if (!m_cfg.extractParameters.connected()) { + throw std::invalid_argument( + "IterativeVertexFinder: " + "No function to extract parameters " + "provided."); + } + + if (!m_cfg.trackLinearizer.connected()) { + throw std::invalid_argument( + "IterativeVertexFinder: " + "No track linearizer provided."); + } + + if (!m_cfg.seedFinder) { + throw std::invalid_argument( + "IterativeVertexFinder: " + "No seed finder provided."); + } + + if (!m_cfg.field) { + throw std::invalid_argument( + "IterativeVertexFinder: " + "No magnetic field provider provided."); + } +} + +auto Acts::IterativeVertexFinder::find( const std::vector& trackVector, const VertexingOptions& vertexingOptions, IVertexFinder::State& anyState) const -> Result> { @@ -29,14 +62,13 @@ auto Acts::IterativeVertexFinder::find( if (!seedRes.ok()) { return seedRes.error(); } + const auto& seedOptional = *seedRes; - const auto& seedVertex = *seedRes; - - if (seedVertex.fullPosition()[eZ] == - vertexingOptions.constraint.position().z()) { + if (!seedOptional.has_value()) { ACTS_DEBUG("No more seed found. Break and stop primary vertex finding."); break; } + const auto& seedVertex = *seedOptional; /// End seeding /// Now take only tracks compatible with current seed @@ -155,10 +187,10 @@ auto Acts::IterativeVertexFinder::find( return vertexCollection; } -template -auto Acts::IterativeVertexFinder::getVertexSeed( +auto Acts::IterativeVertexFinder::getVertexSeed( State& state, const std::vector& seedTracks, - const VertexingOptions& vertexingOptions) const -> Result { + const VertexingOptions& vertexingOptions) const + -> Result> { auto finderState = m_cfg.seedFinder->makeState(state.magContext); auto res = m_cfg.seedFinder->find(seedTracks, vertexingOptions, finderState); @@ -167,20 +199,14 @@ auto Acts::IterativeVertexFinder::getVertexSeed( << seedTracks.size()); return VertexingError::SeedingError; } + const auto& seedVector = *res; - const auto& vertexCollection = *res; + ACTS_DEBUG("Found " << seedVector.size() << " seeds"); - if (vertexCollection.empty()) { - ACTS_ERROR("Empty seed collection was returned. Number of input tracks: " - << seedTracks.size()); - return VertexingError::SeedingError; + if (seedVector.empty()) { + return std::nullopt; } - - ACTS_DEBUG("Found " << vertexCollection.size() << " seeds"); - - // retrieve the seed vertex as the last element in - // the seed vertexCollection - Vertex seedVertex = vertexCollection.back(); + const Vertex& seedVertex = seedVector.back(); ACTS_DEBUG("Use " << seedTracks.size() << " tracks for vertex seed finding.") ACTS_DEBUG( @@ -189,8 +215,7 @@ auto Acts::IterativeVertexFinder::getVertexSeed( return seedVertex; } -template -void Acts::IterativeVertexFinder::removeTracks( +inline void Acts::IterativeVertexFinder::removeTracks( const std::vector& tracksToRemove, std::vector& seedTracks) const { for (const auto& trk : tracksToRemove) { @@ -210,8 +235,7 @@ void Acts::IterativeVertexFinder::removeTracks( } } -template -Acts::Result Acts::IterativeVertexFinder::getCompatibility( +Acts::Result Acts::IterativeVertexFinder::getCompatibility( const BoundTrackParameters& params, const Vertex& vertex, const Surface& perigeeSurface, const VertexingOptions& vertexingOptions, State& state) const { @@ -246,9 +270,7 @@ Acts::Result Acts::IterativeVertexFinder::getCompatibility( return compatibility; } -template -Acts::Result -Acts::IterativeVertexFinder::removeUsedCompatibleTracks( +Acts::Result Acts::IterativeVertexFinder::removeUsedCompatibleTracks( Vertex& vertex, std::vector& tracksToFit, std::vector& seedTracks, const VertexingOptions& vertexingOptions, State& state) const { @@ -339,8 +361,7 @@ Acts::IterativeVertexFinder::removeUsedCompatibleTracks( return {}; } -template -Acts::Result Acts::IterativeVertexFinder::fillTracksToFit( +Acts::Result Acts::IterativeVertexFinder::fillTracksToFit( const std::vector& seedTracks, const Vertex& seedVertex, std::vector& tracksToFitOut, std::vector& tracksToFitSplitVertexOut, @@ -412,9 +433,7 @@ Acts::Result Acts::IterativeVertexFinder::fillTracksToFit( return {}; } -template -Acts::Result -Acts::IterativeVertexFinder::reassignTracksToNewVertex( +Acts::Result Acts::IterativeVertexFinder::reassignTracksToNewVertex( std::vector& vertexCollection, Vertex& currentVertex, std::vector& tracksToFit, std::vector& seedTracks, const std::vector& /* origTracks */, @@ -546,8 +565,7 @@ Acts::IterativeVertexFinder::reassignTracksToNewVertex( return Result::success(isGoodVertex); } -template -int Acts::IterativeVertexFinder::countSignificantTracks( +int Acts::IterativeVertexFinder::countSignificantTracks( const Vertex& vtx) const { return std::count_if(vtx.tracks().begin(), vtx.tracks().end(), [this](const TrackAtVertex& trk) { diff --git a/Core/src/Vertexing/KalmanVertexTrackUpdater.cpp b/Core/src/Vertexing/KalmanVertexTrackUpdater.cpp deleted file mode 100644 index 243335fc5ad..00000000000 --- a/Core/src/Vertexing/KalmanVertexTrackUpdater.cpp +++ /dev/null @@ -1,209 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2023 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "Acts/Vertexing/KalmanVertexTrackUpdater.hpp" - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Surfaces/PerigeeSurface.hpp" -#include "Acts/Vertexing/KalmanVertexUpdater.hpp" -#include "Acts/Vertexing/TrackAtVertex.hpp" - -namespace Acts { - -namespace { - -template -Acts::BoundMatrix calculateTrackCovariance( - const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, - const ActsSquareMatrix& vtxCov, - const BoundVector& newTrkParams) { - // D_k^n - ActsSquareMatrix<3> momCov = - wMat + crossCovVP.transpose() * vtxCov.inverse() * crossCovVP; - - // Full x, y, z, phi, theta, q/p, and, optionally, t covariance matrix. Note - // that we call this set of parameters "free" in the following even though - // that is not quite correct (free parameters correspond to - // x, y, z, t, px, py, pz) - constexpr unsigned int nFreeParams = nDimVertex + 3; - ActsSquareMatrix freeTrkCov( - ActsSquareMatrix::Zero()); - - freeTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); - freeTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); - freeTrkCov.template block<3, 3>(3, 0) = - crossCovVP.template block<3, 3>(0, 0).transpose(); - freeTrkCov.template block<3, 3>(3, 3) = momCov; - - // Fill time (cross-)covariances - if constexpr (nFreeParams == 7) { - freeTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); - freeTrkCov.template block<3, 1>(3, 6) = - crossCovVP.template block<1, 3>(3, 0).transpose(); - freeTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); - freeTrkCov.template block<1, 3>(6, 3) = - crossCovVP.template block<1, 3>(3, 0); - freeTrkCov(6, 6) = vtxCov(3, 3); - } - - // Jacobian relating "free" and bound track parameters - constexpr unsigned int nBoundParams = nDimVertex + 2; - ActsMatrix freeToBoundJac( - ActsMatrix::Zero()); - - // TODO: Jacobian is not quite correct - // First row - freeToBoundJac(0, 0) = -std::sin(newTrkParams[2]); - freeToBoundJac(0, 1) = std::cos(newTrkParams[2]); - - double tanTheta = std::tan(newTrkParams[3]); - - // Second row - freeToBoundJac(1, 0) = -freeToBoundJac(0, 1) / tanTheta; - freeToBoundJac(1, 1) = freeToBoundJac(0, 0) / tanTheta; - - // Dimension of the part of the Jacobian that is an identity matrix - constexpr unsigned int nDimIdentity = nFreeParams - 2; - freeToBoundJac.template block(1, 2) = - ActsMatrix::Identity(); - - // Full perigee track covariance - BoundMatrix boundTrackCov(BoundMatrix::Identity()); - boundTrackCov.block(0, 0) = - (freeToBoundJac * (freeTrkCov * freeToBoundJac.transpose())); - - return boundTrackCov; -} - -template -void updateImpl(TrackAtVertexRef track, const Vector4& vtxPosFull, - const SquareMatrix4& vtxCovFull) { - if constexpr (nDimVertex != 3 && nDimVertex != 4) { - throw std::invalid_argument( - "The vertex dimension must either be 3 (when fitting the spatial " - "coordinates) or 4 (when fitting the spatial coordinates + time)."); - } - - using VertexVector = ActsVector; - using VertexMatrix = ActsSquareMatrix; - constexpr unsigned int nBoundParams = nDimVertex + 2; - using ParameterVector = ActsVector; - using ParameterMatrix = ActsSquareMatrix; - // Check if linearized state exists - if (!track.isLinearized) { - throw std::invalid_argument("TrackAtVertex object must be linearized."); - } - - // Extract vertex position and covariance - // \tilde{x_n} - const VertexVector vtxPos = vtxPosFull.template head(); - // C_n - const VertexMatrix vtxCov = - vtxCovFull.template block(0, 0); - - // Get the linearized track - const LinearizedTrack& linTrack = track.linearizedState; - // Retrieve variables from the track linearization. The comments indicate the - // corresponding symbol used in Ref. (1). - // A_k - const ActsMatrix posJac = - linTrack.positionJacobian.block(0, 0); - // B_k - const ActsMatrix momJac = - linTrack.momentumJacobian.block(0, 0); - // p_k - const ParameterVector trkParams = - linTrack.parametersAtPCA.head(); - // c_k - const ParameterVector constTerm = linTrack.constantTerm.head(); - // TODO we could use `linTrack.weightAtPCA` but only if we would always fit - // time. - // G_k - const ParameterMatrix trkParamWeight = - linTrack.covarianceAtPCA.block(0, 0) - .inverse(); - - // Cache object filled with zeros - KalmanVertexUpdater::Cache cache; - - // Calculate the update of the vertex position when the track is removed. This - // might be unintuitive, but it is needed to compute a symmetric chi2. - KalmanVertexUpdater::calculateUpdate(vtxPosFull, vtxCovFull, linTrack, - track.trackWeight, -1, cache); - - // Refit track momentum with the final vertex position - Vector3 newTrkMomentum = cache.wMat * momJac.transpose() * trkParamWeight * - (trkParams - constTerm - posJac * vtxPos); - - // Track parameters, impact parameters are set to 0 and momentum corresponds - // to newTrkMomentum. TODO: Make transition fitterParams -> fittedMomentum. - BoundVector newTrkParams(BoundVector::Zero()); - - // Get phi and theta and correct for possible periodicity changes - const auto correctedPhiTheta = - Acts::detail::normalizePhiTheta(newTrkMomentum(0), newTrkMomentum(1)); - newTrkParams(BoundIndices::eBoundPhi) = correctedPhiTheta.first; // phi - newTrkParams(BoundIndices::eBoundTheta) = correctedPhiTheta.second; // theta - newTrkParams(BoundIndices::eBoundQOverP) = newTrkMomentum(2); // qOverP - - // E_k^n - const ActsMatrix crossCovVP = - -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; - - // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). - VertexVector posDiff = - vtxPos - cache.newVertexPos.template head(); - - // r_k^n - ParameterVector paramDiff = - trkParams - (constTerm + posJac * vtxPos + momJac * newTrkMomentum); - - // New chi2 to be set later - double chi2 = - posDiff.dot( - cache.newVertexWeight.template block(0, 0) * - posDiff) + - paramDiff.dot(trkParamWeight * paramDiff); - - Acts::BoundMatrix newTrackCov = calculateTrackCovariance( - cache.wMat, crossCovVP, vtxCov, newTrkParams); - - // Create new refitted parameters - std::shared_ptr perigeeSurface = - Surface::makeShared(vtxPos.template head<3>()); - - BoundTrackParameters refittedPerigee = - BoundTrackParameters(perigeeSurface, newTrkParams, std::move(newTrackCov), - track.fittedParams.particleHypothesis()); - - // Set new properties - track.fittedParams = refittedPerigee; - track.chi2Track = chi2; - track.ndf = 2 * track.trackWeight; - - return; -} - -} // namespace - -void Acts::KalmanVertexTrackUpdater::update(TrackAtVertexRef track, - const Vector4& vtxPosFull, - const SquareMatrix4& vtxCovFull, - unsigned int nDimVertex) { - if (nDimVertex == 3) { - updateImpl<3>(track, vtxPosFull, vtxCovFull); - } else if (nDimVertex == 4) { - updateImpl<4>(track, vtxPosFull, vtxCovFull); - } else { - throw std::invalid_argument( - "The vertex dimension must either be 3 (when fitting the spatial " - "coordinates) or 4 (when fitting the spatial coordinates + time)."); - } -} -} // namespace Acts diff --git a/Core/src/Vertexing/KalmanVertexUpdater.cpp b/Core/src/Vertexing/KalmanVertexUpdater.cpp index ac4f0769f4b..05e1471e5d2 100644 --- a/Core/src/Vertexing/KalmanVertexUpdater.cpp +++ b/Core/src/Vertexing/KalmanVertexUpdater.cpp @@ -9,15 +9,50 @@ #include "Acts/Vertexing/KalmanVertexUpdater.hpp" #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Vertexing/LinearizedTrack.hpp" +#include "Acts/Vertexing/Vertex.hpp" namespace Acts { namespace { + +/// Cache object, the comments indicate the names of the variables in Ref. (1) +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +template +struct Cache { + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; + // \tilde{x_k} + VertexVector newVertexPos = VertexVector::Zero(); + // C_k + VertexMatrix newVertexCov = VertexMatrix::Zero(); + // C_k^-1 + VertexMatrix newVertexWeight = VertexMatrix::Zero(); + // C_{k-1}^-1 + VertexMatrix oldVertexWeight = VertexMatrix::Zero(); + // W_k + SquareMatrix3 wMat = SquareMatrix3::Zero(); +}; + +/// @brief Calculates updated vertex position and covariance as well as the +/// updated track momentum when adding/removing linTrack. Saves the result in +/// cache. +/// +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// +/// @param vtx Vertex +/// @param linTrack Linearized track to be added or removed +/// @param trackWeight Track weight +/// @param sign +1 (add track) or -1 (remove track) +/// @note Tracks are removed during the smoothing procedure to compute +/// the chi2 of the track wrt the updated vertex position +/// @param[out] cache A cache to store the results of this function template -void calculateUpdateImpl(const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, - KalmanVertexUpdater::Cache& cache) { +void calculateUpdate(const Vertex& vtx, const Acts::LinearizedTrack& linTrack, + const double trackWeight, const int sign, + Cache& cache) { constexpr unsigned int nBoundParams = nDimVertex + 2; using ParameterVector = ActsVector; using ParameterMatrix = ActsSquareMatrix; @@ -48,10 +83,12 @@ void calculateUpdateImpl(const Vector4& vtxPos, const SquareMatrix4& vtxCov, .inverse(); // Retrieve current position of the vertex and its current weight matrix - const ActsVector oldVtxPos = vtxPos.template head(); + const ActsVector oldVtxPos = + vtx.fullPosition().template head(); // C_{k-1}^-1 cache.oldVertexWeight = - (vtxCov.template block(0, 0)).inverse(); + (vtx.fullCovariance().template block(0, 0)) + .inverse(); // W_k cache.wMat = (momJac.transpose() * (trkParamWeight * momJac)).inverse(); @@ -76,9 +113,8 @@ void calculateUpdateImpl(const Vector4& vtxPos, const SquareMatrix4& vtxCov, } template -double vertexPositionChi2Update( - const Vector4& oldVtxPos, - const KalmanVertexUpdater::Cache& cache) { +double vertexPositionChi2Update(const Vector4& oldVtxPos, + const Cache& cache) { ActsVector posDiff = cache.newVertexPos - oldVtxPos.template head(); @@ -87,9 +123,8 @@ double vertexPositionChi2Update( } template -double trackParametersChi2( - const LinearizedTrack& linTrack, - const KalmanVertexUpdater::Cache& cache) { +double trackParametersChi2(const LinearizedTrack& linTrack, + const Cache& cache) { constexpr unsigned int nBoundParams = nDimVertex + 2; using ParameterVector = ActsVector; using ParameterMatrix = ActsSquareMatrix; @@ -130,9 +165,7 @@ double trackParametersChi2( } template -void update(Vector4& vtxPos, SquareMatrix4& vtxCov, - std::pair& fitQuality, TrackAtVertexRef trk, - int sign) { +void updateVertexWithTrackImpl(Vertex& vtx, TrackAtVertex& trk, int sign) { if constexpr (nDimVertex != 3 && nDimVertex != 4) { throw std::invalid_argument( "The vertex dimension must either be 3 (when fitting the spatial " @@ -142,22 +175,21 @@ void update(Vector4& vtxPos, SquareMatrix4& vtxCov, double trackWeight = trk.trackWeight; // Set up cache where entire content is set to 0 - KalmanVertexUpdater::Cache cache; + Cache cache; // Calculate update and save result in cache - calculateUpdate(vtxPos, vtxCov, trk.linearizedState, trackWeight, sign, - cache); + calculateUpdate(vtx, trk.linearizedState, trackWeight, sign, cache); // Get fit quality parameters wrt to old vertex double chi2 = 0.; double ndf = 0.; - std::tie(chi2, ndf) = fitQuality; + std::tie(chi2, ndf) = vtx.fitQuality(); // Chi2 of the track parameters double trkChi2 = trackParametersChi2(trk.linearizedState, cache); // Update of the chi2 of the vertex position - double vtxPosChi2Update = vertexPositionChi2Update(vtxPos, cache); + double vtxPosChi2Update = vertexPositionChi2Update(vtx.fullPosition(), cache); // Calculate new chi2 chi2 += sign * (vtxPosChi2Update + trackWeight * trkChi2); @@ -167,14 +199,14 @@ void update(Vector4& vtxPos, SquareMatrix4& vtxCov, // Updating the vertex if constexpr (nDimVertex == 3) { - vtxPos.head<3>() = cache.newVertexPos.template head<3>(); - vtxCov.template topLeftCorner<3, 3>() = + vtx.fullPosition().head<3>() = cache.newVertexPos.template head<3>(); + vtx.fullCovariance().template topLeftCorner<3, 3>() = cache.newVertexCov.template topLeftCorner<3, 3>(); } else if constexpr (nDimVertex == 4) { - vtxPos = cache.newVertexPos; - vtxCov = cache.newVertexCov; + vtx.fullPosition() = cache.newVertexPos; + vtx.fullCovariance() = cache.newVertexCov; } - fitQuality = {chi2, ndf}; + vtx.setFitQuality(chi2, ndf); if (sign == 1) { // Update track @@ -190,30 +222,212 @@ void update(Vector4& vtxPos, SquareMatrix4& vtxCov, } } -} // namespace +/// @brief Calculates a covariance matrix for the refitted track parameters +/// +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// +/// @param wMat W_k matrix from Ref. (1) +/// @param crossCovVP Cross-covariance matrix between vertex position and track +/// momentum +/// @param vtxCov Vertex covariance matrix +/// @param newTrkParams Refitted track parameters +template +Acts::BoundMatrix calculateTrackCovariance( + const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, + const ActsSquareMatrix& vtxCov, + const BoundVector& newTrkParams) { + // D_k^n + ActsSquareMatrix<3> momCov = + wMat + crossCovVP.transpose() * vtxCov.inverse() * crossCovVP; + + // Full x, y, z, phi, theta, q/p, and, optionally, t covariance matrix. Note + // that we call this set of parameters "free" in the following even though + // that is not quite correct (free parameters correspond to + // x, y, z, t, px, py, pz) + constexpr unsigned int nFreeParams = nDimVertex + 3; + ActsSquareMatrix freeTrkCov( + ActsSquareMatrix::Zero()); + + freeTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(3, 0) = + crossCovVP.template block<3, 3>(0, 0).transpose(); + freeTrkCov.template block<3, 3>(3, 3) = momCov; + + // Fill time (cross-)covariances + if constexpr (nFreeParams == 7) { + freeTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); + freeTrkCov.template block<3, 1>(3, 6) = + crossCovVP.template block<1, 3>(3, 0).transpose(); + freeTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); + freeTrkCov.template block<1, 3>(6, 3) = + crossCovVP.template block<1, 3>(3, 0); + freeTrkCov(6, 6) = vtxCov(3, 3); + } + + // Jacobian relating "free" and bound track parameters + constexpr unsigned int nBoundParams = nDimVertex + 2; + ActsMatrix freeToBoundJac( + ActsMatrix::Zero()); + + // TODO: Jacobian is not quite correct + // First row + freeToBoundJac(0, 0) = -std::sin(newTrkParams[2]); + freeToBoundJac(0, 1) = std::cos(newTrkParams[2]); + + double tanTheta = std::tan(newTrkParams[3]); + + // Second row + freeToBoundJac(1, 0) = -freeToBoundJac(0, 1) / tanTheta; + freeToBoundJac(1, 1) = freeToBoundJac(0, 0) / tanTheta; + + // Dimension of the part of the Jacobian that is an identity matrix + constexpr unsigned int nDimIdentity = nFreeParams - 2; + freeToBoundJac.template block(1, 2) = + ActsMatrix::Identity(); + + // Full perigee track covariance + BoundMatrix boundTrackCov(BoundMatrix::Identity()); + boundTrackCov.block(0, 0) = + (freeToBoundJac * (freeTrkCov * freeToBoundJac.transpose())); + + return boundTrackCov; +} + +template +void updateTrackWithVertexImpl(TrackAtVertex& track, const Vertex& vtx) { + if constexpr (nDimVertex != 3 && nDimVertex != 4) { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } + + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; + // Check if linearized state exists + if (!track.isLinearized) { + throw std::invalid_argument("TrackAtVertex object must be linearized."); + } + + // Extract vertex position and covariance + // \tilde{x_n} + const VertexVector vtxPos = vtx.fullPosition().template head(); + // C_n + const VertexMatrix vtxCov = + vtx.fullCovariance().template block(0, 0); + + // Get the linearized track + const LinearizedTrack& linTrack = track.linearizedState; + // Retrieve variables from the track linearization. The comments indicate the + // corresponding symbol used in Ref. (1). + // A_k + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); + // B_k + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); + // p_k + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); + // c_k + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time. + // G_k + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); + + // Cache object filled with zeros + Cache cache; + + // Calculate the update of the vertex position when the track is removed. This + // might be unintuitive, but it is needed to compute a symmetric chi2. + calculateUpdate(vtx, linTrack, track.trackWeight, -1, cache); + + // Refit track momentum with the final vertex position + Vector3 newTrkMomentum = cache.wMat * momJac.transpose() * trkParamWeight * + (trkParams - constTerm - posJac * vtxPos); + + // Track parameters, impact parameters are set to 0 and momentum corresponds + // to newTrkMomentum. TODO: Make transition fitterParams -> fittedMomentum. + BoundVector newTrkParams(BoundVector::Zero()); -void KalmanVertexUpdater::detail::calculateUpdate3( - const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, const double trackWeight, - const int sign, Cache<3>& cache) { - calculateUpdateImpl<3>(vtxPos, vtxCov, linTrack, trackWeight, sign, cache); + // Get phi and theta and correct for possible periodicity changes + const auto correctedPhiTheta = + Acts::detail::normalizePhiTheta(newTrkMomentum(0), newTrkMomentum(1)); + newTrkParams(BoundIndices::eBoundPhi) = correctedPhiTheta.first; // phi + newTrkParams(BoundIndices::eBoundTheta) = correctedPhiTheta.second; // theta + newTrkParams(BoundIndices::eBoundQOverP) = newTrkMomentum(2); // qOverP + + // E_k^n + const ActsMatrix crossCovVP = + -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; + + // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). + VertexVector posDiff = + vtxPos - cache.newVertexPos.template head(); + + // r_k^n + ParameterVector paramDiff = + trkParams - (constTerm + posJac * vtxPos + momJac * newTrkMomentum); + + // New chi2 to be set later + double chi2 = + posDiff.dot( + cache.newVertexWeight.template block(0, 0) * + posDiff) + + paramDiff.dot(trkParamWeight * paramDiff); + + Acts::BoundMatrix newTrackCov = calculateTrackCovariance( + cache.wMat, crossCovVP, vtxCov, newTrkParams); + + // Create new refitted parameters + std::shared_ptr perigeeSurface = + Surface::makeShared(vtxPos.template head<3>()); + + BoundTrackParameters refittedPerigee = + BoundTrackParameters(perigeeSurface, newTrkParams, std::move(newTrackCov), + track.fittedParams.particleHypothesis()); + + // Set new properties + track.fittedParams = refittedPerigee; + track.chi2Track = chi2; + track.ndf = 2 * track.trackWeight; + + return; } -void KalmanVertexUpdater::detail::calculateUpdate4( - const Vector4& vtxPos, const SquareMatrix4& vtxCov, - const Acts::LinearizedTrack& linTrack, const double trackWeight, - const int sign, Cache<4>& cache) { - calculateUpdateImpl<4>(vtxPos, vtxCov, linTrack, trackWeight, sign, cache); +} // namespace + +// The two functions don't contain any of the actual update code, they +// only dispatch into templated functions, effectively doing a +// runtime-to-compile time conversion. + +void KalmanVertexUpdater::updateVertexWithTrack(Vertex& vtx, TrackAtVertex& trk, + unsigned int nDimVertex) { + if (nDimVertex == 3) { + updateVertexWithTrackImpl<3>(vtx, trk, 1); + } else if (nDimVertex == 4) { + updateVertexWithTrackImpl<4>(vtx, trk, 1); + } else { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } } -void KalmanVertexUpdater::detail::updateVertexWithTrack( - Vector4& vtxPos, SquareMatrix4& vtxCov, - std::pair& fitQuality, TrackAtVertexRef trk, int sign, - unsigned int nDimVertex) { +void Acts::KalmanVertexUpdater::updateTrackWithVertex(TrackAtVertex& track, + const Vertex& vtx, + unsigned int nDimVertex) { if (nDimVertex == 3) { - update<3>(vtxPos, vtxCov, fitQuality, trk, sign); + updateTrackWithVertexImpl<3>(track, vtx); } else if (nDimVertex == 4) { - update<4>(vtxPos, vtxCov, fitQuality, trk, sign); + updateTrackWithVertexImpl<4>(track, vtx); } else { throw std::invalid_argument( "The vertex dimension must either be 3 (when fitting the spatial " diff --git a/Core/include/Acts/Vertexing/NumericalTrackLinearizer.ipp b/Core/src/Vertexing/NumericalTrackLinearizer.cpp similarity index 97% rename from Core/include/Acts/Vertexing/NumericalTrackLinearizer.ipp rename to Core/src/Vertexing/NumericalTrackLinearizer.cpp index ce40f390827..a6ec9a7d1b8 100644 --- a/Core/include/Acts/Vertexing/NumericalTrackLinearizer.ipp +++ b/Core/src/Vertexing/NumericalTrackLinearizer.cpp @@ -1,16 +1,18 @@ // This file is part of the Acts project. // -// Copyright (C) 2023 CERN for the benefit of the Acts project +// Copyright (C) 2023-2024 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include "Acts/Vertexing/NumericalTrackLinearizer.hpp" + #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Utilities/UnitVectors.hpp" #include "Acts/Vertexing/LinearizerTrackParameters.hpp" -inline Acts::Result +Acts::Result Acts::NumericalTrackLinearizer::linearizeTrack( const BoundTrackParameters& params, double linPointTime, const Surface& perigeeSurface, const Acts::GeometryContext& gctx, diff --git a/Core/include/Acts/Vertexing/TrackDensityVertexFinder.ipp b/Core/src/Vertexing/TrackDensityVertexFinder.cpp similarity index 76% rename from Core/include/Acts/Vertexing/TrackDensityVertexFinder.ipp rename to Core/src/Vertexing/TrackDensityVertexFinder.cpp index 1d1014014ef..fda66cec8cc 100644 --- a/Core/include/Acts/Vertexing/TrackDensityVertexFinder.ipp +++ b/Core/src/Vertexing/TrackDensityVertexFinder.cpp @@ -1,17 +1,18 @@ // This file is part of the Acts project. // -// Copyright (C) 2020 CERN for the benefit of the Acts project +// Copyright (C) 2020-2024 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -template -auto Acts::TrackDensityVertexFinder::find( +#include "Acts/Vertexing/TrackDensityVertexFinder.hpp" + +Acts::Result> Acts::TrackDensityVertexFinder::find( const std::vector& trackVector, const VertexingOptions& vertexingOptions, - IVertexFinder::State& /*state*/) const -> Result> { - typename track_density_t::State densityState(trackVector.size()); + IVertexFinder::State& /*state*/) const { + GaussianTrackDensity::State densityState(trackVector.size()); // Calculate z seed position std::pair zAndWidth = @@ -36,7 +37,5 @@ auto Acts::TrackDensityVertexFinder::find( returnVertex.setFullCovariance(seedCov); - std::vector seedVec{returnVertex}; - - return seedVec; + return std::vector{returnVertex}; } diff --git a/Core/src/Vertexing/Vertex.cpp b/Core/src/Vertexing/Vertex.cpp new file mode 100644 index 00000000000..0d117e4e7fd --- /dev/null +++ b/Core/src/Vertexing/Vertex.cpp @@ -0,0 +1,106 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2019-2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Vertexing/Vertex.hpp" + +Acts::Vertex::Vertex(const Vector3& position) { + m_position[ePos0] = position[ePos0]; + m_position[ePos1] = position[ePos1]; + m_position[ePos2] = position[ePos2]; +} + +Acts::Vertex::Vertex(const Vector4& position) : m_position(position) {} + +Acts::Vertex::Vertex(const Vector3& position, const SquareMatrix3& covariance, + std::vector tracks) + : m_tracksAtVertex(std::move(tracks)) { + m_position[ePos0] = position[ePos0]; + m_position[ePos1] = position[ePos1]; + m_position[ePos2] = position[ePos2]; + m_covariance.block<3, 3>(ePos0, ePos0) = covariance; +} + +Acts::Vertex::Vertex(const Vector4& position, const SquareMatrix4& covariance, + std::vector tracks) + : m_position(position), + m_covariance(covariance), + m_tracksAtVertex(std::move(tracks)) {} + +Acts::Vector3 Acts::Vertex::position() const { + return VectorHelpers::position(m_position); +} + +Acts::ActsScalar Acts::Vertex::time() const { + return m_position[eTime]; +} + +const Acts::Vector4& Acts::Vertex::fullPosition() const { + return m_position; +} + +Acts::Vector4& Acts::Vertex::fullPosition() { + return m_position; +} + +Acts::SquareMatrix3 Acts::Vertex::covariance() const { + return m_covariance.block<3, 3>(ePos0, ePos0); +} + +const Acts::SquareMatrix4& Acts::Vertex::fullCovariance() const { + return m_covariance; +} + +Acts::SquareMatrix4& Acts::Vertex::fullCovariance() { + return m_covariance; +} + +const std::vector& Acts::Vertex::tracks() const { + return m_tracksAtVertex; +} + +std::pair Acts::Vertex::fitQuality() const { + return std::pair(m_chiSquared, m_numberDoF); +} + +void Acts::Vertex::setPosition(const Vector3& position, ActsScalar time) { + m_position[ePos0] = position[ePos0]; + m_position[ePos1] = position[ePos1]; + m_position[ePos2] = position[ePos2]; + m_position[eTime] = time; +} + +void Acts::Vertex::setFullPosition(const Vector4& fullPosition) { + m_position = fullPosition; +} + +void Acts::Vertex::setTime(ActsScalar time) { + m_position[eTime] = time; +} + +void Acts::Vertex::setCovariance(const SquareMatrix3& covariance) { + m_covariance.setZero(); + m_covariance.block<3, 3>(ePos0, ePos0) = covariance; +} + +void Acts::Vertex::setFullCovariance(const SquareMatrix4& covariance) { + m_covariance = covariance; +} + +void Acts::Vertex::setTracksAtVertex(std::vector tracks) { + m_tracksAtVertex = std::move(tracks); +} + +void Acts::Vertex::setFitQuality(double chiSquared, double numberDoF) { + m_chiSquared = chiSquared; + m_numberDoF = numberDoF; +} + +void Acts::Vertex::setFitQuality(std::pair fitQuality) { + m_chiSquared = fitQuality.first; + m_numberDoF = fitQuality.second; +} diff --git a/Core/include/Acts/Vertexing/ZScanVertexFinder.ipp b/Core/src/Vertexing/ZScanVertexFinder.cpp similarity index 83% rename from Core/include/Acts/Vertexing/ZScanVertexFinder.ipp rename to Core/src/Vertexing/ZScanVertexFinder.cpp index 76ea5e4c0ca..1960e7b65dc 100644 --- a/Core/include/Acts/Vertexing/ZScanVertexFinder.ipp +++ b/Core/src/Vertexing/ZScanVertexFinder.cpp @@ -1,15 +1,27 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -inline auto Acts::ZScanVertexFinder::find( +#include "Acts/Vertexing/ZScanVertexFinder.hpp" + +Acts::ZScanVertexFinder::ZScanVertexFinder(const Config& cfg, + std::unique_ptr logger) + : m_cfg(cfg), m_logger(std::move(logger)) { + if (!m_cfg.extractParameters.connected()) { + throw std::invalid_argument( + "ZScanVertexFinder: " + "No track parameter extractor provided."); + } +} + +Acts::Result> Acts::ZScanVertexFinder::find( const std::vector& trackVector, const VertexingOptions& vertexingOptions, - IVertexFinder::State& /*state*/) const -> Result> { + IVertexFinder::State& /*state*/) const { double ZResult = 0.; // Prepare the vector of points, on which the 3d mode has later to be // calculated @@ -95,11 +107,5 @@ inline auto Acts::ZScanVertexFinder::find( vertexingOptions.constraint.time()); Vertex vtxResult = Vertex(output); - // Vector to be filled with one single vertex - std::vector vertexCollection; - - // Add vertex to vertexCollection - vertexCollection.push_back(vtxResult); - - return vertexCollection; + return std::vector{vtxResult}; } diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.cpp index 9b1847cc193..3b28a9acb57 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.cpp @@ -9,50 +9,162 @@ #include "ActsExamples/TruthTracking/TruthVertexFinder.hpp" #include "ActsExamples/EventData/ProtoVertex.hpp" +#include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Utilities/GroupBy.hpp" #include "ActsExamples/Utilities/Range.hpp" +#include "ActsExamples/Validation/TrackClassification.hpp" #include "ActsFatras/EventData/Barcode.hpp" #include +#include #include #include +#include #include +#include namespace ActsExamples { -struct AlgorithmContext; -} // namespace ActsExamples -ActsExamples::TruthVertexFinder::TruthVertexFinder(const Config& config, - Acts::Logging::Level level) +namespace { + +using ConstTrackProxy = ConstTrackContainer::ConstTrackProxy; +using TrackIndex = ConstTrackProxy::IndexType; + +struct TrackMatchEntry { + std::optional particle; + + /// Number of hits on the track that are associated to a particle + /// Sorted by decreasing number of hits + std::vector contributingParticles; +}; + +struct ParticleMatchEntry { + std::optional track; + std::uint32_t duplicates{}; + std::uint32_t fakes{}; +}; + +} // namespace + +TruthVertexFinder::TruthVertexFinder(const Config& config, + Acts::Logging::Level level) : IAlgorithm("TruthVertexFinder", level), m_cfg(config) { + if (m_cfg.inputTracks.empty()) { + throw std::invalid_argument("Missing input tracks collection"); + } if (m_cfg.inputParticles.empty()) { throw std::invalid_argument("Missing input truth particles collection"); } + if (m_cfg.inputMeasurementParticlesMap.empty()) { + throw std::invalid_argument("Missing input hit-particles map collection"); + } if (m_cfg.outputProtoVertices.empty()) { throw std::invalid_argument("Missing output proto vertices collection"); } m_inputParticles.initialize(m_cfg.inputParticles); + m_inputTracks.initialize(m_cfg.inputTracks); + m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap); m_outputProtoVertices.initialize(m_cfg.outputProtoVertices); } -ActsExamples::ProcessCode ActsExamples::TruthVertexFinder::execute( - const AlgorithmContext& ctx) const { +ProcessCode TruthVertexFinder::execute(const AlgorithmContext& ctx) const { // prepare input and output collections - ACTS_VERBOSE("Reading particles from " << m_cfg.inputParticles); + const auto& tracks = m_inputTracks(ctx); const auto& particles = m_inputParticles(ctx); - ProtoVertexContainer protoVertices; + const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx); + ACTS_VERBOSE("Have " << particles.size() << " particles"); + ACTS_VERBOSE("Have " << tracks.size() << " tracks"); + + std::unordered_map trackParticleMatching; + std::unordered_map particleTrackMatching; + + { + // For each particle within a track, how many hits did it contribute + std::vector particleHitCounts; + + for (const auto& track : tracks) { + // Get the majority truth particle to this track + identifyContributingParticles(hitParticlesMap, track, particleHitCounts); + if (particleHitCounts.empty()) { + ACTS_DEBUG( + "No truth particle associated with this trajectory with tip index " + "= " + << track.tipIndex()); + continue; + } + + // Get the majority particleId and majority particle counts + // Note that the majority particle might be not in the truth seeds + // collection + ActsFatras::Barcode majorityParticleId = + particleHitCounts.front().particleId; + std::size_t nMajorityHits = particleHitCounts.front().hitCount; + + if (particles.find(majorityParticleId) == particles.end()) { + ACTS_DEBUG( + "The majority particle is not in the input particle collection, " + "majorityParticleId = " + << majorityParticleId); + continue; + } + + // Check if the trajectory is matched with truth. + // If not, it will be classified as 'fake' + const bool recoMatched = + static_cast(nMajorityHits) / track.nMeasurements() >= + m_cfg.trackMatchingRatio; + + if (recoMatched) { + trackParticleMatching[track.index()] = {majorityParticleId, + particleHitCounts}; + + auto& particleTrackMatch = particleTrackMatching[majorityParticleId]; + if (!particleTrackMatch.track) { + particleTrackMatch.track = track.index(); + } else { + // we already have a track associated with this particle and have to + // resolve the ambiguity. + // we will use the track with more hits and smaller chi2 + const auto& otherTrack = + tracks.getTrack(particleTrackMatch.track.value()); + if (otherTrack.nMeasurements() < track.nMeasurements() || + otherTrack.chi2() > track.chi2()) { + particleTrackMatch.track = track.index(); + } + + ++particleTrackMatch.duplicates; + } + } else { + trackParticleMatching[track.index()] = {std::nullopt, + particleHitCounts}; + + auto& particleTrackMatch = particleTrackMatching[majorityParticleId]; + ++particleTrackMatch.fakes; + } + } + } + + std::unordered_map> protoVertexTrackMap; + + for (const auto& [particle, trackMatch] : particleTrackMatching) { + auto vtxId = + SimBarcode(particle).setParticle(0).setGeneration(0).setSubParticle(0); + + protoVertexTrackMap[vtxId].push_back(trackMatch.track.value()); + } + + ProtoVertexContainer protoVertices; // assumes the begin/end iterator references the particles container - auto addProtoVertex = [&](SimParticleContainer::const_iterator begin, - const SimParticleContainer::const_iterator& end) { + auto addProtoVertex = [&](const std::vector& vtxTracks) { ProtoVertex protoVertex; - protoVertex.reserve(std::distance(begin, end)); - // determine each particle index - for (; begin != end; ++begin) { - protoVertex.push_back(std::distance(particles.begin(), begin)); + protoVertex.reserve(vtxTracks.size()); + for (const auto& track : vtxTracks) { + protoVertex.push_back(track); } protoVertices.push_back(std::move(protoVertex)); }; @@ -60,29 +172,42 @@ ActsExamples::ProcessCode ActsExamples::TruthVertexFinder::execute( if (m_cfg.excludeSecondaries) { // if secondaries are excluded, the `separateSecondaries` flag has no effect // since there will be no secondary vertices to separate - for (auto&& [vtxId, vtxParticles] : groupBySecondaryVertex(particles)) { + for (auto&& [vtxId, vtxTracks] : protoVertexTrackMap) { if (vtxId.vertexSecondary() != 0u) { continue; } - addProtoVertex(vtxParticles.begin(), vtxParticles.end()); + addProtoVertex(vtxTracks); } } else { // particles from secondary vertices should be included if (m_cfg.separateSecondaries) { // secondary particles are added to separate secondary vertices - for (auto&& [vtxId, vtxParticles] : groupBySecondaryVertex(particles)) { - addProtoVertex(vtxParticles.begin(), vtxParticles.end()); + for (auto&& [vtxId, vtxTracks] : protoVertexTrackMap) { + addProtoVertex(vtxTracks); } } else { // secondary particles are included in the primary vertex - for (auto&& [vtxId, vtxParticles] : groupByPrimaryVertex(particles)) { - addProtoVertex(vtxParticles.begin(), vtxParticles.end()); + + std::unordered_map> + protoVertexTrackMap2; + for (auto&& [vtxId, vtxTracks] : protoVertexTrackMap) { + auto vtxId2 = SimBarcode(vtxId).setVertexSecondary(0); + protoVertexTrackMap2[vtxId2].insert(protoVertexTrackMap2[vtxId].end(), + vtxTracks.begin(), vtxTracks.end()); + } + + for (auto&& [vtxId, vtxTracks] : protoVertexTrackMap2) { + addProtoVertex(vtxTracks); } } } ACTS_VERBOSE("Write " << protoVertices.size() << " proto vertex to " << m_cfg.outputProtoVertices); + m_outputProtoVertices(ctx, std::move(protoVertices)); + return ProcessCode::SUCCESS; } + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp index 005f9f2e423..da699caca71 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthVertexFinder.hpp @@ -11,6 +11,7 @@ #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/EventData/ProtoVertex.hpp" #include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IAlgorithm.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" @@ -23,15 +24,23 @@ struct AlgorithmContext; /// Group particles into proto vertices using truth information. class TruthVertexFinder final : public IAlgorithm { public: + using HitParticlesMap = ActsExamples::IndexMultimap; + struct Config { + /// The input tracks that should be used to create proto vertices. + std::string inputTracks; /// The input truth particles that should be used to create proto vertices. std::string inputParticles; + /// Input hit-particles map collection. + std::string inputMeasurementParticlesMap; /// The output proto vertices collection. std::string outputProtoVertices; /// Exclude secondary particles not originating from the primary vertex. - bool excludeSecondaries = false; + bool excludeSecondaries = true; /// Build separate proto vertices for the secondary particles. bool separateSecondaries = false; + /// The minimum number of tracks required to create a proto vertex. + double trackMatchingRatio = 0.5; }; TruthVertexFinder(const Config& config, Acts::Logging::Level level); @@ -44,8 +53,10 @@ class TruthVertexFinder final : public IAlgorithm { private: Config m_cfg; + ReadDataHandle m_inputTracks{this, "InputTracks"}; ReadDataHandle m_inputParticles{this, "InputParticles"}; - + ReadDataHandle m_inputMeasurementParticlesMap{ + this, "InputMeasurementParticlesMap"}; WriteDataHandle m_outputProtoVertices{ this, "OutputProtoVertices"}; }; diff --git a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp index 36a53f1e0f3..7abe177d913 100644 --- a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp +++ b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp @@ -89,14 +89,14 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm { const Config& config() const { return m_cfg; } private: - Acts::AdaptiveMultiVertexFinder makeVertexFinder() const; + Acts::AdaptiveMultiVertexFinder makeVertexFinder() const; Config m_cfg; std::shared_ptr m_propagator; Acts::ImpactPointEstimator m_ipEstimator; Linearizer m_linearizer; - Acts::AdaptiveMultiVertexFinder m_vertexFinder; + Acts::AdaptiveMultiVertexFinder m_vertexFinder; ReadDataHandle m_inputTrackParameters{ this, "InputTrackParameters"}; diff --git a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp index 97732f84ed3..384c16051db 100644 --- a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp +++ b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp @@ -54,8 +54,8 @@ class IterativeVertexFinderAlgorithm final : public IAlgorithm { using Propagator = Acts::Propagator>; using Linearizer = Acts::HelicalTrackLinearizer; using Fitter = Acts::FullBilloirVertexFitter; - using Seeder = Acts::TrackDensityVertexFinder; - using Finder = Acts::IterativeVertexFinder; + using Seeder = Acts::TrackDensityVertexFinder; + using Finder = Acts::IterativeVertexFinder; using Options = Acts::VertexingOptions; using VertexCollection = std::vector; diff --git a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp index fbc420c27e1..38e003576db 100644 --- a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp +++ b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp @@ -80,11 +80,10 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm:: } auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const - -> Acts::AdaptiveMultiVertexFinder { + -> Acts::AdaptiveMultiVertexFinder { std::shared_ptr seedFinder; if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) { - using Seeder = Acts::TrackDensityVertexFinder; - using Finder = Acts::AdaptiveMultiVertexFinder; + using Seeder = Acts::TrackDensityVertexFinder; Acts::GaussianTrackDensity::Config trkDensityCfg; trkDensityCfg.extractParameters .connect<&Acts::InputTrack::extractParameters>(); @@ -101,7 +100,6 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const // Set up vertex seeder and finder using Seeder = Acts::AdaptiveGridDensityVertexFinder; - using Finder = Acts::AdaptiveMultiVertexFinder; Seeder::Config seederConfig(trkDensity); seederConfig.extractParameters .connect<&Acts::InputTrack::extractParameters>(); @@ -110,8 +108,6 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const throw std::invalid_argument("Unknown seed finder"); } - using Finder = Acts::AdaptiveMultiVertexFinder; - // Set up deterministic annealing with user-defined temperatures Acts::AnnealingUtility::Config annealingConfig; annealingConfig.setOfTemperatures = {1.}; @@ -128,8 +124,8 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const Fitter fitter(std::move(fitterCfg), logger().cloneWithSuffix("AdaptiveMultiVertexFitter")); - Finder::Config finderConfig(std::move(fitter), seedFinder, m_ipEstimator, - m_cfg.bField); + Acts::AdaptiveMultiVertexFinder::Config finderConfig( + std::move(fitter), seedFinder, m_ipEstimator, m_cfg.bField); // Set the initial variance of the 4D vertex position. Since time is on a // numerical scale, we have to provide a greater value in the corresponding // dimension. @@ -155,7 +151,8 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const .template connect<&Acts::InputTrack::extractParameters>(); // Instantiate the finder - return Finder(std::move(finderConfig), logger().clone()); + return Acts::AdaptiveMultiVertexFinder(std::move(finderConfig), + logger().clone()); } ActsExamples::ProcessCode diff --git a/Examples/Io/Json/src/JsonDigitizationConfig.cpp b/Examples/Io/Json/src/JsonDigitizationConfig.cpp index 06462a461e0..f9aa0ba9249 100644 --- a/Examples/Io/Json/src/JsonDigitizationConfig.cpp +++ b/Examples/Io/Json/src/JsonDigitizationConfig.cpp @@ -100,7 +100,7 @@ void from_json( } else if (sType == "Digitial") { Acts::BinningData bd; from_json(j["bindata"], bd); - f = Digitization::Uniform(std::move(bd)); + f = Digitization::Digital(std::move(bd)); } else if (sType == "Exact") { f = Digitization::Exact{}; } else { diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp index dfeea124ae7..47f2fce0368 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp @@ -232,51 +232,31 @@ ActsExamples::ProcessCode ActsExamples::CKFPerformanceWriter::writeT( // Fill fake rate plots m_fakeRatePlotTool.fill(m_fakeRatePlotCache, fittedParameters, isFake); - // Use neural network classification for duplication rate plots - // Currently, the network used for this example can only handle - // good/duplicate classification, so need to manually exclude fake tracks - if (m_cfg.duplicatedPredictor && !isFake) { - inputFeatures[0] = track.nMeasurements(); - inputFeatures[1] = track.nOutliers(); - inputFeatures[2] = track.chi2() * 1.0 / track.nDoF(); - // predict if current trajectory is 'duplicate' - bool isDuplicated = m_cfg.duplicatedPredictor(inputFeatures); - // Add to number of duplicated particles - if (isDuplicated) { - m_nTotalDuplicateTracks++; - } - // Fill the duplication rate - m_duplicationPlotTool.fill(m_duplicationPlotCache, fittedParameters, - isDuplicated); - } // Counting number of total trajectories m_nTotalTracks++; } // Use truth-based classification for duplication rate plots - if (!m_cfg.duplicatedPredictor) { - // Loop over all truth-matched reco tracks for duplication rate plots - for (auto& [particleId, matchedTracks] : matched) { - // Sort the reco tracks matched to this particle by the number of majority - // hits - std::sort(matchedTracks.begin(), matchedTracks.end(), - [](const RecoTrackInfo& lhs, const RecoTrackInfo& rhs) { - return lhs.first > rhs.first; - }); - for (std::size_t itrack = 0; itrack < matchedTracks.size(); itrack++) { - const auto& [nMajorityHits, fittedParameters] = - matchedTracks.at(itrack); - // The tracks with maximum number of majority hits is taken as the - // 'real' track; others are as 'duplicated' - bool isDuplicated = (itrack != 0); - // the track is associated to the same particle - if (isDuplicated) { - m_nTotalDuplicateTracks++; - } - // Fill the duplication rate - m_duplicationPlotTool.fill(m_duplicationPlotCache, fittedParameters, - isDuplicated); + // Loop over all truth-matched reco tracks for duplication rate plots + for (auto& [particleId, matchedTracks] : matched) { + // Sort the reco tracks matched to this particle by the number of majority + // hits + std::sort(matchedTracks.begin(), matchedTracks.end(), + [](const RecoTrackInfo& lhs, const RecoTrackInfo& rhs) { + return lhs.first > rhs.first; + }); + for (std::size_t itrack = 0; itrack < matchedTracks.size(); itrack++) { + const auto& [nMajorityHits, fittedParameters] = matchedTracks.at(itrack); + // The tracks with maximum number of majority hits is taken as the + // 'real' track; others are as 'duplicated' + bool isDuplicated = (itrack != 0); + // the track is associated to the same particle + if (isDuplicated) { + m_nTotalDuplicateTracks++; } + // Fill the duplication rate + m_duplicationPlotTool.fill(m_duplicationPlotCache, fittedParameters, + isDuplicated); } } diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp index 8ab6cc44741..7d2cdd4d4f3 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp @@ -75,10 +75,6 @@ class CKFPerformanceWriter final : public WriterT { /// Write additional matching details to a TTree bool writeMatchingDetails = false; - - /// function to check if neural network predicted track label is - /// duplicate - std::function&)> duplicatedPredictor = nullptr; }; /// Construct from configuration and log level. diff --git a/Examples/Io/Root/src/RootParticleWriter.cpp b/Examples/Io/Root/src/RootParticleWriter.cpp index 96b27b9cb25..6be50871d1d 100644 --- a/Examples/Io/Root/src/RootParticleWriter.cpp +++ b/Examples/Io/Root/src/RootParticleWriter.cpp @@ -186,8 +186,9 @@ ActsExamples::ProcessCode ActsExamples::RootParticleWriter::writeT( // get the sim hits auto it = hitsPerParticle.find(particle.particleId()); if (it == hitsPerParticle.end()) { - ACTS_INFO("Could not find sim hits for " - << particle.particleId() << " in event " << ctx.eventNumber); + ACTS_DEBUG("Could not find sim hits for " + << particle.particleId() << " in event " << ctx.eventNumber); + m_numberOfHits.push_back(0); } else { // get the number of hits m_numberOfHits.push_back(it->second); diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 86e29fc9666..be188803890 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1803,7 +1803,9 @@ def addVertexFitting( if vertexFinder == VertexFinder.Truth: findVertices = TruthVertexFinder( level=customLogLevel(), + inputTracks=tracks, inputParticles=selectedParticles, + inputMeasurementParticlesMap="measurement_particles_map", outputProtoVertices=outputProtoVertices, excludeSecondaries=True, ) diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp index 194816da8ad..4227e478f57 100644 --- a/Examples/Python/src/Output.cpp +++ b/Examples/Python/src/Output.cpp @@ -380,9 +380,8 @@ void addOutput(Context& ctx) { ActsExamples::CKFPerformanceWriter, mex, "CKFPerformanceWriter", inputTracks, inputParticles, inputMeasurementParticlesMap, filePath, fileMode, effPlotToolConfig, fakeRatePlotToolConfig, - duplicationPlotToolConfig, trackSummaryPlotToolConfig, - duplicatedPredictor, truthMatchProbMin, doubleMatching, - writeMatchingDetails); + duplicationPlotToolConfig, trackSummaryPlotToolConfig, truthMatchProbMin, + doubleMatching, writeMatchingDetails); ACTS_PYTHON_DECLARE_WRITER( ActsExamples::RootNuclearInteractionParametersWriter, mex, diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index dc360d3bf9a..fe1e1bf162d 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -185,8 +185,9 @@ void addTruthTracking(Context& ctx) { } ACTS_PYTHON_DECLARE_ALGORITHM( - ActsExamples::TruthVertexFinder, mex, "TruthVertexFinder", inputParticles, - outputProtoVertices, excludeSecondaries, separateSecondaries); + ActsExamples::TruthVertexFinder, mex, "TruthVertexFinder", inputTracks, + inputParticles, inputMeasurementParticlesMap, outputProtoVertices, + excludeSecondaries, separateSecondaries, trackMatchingRatio); ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::TrackModifier, mex, "TrackModifier", inputTracks, outputTracks, diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 0afa67afe32..d926bbe28f8 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -1,7 +1,7 @@ test_pythia8__pythia8_particles.root: 082eb3dbf142929df28acfec48d61f260e132cb103966d143f430eeeda89daa1 test_fatras__particles_simulation.root: 4ce0feb0ecb234143f418789a7a0d377f4a653d529c49bf999dda1878c50cee1 test_fatras__hits.root: 2e47d9ba55fa1b377f70c361107fe811e9880d14c42cb3d7a9cd4616a6f33a54 -test_geant4__particles_simulation.root: 366281a37f60392aa8f27a33554d6c7b7f7fdc967a1f8c4c4d2e7426b0dbf78d +test_geant4__particles_simulation.root: 29fb3ed0c9ea48bc64cb6e4a83f4f6ca535415ab7a71051ca385943ce3dea31f test_geant4__hits.root: 1ed8ad8d3081980d8eb5b3cdd940f862f03c7683b7dda1cf6fe6ac900ac510be test_seeding__estimatedparams.root: 7aaf6f92aeced00aba73ffa7b44db76e1f6be9820f864d8a27646feee9e64544 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880 diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp index e0fefc7b364..3eef1364f7f 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp @@ -112,20 +112,16 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_test) { Fitter fitter(fitterCfg); - using SeedFinder = TrackDensityVertexFinder; - GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect<&InputTrack::extractParameters>(); - auto seedFinder = - std::make_shared(SeedFinder::Config{densityCfg}); - - using Finder = AdaptiveMultiVertexFinder; + auto seedFinder = std::make_shared( + TrackDensityVertexFinder::Config{densityCfg}); - Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator, - bField); + AdaptiveMultiVertexFinder::Config finderConfig(std::move(fitter), seedFinder, + ipEstimator, bField); finderConfig.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder(std::move(finderConfig)); + AdaptiveMultiVertexFinder finder(std::move(finderConfig)); IVertexFinder::State state = finder.makeState(magFieldContext); auto csvData = readTracksAndVertexCSV(toolString); @@ -274,20 +270,16 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_usertype_test) { Fitter fitter(fitterCfg); - using SeedFinder = TrackDensityVertexFinder; - GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect(extractParameters); - auto seedFinder = - std::make_shared(SeedFinder::Config{densityCfg}); + auto seedFinder = std::make_shared( + TrackDensityVertexFinder::Config{densityCfg}); - using Finder = AdaptiveMultiVertexFinder; - - Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator, - bField); + AdaptiveMultiVertexFinder::Config finderConfig( + std::move(fitter), std::move(seedFinder), ipEstimator, bField); finderConfig.extractParameters.connect(extractParameters); - Finder finder(std::move(finderConfig)); + AdaptiveMultiVertexFinder finder(std::move(finderConfig)); IVertexFinder::State state = finder.makeState(magFieldContext); auto csvData = readTracksAndVertexCSV(toolString); @@ -428,13 +420,11 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_grid_seed_finder_test) { auto seedFinder = std::make_shared(seedFinderCfg); - using Finder = AdaptiveMultiVertexFinder; - - Finder::Config finderConfig(std::move(fitter), std::move(seedFinder), ipEst, - bField); + AdaptiveMultiVertexFinder::Config finderConfig( + std::move(fitter), std::move(seedFinder), ipEst, bField); finderConfig.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder(std::move(finderConfig)); + AdaptiveMultiVertexFinder finder(std::move(finderConfig)); IVertexFinder::State state = finder.makeState(magFieldContext); auto csvData = readTracksAndVertexCSV(toolString); @@ -587,13 +577,11 @@ BOOST_AUTO_TEST_CASE( auto seedFinder = std::make_shared(seedFinderCfg); - using Finder = AdaptiveMultiVertexFinder; - - Finder::Config finderConfig(std::move(fitter), std::move(seedFinder), ipEst, - bField); + AdaptiveMultiVertexFinder::Config finderConfig( + std::move(fitter), std::move(seedFinder), ipEst, bField); finderConfig.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder(std::move(finderConfig)); + AdaptiveMultiVertexFinder finder(std::move(finderConfig)); IVertexFinder::State state = finder.makeState(magFieldContext); auto csvData = readTracksAndVertexCSV(toolString); diff --git a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt index c0a7ac4fa0c..ab7a2873500 100644 --- a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt +++ b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt @@ -1,7 +1,6 @@ add_unittest(FullBilloirVertexFitter FullBilloirVertexFitterTests.cpp) add_unittest(ImpactPointEstimator ImpactPointEstimatorTests.cpp) add_unittest(IterativeVertexFinder IterativeVertexFinderTests.cpp) -add_unittest(KalmanVertexTrackUpdater KalmanVertexTrackUpdaterTests.cpp) add_unittest(KalmanVertexUpdater KalmanVertexUpdaterTests.cpp) add_unittest(LinearizedTrackFactory LinearizedTrackFactoryTests.cpp) add_unittest(AdaptiveMultiVertexFitter AdaptiveMultiVertexFitterTests.cpp) diff --git a/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp index eed6c0fbfe9..f14fe5b180b 100644 --- a/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp @@ -161,18 +161,18 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test) { auto sFinder = std::make_shared(seedFinderCfg); // Vertex Finder - using VertexFinder = IterativeVertexFinder; - VertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), - ipEstimator); + IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), + ipEstimator); cfg.field = bField; cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer); cfg.reassignTracksAfterFirstFit = true; cfg.extractParameters.connect<&InputTrack::extractParameters>(); - VertexFinder finder(std::move(cfg)); - IVertexFinder::State state{VertexFinder::State(*bField, magFieldContext)}; + IterativeVertexFinder finder(std::move(cfg)); + IVertexFinder::State state{ + IterativeVertexFinder::State(*bField, magFieldContext)}; // Vector to be filled with all tracks in current event std::vector> tracks; @@ -383,15 +383,16 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test_user_track_type) { auto sFinder = std::make_shared(seedFinderCfg); // Vertex Finder - using VertexFinder = IterativeVertexFinder; - VertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), - ipEstimator); + IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), + ipEstimator); cfg.reassignTracksAfterFirstFit = true; cfg.extractParameters.connect(extractParameters); cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer); cfg.field = bField; - VertexFinder finder(std::move(cfg)); - IVertexFinder::State state{VertexFinder::State(*bField, magFieldContext)}; + + IterativeVertexFinder finder(std::move(cfg)); + IVertexFinder::State state{ + IterativeVertexFinder::State(*bField, magFieldContext)}; // Same for user track type tracks std::vector> tracks; @@ -585,10 +586,8 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test_athena_reference) { auto sFinder = std::make_shared(seedFinderCfg); - // Vertex Finder - using VertexFinder = IterativeVertexFinder; - - VertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), ipEstimator); + IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder), + ipEstimator); cfg.maxVertices = 200; cfg.maximumChi2cutForSeeding = 49; cfg.significanceCutSeeding = 12; @@ -596,8 +595,9 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test_athena_reference) { cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer); cfg.field = bField; - VertexFinder finder(std::move(cfg)); - IVertexFinder::State state{VertexFinder::State(*bField, magFieldContext)}; + IterativeVertexFinder finder(std::move(cfg)); + IVertexFinder::State state{ + IterativeVertexFinder::State(*bField, magFieldContext)}; auto csvData = readTracksAndVertexCSV(toolString); auto tracks = std::get(csvData); diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp deleted file mode 100644 index baa2828391c..00000000000 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp +++ /dev/null @@ -1,211 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include -#include -#include - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/Direction.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/GenericBoundTrackParameters.hpp" -#include "Acts/EventData/TrackParameters.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Geometry/GeometryIdentifier.hpp" -#include "Acts/MagneticField/ConstantBField.hpp" -#include "Acts/MagneticField/MagneticFieldContext.hpp" -#include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/Propagator.hpp" -#include "Acts/Surfaces/PerigeeSurface.hpp" -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/Result.hpp" -#include "Acts/Vertexing/HelicalTrackLinearizer.hpp" -#include "Acts/Vertexing/ImpactPointEstimator.hpp" -#include "Acts/Vertexing/KalmanVertexTrackUpdater.hpp" -#include "Acts/Vertexing/LinearizedTrack.hpp" -#include "Acts/Vertexing/TrackAtVertex.hpp" -#include "Acts/Vertexing/Vertex.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace bdata = boost::unit_test::data; -using namespace Acts::UnitLiterals; - -namespace Acts { -namespace Test { - -using Covariance = BoundSquareMatrix; -using Propagator = Acts::Propagator>; -using Linearizer = HelicalTrackLinearizer; - -// Create a test context -GeometryContext geoContext = GeometryContext(); -MagneticFieldContext magFieldContext = MagneticFieldContext(); - -// Vertex x/y position distribution -std::uniform_real_distribution vXYDist(-0.1_mm, 0.1_mm); -// Vertex z position distribution -std::uniform_real_distribution vZDist(-20_mm, 20_mm); -// Track d0 distribution -std::uniform_real_distribution d0Dist(-0.01_mm, 0.01_mm); -// Track z0 distribution -std::uniform_real_distribution z0Dist(-0.2_mm, 0.2_mm); -// Track pT distribution -std::uniform_real_distribution pTDist(0.4_GeV, 10_GeV); -// Track phi distribution -std::uniform_real_distribution phiDist(-M_PI, M_PI); -// Track theta distribution -std::uniform_real_distribution thetaDist(1.0, M_PI - 1.0); -// Track charge helper distribution -std::uniform_real_distribution qDist(-1, 1); -// Track IP resolution distribution -std::uniform_real_distribution resIPDist(0., 100_um); -// Track angular distribution -std::uniform_real_distribution resAngDist(0., 0.1); -// Track q/p resolution distribution -std::uniform_real_distribution resQoPDist(-0.01, 0.01); - -/// -/// @brief Unit test for KalmanVertexTrackUpdater -/// -BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { - bool debug = true; - - // Number of tests - unsigned int nTests = 10; - - // Set up RNG - int mySeed = 31415; - std::mt19937 gen(mySeed); - - // Set up constant B-Field - auto bField = std::make_shared(Vector3{0.0, 0.0, 1_T}); - - // Set up Eigenstepper - EigenStepper<> stepper(bField); - - // Set up propagator with void navigator - auto propagator = std::make_shared(stepper); - - // Set up ImpactPointEstimator, used for comparisons later - ImpactPointEstimator::Config ip3dEstConfig(bField, propagator); - ImpactPointEstimator ip3dEst(ip3dEstConfig); - ImpactPointEstimator::State state(bField->makeCache(magFieldContext)); - - // Set up HelicalTrackLinearizer, needed for linearizing the tracks - // Linearizer for BoundTrackParameters type test - Linearizer::Config ltConfig; - ltConfig.bField = bField; - ltConfig.propagator = propagator; - Linearizer linearizer(ltConfig); - auto fieldCache = bField->makeCache(magFieldContext); - - // Create perigee surface at origin - std::shared_ptr perigeeSurface = - Surface::makeShared(Vector3(0., 0., 0.)); - - // Create random tracks around origin and a random vertex. - // Update tracks with the assumption that they originate from - // the vertex position and check if they are closer to the - // vertex after the update process - for (unsigned int i = 0; i < nTests; ++i) { - // Construct positive or negative charge randomly - double q = qDist(gen) < 0 ? -1. : 1.; - - // Construct random track parameters - BoundTrackParameters::ParametersVector paramVec; - - paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen), - q / pTDist(gen), 0.; - - if (debug) { - std::cout << "Creating track parameters: " << paramVec << std::endl; - } - - // Fill vector of track objects with simple covariance matrix - Covariance covMat; - - // Resolutions - double res_d0 = resIPDist(gen); - double res_z0 = resIPDist(gen); - double res_ph = resAngDist(gen); - double res_th = resAngDist(gen); - double res_qp = resQoPDist(gen); - - covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0., 0., - 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0., - res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0., 0., - 0., 0., 0., 1.; - BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat), - ParticleHypothesis::pion()); - - std::shared_ptr perigee = - Surface::makeShared(Vector3::Zero()); - - // Linearized state of the track - LinearizedTrack linTrack = - linearizer - .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext, - fieldCache) - .value(); - - // Create TrackAtVertex - TrackAtVertex trkAtVtx(0., params, InputTrack{¶ms}); - - // Set linearized state of trackAtVertex - trkAtVtx.linearizedState = linTrack; - trkAtVtx.isLinearized = true; - - // Copy parameters for later comparison of old and new version - auto fittedParamsCopy = trkAtVtx.fittedParams; - - // Create a vertex - Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen)); - Vertex vtx(vtxPos); - - // Update trkAtVertex with assumption of originating from vtx - KalmanVertexTrackUpdater::update(trkAtVtx, vtx.fullPosition(), - vtx.fullCovariance(), 3); - - // The old distance - double oldDistance = - ip3dEst.calculateDistance(geoContext, fittedParamsCopy, vtxPos, state) - .value(); - - // The new distance after update - double newDistance = - ip3dEst - .calculateDistance(geoContext, trkAtVtx.fittedParams, vtxPos, state) - .value(); - if (debug) { - std::cout << "Old distance: " << oldDistance << std::endl; - std::cout << "New distance: " << newDistance << std::endl; - } - - // Parameters should have changed - BOOST_CHECK_NE(fittedParamsCopy, trkAtVtx.fittedParams); - - // After update, track should be closer to the vertex - BOOST_CHECK_LT(newDistance, oldDistance); - - } // end for loop - -} // end test case - -} // namespace Test -} // namespace Acts diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp index 5dbf8016eae..0a50341fd37 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp @@ -26,6 +26,7 @@ #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/HelicalTrackLinearizer.hpp" +#include "Acts/Vertexing/ImpactPointEstimator.hpp" #include "Acts/Vertexing/KalmanVertexUpdater.hpp" #include "Acts/Vertexing/LinearizedTrack.hpp" #include "Acts/Vertexing/TrackAtVertex.hpp" @@ -174,7 +175,7 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) { vtx.setFullCovariance(SquareMatrix4::Identity() * 0.01); // Update trkAtVertex with assumption of originating from vtx - KalmanVertexUpdater::updateVertexWithTrack<3>(vtx, trkAtVtx); + KalmanVertexUpdater::updateVertexWithTrack(vtx, trkAtVtx, 3); if (debug) { std::cout << "Old vertex position: " << vtxPos << std::endl; @@ -202,5 +203,131 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) { } // end test case +/// +/// @brief Unit test for KalmanVertexTrackUpdater +/// +BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { + bool debug = true; + + // Number of tests + unsigned int nTests = 10; + + // Set up RNG + int mySeed = 31415; + std::mt19937 gen(mySeed); + + // Set up constant B-Field + auto bField = std::make_shared(Vector3{0.0, 0.0, 1_T}); + + // Set up Eigenstepper + EigenStepper<> stepper(bField); + + // Set up propagator with void navigator + auto propagator = std::make_shared(stepper); + + // Set up ImpactPointEstimator, used for comparisons later + ImpactPointEstimator::Config ip3dEstConfig(bField, propagator); + ImpactPointEstimator ip3dEst(ip3dEstConfig); + ImpactPointEstimator::State state(bField->makeCache(magFieldContext)); + + // Set up HelicalTrackLinearizer, needed for linearizing the tracks + // Linearizer for BoundTrackParameters type test + Linearizer::Config ltConfig; + ltConfig.bField = bField; + ltConfig.propagator = propagator; + Linearizer linearizer(ltConfig); + auto fieldCache = bField->makeCache(magFieldContext); + + // Create perigee surface at origin + std::shared_ptr perigeeSurface = + Surface::makeShared(Vector3(0., 0., 0.)); + + // Create random tracks around origin and a random vertex. + // Update tracks with the assumption that they originate from + // the vertex position and check if they are closer to the + // vertex after the update process + for (unsigned int i = 0; i < nTests; ++i) { + // Construct positive or negative charge randomly + double q = qDist(gen) < 0 ? -1. : 1.; + + // Construct random track parameters + BoundTrackParameters::ParametersVector paramVec; + + paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen), + q / pTDist(gen), 0.; + + if (debug) { + std::cout << "Creating track parameters: " << paramVec << std::endl; + } + + // Fill vector of track objects with simple covariance matrix + Covariance covMat; + + // Resolutions + double res_d0 = resIPDist(gen); + double res_z0 = resIPDist(gen); + double res_ph = resAngDist(gen); + double res_th = resAngDist(gen); + double res_qp = resQoPDist(gen); + + covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0., 0., + 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0., + res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0., 0., + 0., 0., 0., 1.; + BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat), + ParticleHypothesis::pion()); + + std::shared_ptr perigee = + Surface::makeShared(Vector3::Zero()); + + // Linearized state of the track + LinearizedTrack linTrack = + linearizer + .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext, + fieldCache) + .value(); + + // Create TrackAtVertex + TrackAtVertex trkAtVtx(0., params, InputTrack{¶ms}); + + // Set linearized state of trackAtVertex + trkAtVtx.linearizedState = linTrack; + trkAtVtx.isLinearized = true; + + // Copy parameters for later comparison of old and new version + auto fittedParamsCopy = trkAtVtx.fittedParams; + + // Create a vertex + Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen)); + Vertex vtx(vtxPos); + + // Update trkAtVertex with assumption of originating from vtx + KalmanVertexUpdater::updateTrackWithVertex(trkAtVtx, vtx, 3); + + // The old distance + double oldDistance = + ip3dEst.calculateDistance(geoContext, fittedParamsCopy, vtxPos, state) + .value(); + + // The new distance after update + double newDistance = + ip3dEst + .calculateDistance(geoContext, trkAtVtx.fittedParams, vtxPos, state) + .value(); + if (debug) { + std::cout << "Old distance: " << oldDistance << std::endl; + std::cout << "New distance: " << newDistance << std::endl; + } + + // Parameters should have changed + BOOST_CHECK_NE(fittedParamsCopy, trkAtVtx.fittedParams); + + // After update, track should be closer to the vertex + BOOST_CHECK_LT(newDistance, oldDistance); + + } // end for loop + +} // end test case + } // namespace Test } // namespace Acts diff --git a/Tests/UnitTests/Core/Vertexing/TrackDensityVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/TrackDensityVertexFinderTests.cpp index f144b7f95cc..44249fb1cc6 100644 --- a/Tests/UnitTests/Core/Vertexing/TrackDensityVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/TrackDensityVertexFinderTests.cpp @@ -70,10 +70,9 @@ BOOST_AUTO_TEST_CASE(track_density_finder_test) { Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV}; VertexingOptions vertexingOptions(geoContext, magFieldContext); - using Finder = TrackDensityVertexFinder; GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder{{{densityCfg}}}; + TrackDensityVertexFinder finder{{{densityCfg}}}; auto state = finder.makeState(magFieldContext); // Start creating some track parameters @@ -150,10 +149,9 @@ BOOST_AUTO_TEST_CASE(track_density_finder_constr_test) { // Finder options VertexingOptions vertexingOptions(geoContext, magFieldContext, constraint); - using Finder = TrackDensityVertexFinder; GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder{{{densityCfg}}}; + TrackDensityVertexFinder finder{{{densityCfg}}}; auto state = finder.makeState(magFieldContext); // Start creating some track parameters @@ -227,10 +225,9 @@ BOOST_AUTO_TEST_CASE(track_density_finder_random_test) { Surface::makeShared(pos0); VertexingOptions vertexingOptions(geoContext, magFieldContext); - using Finder = TrackDensityVertexFinder; GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect<&InputTrack::extractParameters>(); - Finder finder{{{densityCfg}}}; + TrackDensityVertexFinder finder{{{densityCfg}}}; auto state = finder.makeState(magFieldContext); int mySeed = 31415; @@ -328,11 +325,9 @@ BOOST_AUTO_TEST_CASE(track_density_finder_usertrack_test) { return params.as()->parameters(); }; - using Finder = TrackDensityVertexFinder; - GaussianTrackDensity::Config densityCfg; densityCfg.extractParameters.connect(extractParameters); - Finder finder{{{densityCfg}}}; + TrackDensityVertexFinder finder{{{densityCfg}}}; auto state = finder.makeState(magFieldContext); // Start creating some track parameters diff --git a/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp b/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp index 2f0a5dfa911..89f93b5ac18 100644 --- a/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp +++ b/Tests/UnitTests/Core/Vertexing/VertexingDataHelper.hpp @@ -9,6 +9,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Units.hpp" +#include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Tests/CommonHelpers/DataDirectory.hpp" #include "Acts/Vertexing/Vertex.hpp" diff --git a/docs/getting_started.md b/docs/getting_started.md index 56b550d6458..5bd2effa723 100644 --- a/docs/getting_started.md +++ b/docs/getting_started.md @@ -267,6 +267,8 @@ components. | ACTS_BUILD_PLUGIN_PODIO | Build Podio plugin
type: `bool`, default: `OFF` | | ACTS_BUILD_PLUGIN_EDM4HEP | Build EDM4hep plugin
type: `bool`, default: `OFF` | | ACTS_BUILD_PLUGIN_FPEMON | Build FPE monitoring plugin
type: `bool`, default: `OFF` | +| ACTS_BUILD_PLUGIN_GEOMODEL | Build GeoModel plugin
type: `bool`, default: `OFF` | +| ACTS_USE_SYSTEM_GEOMODEL | Use a system-provided GeoModel
installation
type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF` | | ACTS_BUILD_PLUGIN_GEANT4 | Build Geant4 plugin
type: `bool`, default: `OFF` | | ACTS_BUILD_PLUGIN_EXATRKX | Build the Exa.TrkX plugin
type: `bool`, default: `OFF` | | ACTS_EXATRKX_ENABLE_ONNX | Build the Onnx backend for the exatrkx
plugin
type: `bool`, default: `OFF` | diff --git a/thirdparty/GeoModel/CMakeLists.txt b/thirdparty/GeoModel/CMakeLists.txt new file mode 100644 index 00000000000..87726a3ce12 --- /dev/null +++ b/thirdparty/GeoModel/CMakeLists.txt @@ -0,0 +1,28 @@ +# This file is part of the Acts project. +# +# Copyright (C) 2024 CERN for the benefit of the Acts project +# +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. + +# CMake include(s). +cmake_minimum_required( VERSION 3.11 ) +include( FetchContent ) + +# Tell the user what's happening. +message( STATUS "Building GeoModel as part of the ACTS project" ) + +set( GEOMODEL_VERSION "${_acts_geomodel_version}") + +# Declare where to get VecMem from. +set( ACTS_GEOMODEL_GIT_REPOSITORY "https://gitlab.cern.ch/GeoModelDev/GeoModel" + CACHE STRING "Git repository to take GeoModel from" ) +set( ACTS_GEOMODEL_GIT_TAG "${GEOMODEL_VERSION}" CACHE STRING "Version of GeoModel to build" ) +mark_as_advanced( ACTS_GEOMODEL_GIT_REPOSITORY ACTS_GEOMODEL_GIT_TAG ) +FetchContent_Declare( geomodel + GIT_REPOSITORY "${ACTS_GEOMODEL_GIT_REPOSITORY}" + GIT_TAG "${ACTS_GEOMODEL_GIT_TAG}" ) + +# Now set up its build. +FetchContent_MakeAvailable( geomodel ) diff --git a/thirdparty/GeoModel/README.md b/thirdparty/GeoModel/README.md new file mode 100644 index 00000000000..663446149eb --- /dev/null +++ b/thirdparty/GeoModel/README.md @@ -0,0 +1,5 @@ +# Build Recipe for GeoModel + +This directory holds a simple build recipe for the +[GeoModel](https://gitlab.cern.ch/GeoModelDev/GeoModel) project. It is used + in case `ACTS_USE_SYSTEM_GEOMODEL` is set to `FALSE` for the build. \ No newline at end of file