From 82b7a2d467a36c23a7fe5a54e2087ed62a93e9f3 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 13 Nov 2023 21:39:21 +0100 Subject: [PATCH 01/16] refactor --- .../Acts/Propagator/MultiEigenStepperLoop.hpp | 1 - .../Acts/TrackFitting/GaussianSumFitter.hpp | 4 +-- .../Acts/TrackFitting/GsfMixtureReduction.hpp | 8 ++++- Core/include/Acts/TrackFitting/GsfOptions.hpp | 11 +++++-- .../Acts/TrackFitting/detail/GsfActor.hpp | 19 ++++++------ .../detail/GsfComponentMerging.hpp} | 21 +++++-------- .../detail/SymmetricKlDistanceMatrix.hpp | 2 +- .../TrackFitting/TrackFitterFunction.hpp | 19 +++++++----- .../TrackFitting/src/GsfFitterFunction.cpp | 28 +++++++++-------- .../python/acts/examples/reconstruction.py | 5 ++-- Examples/Python/src/TrackFitting.cpp | 30 ++++++++++--------- .../TrackFitting/GsfComponentMergingTests.cpp | 2 +- 12 files changed, 83 insertions(+), 67 deletions(-) rename Core/include/Acts/{Utilities/GaussianMixtureReduction.hpp => TrackFitting/detail/GsfComponentMerging.hpp} (94%) diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp index 50323b77fb9..7fd4eb4766e 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp @@ -25,7 +25,6 @@ #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/detail/LoopStepperUtils.hpp" #include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/GaussianMixtureReduction.hpp" #include "Acts/Utilities/Intersection.hpp" #include "Acts/Utilities/Result.hpp" diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index d2bdb76a49d..030331650f8 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -449,9 +449,9 @@ struct GaussianSumFitter { if (options.referenceSurface) { const auto& params = *bwdResult->endParameters; - const auto [finalPars, finalCov] = Acts::reduceGaussianMixture( + const auto [finalPars, finalCov] = detail::mergeGaussianMixture( params.components(), params.referenceSurface(), - options.stateReductionMethod, [](auto& t) { + options.componentMergeMethod, [](auto& t) { return std::tie(std::get<0>(t), std::get<1>(t), *std::get<2>(t)); }); diff --git a/Core/include/Acts/TrackFitting/GsfMixtureReduction.hpp b/Core/include/Acts/TrackFitting/GsfMixtureReduction.hpp index 3b0a3d9704f..0624c948776 100644 --- a/Core/include/Acts/TrackFitting/GsfMixtureReduction.hpp +++ b/Core/include/Acts/TrackFitting/GsfMixtureReduction.hpp @@ -13,8 +13,14 @@ namespace Acts { +/// Greedy component reduction algorithm. Reduces the components with the +/// minimal symmetric KL-distance (applied only to the q/p-dimension) until the +/// required number of components is reached. +/// @param cmpCache the component collection +/// @param maxCmpsAfterMerge the number of components we want to reach +/// @param surface the surface type on which the components are void reduceMixtureWithKLDistance(std::vector &cmpCache, std::size_t maxCmpsAfterMerge, const Surface &surface); -} +} // namespace Acts diff --git a/Core/include/Acts/TrackFitting/GsfOptions.hpp b/Core/include/Acts/TrackFitting/GsfOptions.hpp index 9930c534d0c..0cf35fb11e6 100644 --- a/Core/include/Acts/TrackFitting/GsfOptions.hpp +++ b/Core/include/Acts/TrackFitting/GsfOptions.hpp @@ -20,6 +20,14 @@ namespace Acts { +/// @enum ComponentMergeMethod +/// +/// Available reduction methods for the reduction of a Gaussian mixture +enum class ComponentMergeMethod { eMean, eMaxWeight }; + +/// @struct GsfComponent +/// +/// Encapsulates a component of a Gaussian mixture as used by the GSF struct GsfComponent { ActsScalar weight = 0; BoundVector boundPars = BoundVector::Zero(); @@ -99,8 +107,7 @@ struct GsfOptions { std::string_view finalMultiComponentStateColumn = ""; - MixtureReductionMethod stateReductionMethod = - MixtureReductionMethod::eMaxWeight; + ComponentMergeMethod componentMergeMethod = ComponentMergeMethod::eMaxWeight; #if __cplusplus < 202002L GsfOptions() = delete; diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 0ec658a0ad6..3800d3d63be 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -20,6 +20,7 @@ #include "Acts/TrackFitting/GsfError.hpp" #include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/TrackFitting/KalmanFitter.hpp" +#include "Acts/TrackFitting/detail/GsfComponentMerging.hpp" #include "Acts/TrackFitting/detail/GsfUtils.hpp" #include "Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp" #include "Acts/Utilities/Zip.hpp" @@ -114,7 +115,7 @@ struct GsfActor { bool inReversePass = false; /// How to reduce the states that are stored in the multi trajectory - MixtureReductionMethod reductionMethod = MixtureReductionMethod::eMaxWeight; + ComponentMergeMethod mergeMethod = ComponentMergeMethod::eMaxWeight; const Logger* logger{nullptr}; @@ -714,15 +715,15 @@ struct GsfActor { proxy.setReferenceSurface(surface.getSharedPtr()); proxy.copyFrom(firstCmpProxy, mask); - auto [prtMean, prtCov] = reduceGaussianMixture( - tmpStates.tips, surface, m_cfg.reductionMethod, - PrtProjector{tmpStates.traj, tmpStates.weights}); + auto [prtMean, prtCov] = + mergeGaussianMixture(tmpStates.tips, surface, m_cfg.mergeMethod, + PrtProjector{tmpStates.traj, tmpStates.weights}); proxy.predicted() = prtMean; proxy.predictedCovariance() = prtCov; if (isMeasurement) { - auto [fltMean, fltCov] = reduceGaussianMixture( - tmpStates.tips, surface, m_cfg.reductionMethod, + auto [fltMean, fltCov] = mergeGaussianMixture( + tmpStates.tips, surface, m_cfg.mergeMethod, FltProjector{tmpStates.traj, tmpStates.weights}); proxy.filtered() = fltMean; proxy.filteredCovariance() = fltCov; @@ -744,8 +745,8 @@ struct GsfActor { result.surfacesVisitedBwdAgain.push_back(&surface); if (trackState.hasSmoothed()) { - const auto [smtMean, smtCov] = reduceGaussianMixture( - tmpStates.tips, surface, m_cfg.reductionMethod, + const auto [smtMean, smtCov] = mergeGaussianMixture( + tmpStates.tips, surface, m_cfg.mergeMethod, FltProjector{tmpStates.traj, tmpStates.weights}); trackState.smoothed() = smtMean; @@ -766,7 +767,7 @@ struct GsfActor { m_cfg.abortOnError = options.abortOnError; m_cfg.disableAllMaterialHandling = options.disableAllMaterialHandling; m_cfg.weightCutoff = options.weightCutoff; - m_cfg.reductionMethod = options.stateReductionMethod; + m_cfg.mergeMethod = options.componentMergeMethod; m_cfg.calibrationContext = &options.calibrationContext.get(); } }; diff --git a/Core/include/Acts/Utilities/GaussianMixtureReduction.hpp b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp similarity index 94% rename from Core/include/Acts/Utilities/GaussianMixtureReduction.hpp rename to Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp index 6811ebcfb14..cf5bc4bf1b0 100644 --- a/Core/include/Acts/Utilities/GaussianMixtureReduction.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/Utilities/Identity.hpp" #include "Acts/Utilities/detail/periodic.hpp" @@ -18,8 +19,7 @@ #include #include -namespace Acts { -namespace detail { +namespace Acts::detail { /// Angle descriptions for the combineBoundGaussianMixture function template @@ -203,13 +203,6 @@ auto gaussianMixtureMeanCov(const components_t components, return RetType{mean, cov}; } -} // namespace detail - -/// @enum MixtureReductionMethod -/// -/// Available reduction methods for the reduction of a Gaussian mixture -enum class MixtureReductionMethod { eMean, eMaxWeight }; - /// Reduce Gaussian mixture /// /// @param mixture The mixture iterable @@ -220,16 +213,16 @@ enum class MixtureReductionMethod { eMean, eMaxWeight }; /// /// @return parameters and covariance as std::tuple< BoundVector, BoundMatrix > template -auto reduceGaussianMixture(const mixture_t &mixture, const Surface &surface, - MixtureReductionMethod method, - projector_t &&projector = projector_t{}) { +auto mergeGaussianMixture(const mixture_t &mixture, const Surface &surface, + ComponentMergeMethod method, + projector_t &&projector = projector_t{}) { using R = std::tuple; const auto [mean, cov] = detail::angleDescriptionSwitch(surface, [&](const auto &desc) { return detail::gaussianMixtureMeanCov(mixture, projector, desc); }); - if (method == MixtureReductionMethod::eMean) { + if (method == ComponentMergeMethod::eMean) { return R{mean, cov}; } else { const auto maxWeightIt = std::max_element( @@ -242,4 +235,4 @@ auto reduceGaussianMixture(const mixture_t &mixture, const Surface &surface, } } -} // namespace Acts +} // namespace Acts::detail diff --git a/Core/include/Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp b/Core/include/Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp index 3c511091e95..a4ab04e4ffa 100644 --- a/Core/include/Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp +++ b/Core/include/Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp @@ -11,8 +11,8 @@ #include "Acts/EventData/MultiComponentTrackParameters.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/TrackParameters.hpp" +#include "Acts/TrackFitting/detail/GsfComponentMerging.hpp" #include "Acts/TrackFitting/detail/GsfUtils.hpp" -#include "Acts/Utilities/GaussianMixtureReduction.hpp" namespace Acts::detail { diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp index 0e8d4957695..f126c55a918 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp @@ -16,9 +16,9 @@ #include "Acts/Geometry/TrackingGeometry.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" #include "Acts/MagneticField/MagneticFieldProvider.hpp" -#include "Acts/Propagator/MultiEigenStepperLoop.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/TrackFitting/BetheHeitlerApprox.hpp" +#include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/Utilities/CalibrationContext.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/EventData/MeasurementCalibration.hpp" @@ -74,22 +74,27 @@ std::shared_ptr makeKalmanFitterFunction( /// approximation using BetheHeitlerApprox = Acts::AtlasBetheHeitlerApprox<6, 5>; +/// Available algorithms for the mixture reduction +enum class MixtureReductionAlgorithm { KLDistance }; + /// Makes a fitter function object for the GSF /// /// @param trackingGeometry the trackingGeometry for the propagator /// @param magneticField the magnetic field for the propagator /// @param betheHeitlerApprox The object that encapsulates the approximation. /// @param maxComponents number of maximum components in the track state -/// @param abortOnError whether to call std::abort on an error -/// @param disableAllMaterialHandling run the GSF like a KF (no energy loss, -/// always 1 component, ...) +/// @param weightCutoff when to drop components +/// @param componentMergeMethod How to merge a mixture to a single set of +/// parameters and covariance +/// @param mixtureReductionAlgorithm How to reduce the number of components +/// in a mixture /// @param logger a logger instance std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - BetheHeitlerApprox betheHeitlerApprox, size_t maxComponents, - double weightCutoff, Acts::MixtureReductionMethod finalReductionMethod, - bool abortOnError, bool disableAllMaterialHandling, + BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, + double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger); /// Makes a fitter function object for the Global Chi Square Fitter (GX2F) diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index d71d7baffb9..b90b2e2377d 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -24,7 +24,6 @@ #include "Acts/TrackFitting/GsfMixtureReduction.hpp" #include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/Utilities/Delegate.hpp" -#include "Acts/Utilities/GaussianMixtureReduction.hpp" #include "Acts/Utilities/HashedString.hpp" #include "Acts/Utilities/Intersection.hpp" #include "Acts/Utilities/Logger.hpp" @@ -81,8 +80,10 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { double weightCutoff = 0; bool abortOnError = false; bool disableAllMaterialHandling = false; - Acts::MixtureReductionMethod reductionMethod = - Acts::MixtureReductionMethod::eMaxWeight; + MixtureReductionAlgorithm reductionAlg = + MixtureReductionAlgorithm::KLDistance; + Acts::ComponentMergeMethod mergeMethod = + Acts::ComponentMergeMethod::eMaxWeight; IndexSourceLink::SurfaceAccessor m_slSurfaceAccessor; @@ -111,15 +112,19 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { weightCutoff, abortOnError, disableAllMaterialHandling}; - gsfOptions.stateReductionMethod = reductionMethod; + gsfOptions.componentMergeMethod = mergeMethod; gsfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>( &calibrator); gsfOptions.extensions.surfaceAccessor .connect<&IndexSourceLink::SurfaceAccessor::operator()>( &m_slSurfaceAccessor); - gsfOptions.extensions.mixtureReducer - .connect<&Acts::reduceMixtureWithKLDistance>(); + switch (reductionAlg) { + case MixtureReductionAlgorithm::KLDistance: { + gsfOptions.extensions.mixtureReducer + .connect<&Acts::reduceMixtureWithKLDistance>(); + } break; + } return gsfOptions; } @@ -167,9 +172,9 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - BetheHeitlerApprox betheHeitlerApprox, size_t maxComponents, - double weightCutoff, Acts::MixtureReductionMethod finalReductionMethod, - bool abortOnError, bool disableAllMaterialHandling, + BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, + double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger) { // Standard fitter MultiStepper stepper(magneticField, logger.cloneWithSuffix("Step")); @@ -202,9 +207,8 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::move(trackFitter), std::move(directTrackFitter), geo); fitterFunction->maxComponents = maxComponents; fitterFunction->weightCutoff = weightCutoff; - fitterFunction->abortOnError = abortOnError; - fitterFunction->disableAllMaterialHandling = disableAllMaterialHandling; - fitterFunction->reductionMethod = finalReductionMethod; + fitterFunction->mergeMethod = componentMergeMethod; + fitterFunction->reductionAlg = mixtureReductionAlgorithm; return fitterFunction; } diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 36b311a6781..8dbda185a9d 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -981,9 +981,8 @@ def addTruthTrackingGsf( gsfOptions = { "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(), "maxComponents": 12, - "abortOnError": False, - "disableAllMaterialHandling": False, - "finalReductionMethod": acts.examples.FinalReductionMethod.maxWeight, + "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight, + "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance, "weightCutoff": 1.0e-4, "level": customLogLevel(), } diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index cdba4a4cb0c..4fc7ac46f44 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -10,7 +10,7 @@ #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" #include "Acts/Plugins/Python/Utilities.hpp" #include "Acts/TrackFitting/BetheHeitlerApprox.hpp" -#include "Acts/Utilities/GaussianMixtureReduction.hpp" +#include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/MeasurementCalibration.hpp" @@ -97,9 +97,13 @@ void addTrackFitting(Context& ctx) { }, py::arg("path")); - py::enum_(mex, "FinalReductionMethod") - .value("mean", Acts::MixtureReductionMethod::eMean) - .value("maxWeight", Acts::MixtureReductionMethod::eMaxWeight); + py::enum_(mex, "ComponentMergeMethod") + .value("mean", Acts::ComponentMergeMethod::eMean) + .value("maxWeight", Acts::ComponentMergeMethod::eMaxWeight); + + py::enum_( + mex, "MixtureReductionAlgorithm") + .value("KLDistance", MixtureReductionAlgorithm::KLDistance); py::class_(mex, "AtlasBetheHeitlerApprox") .def_static("loadFromFiles", @@ -107,26 +111,24 @@ void addTrackFitting(Context& ctx) { py::arg("lowParametersPath"), py::arg("lowParametersPath")) .def_static("makeDefault", []() { return Acts::makeDefaultBetheHeitlerApprox(); }); - mex.def( "makeGsfFitterFunction", [](std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - BetheHeitlerApprox betheHeitlerApprox, size_t maxComponents, - double weightCutoff, - Acts::MixtureReductionMethod finalReductionMethod, bool abortOnError, - bool disableAllMaterialHandling, Logging::Level level) { + BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, + double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + ActsExamples::MixtureReductionAlgorithm mixtureReductionAlgorithm, + Logging::Level level) { return ActsExamples::makeGsfFitterFunction( trackingGeometry, magneticField, betheHeitlerApprox, - maxComponents, weightCutoff, finalReductionMethod, abortOnError, - disableAllMaterialHandling, + maxComponents, weightCutoff, componentMergeMethod, + mixtureReductionAlgorithm, *Acts::getDefaultLogger("GSFFunc", level)); }, py::arg("trackingGeometry"), py::arg("magneticField"), py::arg("betheHeitlerApprox"), py::arg("maxComponents"), - py::arg("weightCutoff"), py::arg("finalReductionMethod"), - py::arg("abortOnError"), py::arg("disableAllMaterialHandling"), - py::arg("level")); + py::arg("weightCutoff"), py::arg("componentMergeMethod"), + py::arg("mixtureReductionAlgorithm"), py::arg("level")); mex.def( "makeGlobalChiSquareFitterFunction", diff --git a/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp b/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp index 12b2c9cd746..ddac64213f4 100644 --- a/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp @@ -21,7 +21,7 @@ #include "Acts/Surfaces/PlaneSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Surfaces/SurfaceBounds.hpp" -#include "Acts/Utilities/GaussianMixtureReduction.hpp" +#include "Acts/TrackFitting/detail/GsfComponentMerging.hpp" #include "Acts/Utilities/Identity.hpp" #include "Acts/Utilities/Intersection.hpp" #include "Acts/Utilities/Result.hpp" From 1fa2f866c01170bb734b32b38586deed9c26ebe4 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 13 Nov 2023 21:47:44 +0100 Subject: [PATCH 02/16] update docs --- docs/core/reconstruction/track_fitting.md | 25 +++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/docs/core/reconstruction/track_fitting.md b/docs/core/reconstruction/track_fitting.md index 9520d20297a..59a627c993f 100644 --- a/docs/core/reconstruction/track_fitting.md +++ b/docs/core/reconstruction/track_fitting.md @@ -67,12 +67,6 @@ Simplified overview of the GSF algorithm. ### The Multi-Stepper To implement the GSF, a special stepper is needed, that can handle a multi-component state internally: The {class}`Acts::MultiEigenStepperLoop`, which is based on the {class}`Acts::EigenStepper` and thus shares a lot of code with it. It interfaces to the navigation as one aggregate state to limit the navigation overhead, but internally processes a multi-component state. How this aggregation is performed can be configured via a template parameter, by default weighted average is used ({struct}`Acts::WeightedComponentReducerLoop`). -At the end of the fit the multi-component state must be reduced to a single set of parameters with a corresponding covariance matrix. This is supported by the {class}`Acts::MultiEigenStepperLoop` in two different ways currently: The *mean* method computes the mean and the covariance matrix of the multi-component state, whereas the *maximum weight* method just returns the component with the maximum weight. This can be configured in the constructor of the {class}`Acts::MultiEigenStepperLoop` with the {enum}`Acts::MixtureReductionMethod`. In the future there is planned to add a *mode* finding method as well. - -:::{note} -In practice it turned out that the *maximum weight* method leads to better results so far. -::: - Even though the multi-stepper interface exposes only one aggregate state and thus is compatible with most standard tools, there is a special aborter is required to stop the navigation when the surface is reached, the {struct}`Acts::MultiStepperSurfaceReached`. It checks if all components have reached the target surface already and updates their state accordingly. Optionally, it also can stop the propagation when the aggregate state reaches the surface. @@ -87,9 +81,24 @@ outline: --- ``` -The fit can be customized with several options, e.g., the maximum number of components. All options can be found in the {struct}`Acts::GsfOptions`. +The fit can be customized with several options. Important ones are: +* *maximum components*: How many components at maximum should be kept. +* *weight cut*: When to drop components. +* *component merging*: How a multi-component state is reduced to a single set of parameters and covariance. The method can be chosen with the enum {enum}`Acts::ComponentMergeMethod`. Two methods are supported currently: + * The *mean* computes the mean and the covariance of the mean. + * *max weight* takes the parameters of component with the maximum weight and computes the variance around these. This is a cheap approximation of the mode, which is not implemented currently. + +:::{note} +A good starting configuration is to use 12 components, the *max weight* merging and the *KL distance* reduction. +::: + +All options can be found in the {struct}`Acts::GsfOptions`: -To simplify integration, the GSF returns an {struct}`Acts::KalmanFitterResult` object, the same as the {class}`Acts::KalmanFitter`. This allows to use the same analysis tools for both fitters. +```{doxygenstruct} Acts::GsfOptions +--- +outline: +--- +``` If the GSF finds the column with the string identifier *"gsf-final-multi-component-state"* (defined in `Acts::GsfConstants::kFinalMultiComponentStateColumn`) in the track container, it adds the final multi-component state to the track as a `std::optional>` object. From 22ec41dcccbdb635a8ba553bd546f7dc5dd23d70 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 13 Nov 2023 22:26:19 +0100 Subject: [PATCH 03/16] fix GSF early stopping --- .../Acts/TrackFitting/GaussianSumFitter.hpp | 18 ++++++++++++++++-- .../Acts/TrackFitting/detail/GsfActor.hpp | 3 ++- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index 030331650f8..6e17181dfa3 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -76,6 +76,20 @@ struct GaussianSumFitter { /// The actor type using GsfActor = detail::GsfActor; + /// This allows to break the propagation by setting the navigationBreak + /// TODO refactor once we can do this more elegantly + struct NavigationBreakAborter { + NavigationBreakAborter() = default; + + template + bool operator()(propagator_state_t& state, const stepper_t& /*stepper*/, + const navigator_t& navigator, + const Logger& /*logger*/) const { + return navigator.navigationBreak(state.navigation); + } + }; + /// @brief The fit function for the Direct navigator template class holder_t> @@ -92,7 +106,7 @@ struct GaussianSumFitter { // Initialize the forward propagation with the DirectNavigator auto fwdPropInitializer = [&sSequence, this](const auto& opts) { using Actors = ActionList; - using Aborters = AbortList<>; + using Aborters = AbortList; PropagatorOptions propOptions(opts.geoContext, opts.magFieldContext); @@ -147,7 +161,7 @@ struct GaussianSumFitter { // Initialize the forward propagation with the DirectNavigator auto fwdPropInitializer = [this](const auto& opts) { using Actors = ActionList; - using Aborters = AbortList; + using Aborters = AbortList; PropagatorOptions propOptions(opts.geoContext, opts.magFieldContext); diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 3800d3d63be..ea519b51398 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -317,7 +317,8 @@ struct GsfActor { // Break the navigation if we found all measurements if (m_cfg.numberMeasurements && result.measurementStates == m_cfg.numberMeasurements) { - navigator.targetReached(state.navigation, true); + ACTS_VERBOSE("Stop navigation because all measurements are found"); + navigator.navigationBreak(state.navigation, true); } } From da319cb1afefd5d3a9034636c8879c093126073f Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 13 Nov 2023 22:54:03 +0100 Subject: [PATCH 04/16] add momentum cutoff --- .../Acts/TrackFitting/GaussianSumFitter.hpp | 7 +- Core/include/Acts/TrackFitting/GsfOptions.hpp | 2 + .../Acts/TrackFitting/detail/GsfActor.hpp | 157 ++++++++++-------- .../TrackFitting/src/GsfFitterFunction.cpp | 3 + docs/core/reconstruction/track_fitting.md | 3 +- 5 files changed, 97 insertions(+), 75 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index 6e17181dfa3..28ba254be38 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -384,7 +384,12 @@ struct GaussianSumFitter { r.measurementStates++; r.processedStates++; - const auto& params = *fwdGsfResult.lastMeasurementState; + assert(!fwdGsfResult.lastMeasurementComponents.empty()); + assert(fwdGsfResult.lastMeasurementSurface != nullptr); + MultiComponentBoundTrackParameters params( + fwdGsfResult.lastMeasurementSurface->getSharedPtr(), + fwdGsfResult.lastMeasurementComponents, + sParameters.particleHypothesis()); return m_propagator.template propagate, decltype(bwdPropOptions), diff --git a/Core/include/Acts/TrackFitting/GsfOptions.hpp b/Core/include/Acts/TrackFitting/GsfOptions.hpp index 0cf35fb11e6..210e26c546e 100644 --- a/Core/include/Acts/TrackFitting/GsfOptions.hpp +++ b/Core/include/Acts/TrackFitting/GsfOptions.hpp @@ -101,6 +101,8 @@ struct GsfOptions { double weightCutoff = 1.e-4; + double momentumCutoff = 500_MeV; + bool abortOnError = false; bool disableAllMaterialHandling = false; diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index ea519b51398..ed4098ca8ea 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -46,7 +46,11 @@ struct GsfResult { /// The last multi-component measurement state. Used to initialize the /// backward pass. - std::optional lastMeasurementState; + std::vector> + lastMeasurementComponents; + + /// The last measurement surface. Used to initialize the backward pass. + const Acts::Surface* lastMeasurementSurface = nullptr; /// Some counting std::size_t measurementStates = 0; @@ -94,6 +98,9 @@ struct GsfActor { /// When to discard components double weightCutoff = 1.0e-4; + /// When to discard components + double momentumCutoff = 500_MeV; + /// When this option is enabled, material information on all surfaces is /// ignored. This disables the component convolution as well as the handling /// of energy. This may be useful for debugging. @@ -132,6 +139,8 @@ struct GsfActor { std::map weights; }; + using FiltProjector = MultiTrajectoryProjector; + /// @brief GSF actor operation /// /// @tparam propagator_state_t is the type of Propagator state @@ -260,7 +269,8 @@ struct GsfActor { return; } - updateStepper(state, stepper, tmpStates); + FiltProjector proj{tmpStates.traj, tmpStates.weights}; + updateStepper(state, stepper, navigator, tmpStates.tips, proj); } // We have material, we thus need a component cache since we will // convolute the components and later reduce them again before updating @@ -305,7 +315,8 @@ struct GsfActor { removeLowWeightComponents(componentCache); - updateStepper(state, stepper, navigator, componentCache); + auto proj = [](const auto& a) -> decltype(a) { return a; }; + updateStepper(state, stepper, navigator, componentCache, proj); } // If we have only done preUpdate before, now do postUpdate @@ -405,6 +416,13 @@ struct GsfActor { assert(p_prev + delta_p > 0. && "new momentum must be > 0"); new_pars[eBoundQOverP] = old_bound.charge() / (p_prev + delta_p); + const auto p_new = state.stepping.particleHypothesis.extractMomentum( + new_pars[eBoundQOverP]); + if (p_new < m_cfg.momentumCutoff) { + ACTS_VERBOSE("Skip new component with p=" << p_new << " GeV"); + continue; + } + // compute inverse variance of p from mixture and update covariance auto new_cov = old_bound.covariance().value(); @@ -451,50 +469,19 @@ struct GsfActor { } } - /// Function that updates the stepper from the MultiTrajectory - template - void updateStepper(propagator_state_t& state, const stepper_t& stepper, - const TemporaryStates& tmpStates) const { - auto cmps = stepper.componentIterable(state.stepping); - - for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) { - // we set ignored components to missed, so we can remove them after - // the loop - if (tmpStates.weights.at(idx) < m_cfg.weightCutoff) { - cmp.status() = Intersection3D::Status::missed; - continue; - } - - auto proxy = tmpStates.traj.getTrackState(idx); - - cmp.pars() = - MultiTrajectoryHelpers::freeFiltered(state.options.geoContext, proxy); - cmp.cov() = proxy.filteredCovariance(); - cmp.weight() = tmpStates.weights.at(idx); - } - - stepper.removeMissedComponents(state.stepping); - - // TODO we have two normalization passes here now, this can probably be - // optimized - detail::normalizeWeights(cmps, - [&](auto cmp) -> double& { return cmp.weight(); }); - } - /// Function that updates the stepper from the ComponentCache template + typename navigator_t, typename range_t, typename proj_t> void updateStepper(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& navigator, - const std::vector& componentCache) const { + const navigator_t& navigator, const range_t& range, + const proj_t& proj) const { const auto& surface = *navigator.currentSurface(state.navigation); // Clear components before adding new ones stepper.clearComponents(state.stepping); - // Finally loop over components - for (const auto& [weight, pars, cov] : componentCache) { - // Add the component to the stepper + for (const auto& cmp : range) { + const auto& [weight, pars, cov] = proj(cmp); BoundTrackParameters bound(surface.getSharedPtr(), pars, cov, stepper.particleHypothesis(state.stepping)); @@ -505,12 +492,12 @@ struct GsfActor { continue; } - auto& cmp = *res; - cmp.jacToGlobal() = surface.boundToFreeJacobian(state.geoContext, pars); - cmp.pathAccumulated() = state.stepping.pathAccumulated; - cmp.jacobian() = Acts::BoundMatrix::Identity(); - cmp.derivative() = Acts::FreeVector::Zero(); - cmp.jacTransport() = Acts::FreeMatrix::Identity(); + auto& proxy = *res; + proxy.jacToGlobal() = surface.boundToFreeJacobian(state.geoContext, pars); + proxy.pathAccumulated() = state.stepping.pathAccumulated; + proxy.jacobian() = Acts::BoundMatrix::Identity(); + proxy.derivative() = Acts::FreeVector::Zero(); + proxy.jacTransport() = Acts::FreeMatrix::Identity(); } } @@ -524,11 +511,14 @@ struct GsfActor { const SourceLink& source_link) const { const auto& surface = *navigator.currentSurface(state.navigation); + // This allows to easily project the to weight, filtered pars, filtered cov + const FiltProjector proj{tmpStates.traj, tmpStates.weights}; + // Boolean flag, to distinguish measurement and outlier states. This flag // is only modified by the valid-measurement-branch, so only if there // isn't any valid measurement state, the flag stays false and the state // is thus counted as an outlier - bool is_valid_measurement = false; + bool isValidMeasurement = false; auto cmps = stepper.componentIterable(state.stepping); for (auto cmp : cmps) { @@ -546,54 +536,73 @@ struct GsfActor { const auto& trackStateProxy = *trackStateProxyRes; + const auto p = state.stepping.particleHypothesis.extractMomentum( + trackStateProxy.filtered()[eBoundQOverP]); + if (p < m_cfg.momentumCutoff) { + ACTS_VERBOSE("Discard component with momentum " + << p << " GeV after Kalman update"); + continue; + } + // If at least one component is no outlier, we consider the whole thing // as a measurementState if (trackStateProxy.typeFlags().test( Acts::TrackStateFlag::MeasurementFlag)) { - is_valid_measurement = true; + isValidMeasurement = true; } tmpStates.tips.push_back(trackStateProxy.index()); tmpStates.weights[tmpStates.tips.back()] = cmp.weight(); } - computePosteriorWeights(tmpStates.traj, tmpStates.tips, tmpStates.weights); + // compute the posterior weights + if (!tmpStates.tips.empty()) { + computePosteriorWeights(tmpStates.traj, tmpStates.tips, + tmpStates.weights); - detail::normalizeWeights(tmpStates.tips, [&](auto idx) -> double& { - return tmpStates.weights.at(idx); - }); + detail::normalizeWeights(tmpStates.tips, [&](auto idx) -> double& { + return tmpStates.weights.at(idx); + }); + } + + // Remove low weight components + auto newEnd = std::remove_if( + tmpStates.tips.begin(), tmpStates.tips.end(), + [&](auto t) { return tmpStates.weights[t] < m_cfg.weightCutoff; }); + + if (newEnd != tmpStates.tips.end()) { + tmpStates.tips.erase(newEnd, tmpStates.tips.end()); + detail::normalizeWeights(tmpStates.tips, [&](auto idx) -> double& { + return tmpStates.weights.at(idx); + }); + } + + // Break navigation if we have no states left + if (tmpStates.tips.empty()) { + ACTS_DEBUG("no components left after Kalman update, break navigation!"); + navigator.navigationBreak(state.navigation, true); + return Acts::Result::success(); + } // Do the statistics ++result.processedStates; // TODO should outlier states also be counted here? - if (is_valid_measurement) { + if (isValidMeasurement) { ++result.measurementStates; } - addCombinedState(result, tmpStates, surface); - result.lastMeasurementTip = result.currentTip; - - using FiltProjector = - MultiTrajectoryProjector; - FiltProjector proj{tmpStates.traj, tmpStates.weights}; - - std::vector> v; + updateMultiTrajectory(result, tmpStates, surface); - // TODO Check why can zero weights can occur + result.lastMeasurementTip = result.currentTip; + result.lastMeasurementSurface = &surface; + result.lastMeasurementComponents.clear(); for (const auto& idx : tmpStates.tips) { const auto [w, p, c] = proj(idx); - if (w > 0.0) { - v.push_back({w, p, c}); - } + assert(w > 0.0); + result.lastMeasurementComponents.push_back({w, p, c}); } - normalizeWeights(v, [](auto& c) -> double& { return std::get(c); }); - - result.lastMeasurementState = MultiComponentBoundTrackParameters( - surface.getSharedPtr(), std::move(v), - stepper.particleHypothesis(state.stepping)); - // Return success return Acts::Result::success(); } @@ -644,7 +653,7 @@ struct GsfActor { ++result.processedStates; - addCombinedState(result, tmpStates, surface); + updateMultiTrajectory(result, tmpStates, surface); return Result::success(); } @@ -690,8 +699,9 @@ struct GsfActor { } } - void addCombinedState(result_type& result, const TemporaryStates& tmpStates, - const Surface& surface) const { + void updateMultiTrajectory(result_type& result, + const TemporaryStates& tmpStates, + const Surface& surface) const { using PrtProjector = MultiTrajectoryProjector; using FltProjector = @@ -768,6 +778,7 @@ struct GsfActor { m_cfg.abortOnError = options.abortOnError; m_cfg.disableAllMaterialHandling = options.disableAllMaterialHandling; m_cfg.weightCutoff = options.weightCutoff; + m_cfg.momentumCutoff = options.momentumCutoff; m_cfg.mergeMethod = options.componentMergeMethod; m_cfg.calibrationContext = &options.calibrationContext.get(); } diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index b90b2e2377d..b200b19b0d9 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -54,6 +54,7 @@ class TrackingGeometry; } // namespace Acts using namespace ActsExamples; +using namespace Acts::UnitLiterals; namespace { @@ -101,6 +102,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { &Acts::GainMatrixUpdater::operator()>( &updater); + const double momentumCutoff = 500_MeV; Acts::GsfOptions gsfOptions{ options.geoContext, options.magFieldContext, @@ -110,6 +112,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { &(*options.referenceSurface), maxComponents, weightCutoff, + momentumCutoff, abortOnError, disableAllMaterialHandling}; gsfOptions.componentMergeMethod = mergeMethod; diff --git a/docs/core/reconstruction/track_fitting.md b/docs/core/reconstruction/track_fitting.md index 59a627c993f..f52935c3c23 100644 --- a/docs/core/reconstruction/track_fitting.md +++ b/docs/core/reconstruction/track_fitting.md @@ -83,7 +83,8 @@ outline: The fit can be customized with several options. Important ones are: * *maximum components*: How many components at maximum should be kept. -* *weight cut*: When to drop components. +* *weight cut*: When to drop components because of too little weight. +* *momentum cut*: When to drop components because of too low momentum. This can help stabilizing the fit, as low momenta tend to disturb the navigation. * *component merging*: How a multi-component state is reduced to a single set of parameters and covariance. The method can be chosen with the enum {enum}`Acts::ComponentMergeMethod`. Two methods are supported currently: * The *mean* computes the mean and the covariance of the mean. * *max weight* takes the parameters of component with the maximum weight and computes the variance around these. This is a cheap approximation of the mode, which is not implemented currently. From e38b7e45f65773f852f1ceb23742ec3fc71c4147 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Tue, 14 Nov 2023 10:12:57 +0100 Subject: [PATCH 05/16] update --- .../include/Acts/TrackFitting/detail/GsfActor.hpp | 15 ++++++--------- .../include/Acts/TrackFitting/detail/GsfUtils.hpp | 1 - .../TrackFitting/src/GsfFitterFunction.cpp | 2 +- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index ed4098ca8ea..347d678bb5b 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -300,11 +300,9 @@ struct GsfActor { result.nInvalidBetheHeitler); if (componentCache.empty()) { - ACTS_WARNING( - "No components left after applying energy loss. " - "Is the weight cutoff " - << m_cfg.weightCutoff << " too high?"); - ACTS_WARNING("Return to propagator without applying energy loss"); + ACTS_DEBUG( + "No components left after applying energy loss, stop propagation"); + navigator.navigationBreak(state.navigation, true); return; } @@ -416,10 +414,9 @@ struct GsfActor { assert(p_prev + delta_p > 0. && "new momentum must be > 0"); new_pars[eBoundQOverP] = old_bound.charge() / (p_prev + delta_p); - const auto p_new = state.stepping.particleHypothesis.extractMomentum( - new_pars[eBoundQOverP]); - if (p_new < m_cfg.momentumCutoff) { - ACTS_VERBOSE("Skip new component with p=" << p_new << " GeV"); + if (p_prev + delta_p < m_cfg.momentumCutoff) { + ACTS_VERBOSE("Skip new component with p=" << p_prev + delta_p + << " GeV"); continue; } diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp index b52dc2d0c42..6bf15c5b0e5 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp @@ -109,7 +109,6 @@ class ScopedGsfInfoPrinterAndChecker { assert(allFinite && "weights not finite at the start"); assert(allNormalized && "not normalized at the start"); } else { - assert(!zeroComponents && "no cmps at the end"); assert(allFinite && "weights not finite at the end"); assert(allNormalized && "not normalized at the end"); } diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index b200b19b0d9..847f4976699 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -79,6 +79,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { size_t maxComponents = 0; double weightCutoff = 0; + const double momentumCutoff = 500_MeV; bool abortOnError = false; bool disableAllMaterialHandling = false; MixtureReductionAlgorithm reductionAlg = @@ -102,7 +103,6 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { &Acts::GainMatrixUpdater::operator()>( &updater); - const double momentumCutoff = 500_MeV; Acts::GsfOptions gsfOptions{ options.geoContext, options.magFieldContext, From d6346fbf49a859e096422b3e4e9ab22f1676d557 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Tue, 14 Nov 2023 14:26:29 +0100 Subject: [PATCH 06/16] Apply suggestions from code review --- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 347d678bb5b..baa925435be 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -508,7 +508,7 @@ struct GsfActor { const SourceLink& source_link) const { const auto& surface = *navigator.currentSurface(state.navigation); - // This allows to easily project the to weight, filtered pars, filtered cov + // This allows to easily project the state to weight, filtered pars, filtered cov const FiltProjector proj{tmpStates.traj, tmpStates.weights}; // Boolean flag, to distinguish measurement and outlier states. This flag From cc81762d44d2725212a110e5fa05515eafcfe8a0 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Tue, 14 Nov 2023 16:13:28 +0100 Subject: [PATCH 07/16] format --- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index baa925435be..20264466ba3 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -508,7 +508,8 @@ struct GsfActor { const SourceLink& source_link) const { const auto& surface = *navigator.currentSurface(state.navigation); - // This allows to easily project the state to weight, filtered pars, filtered cov + // This allows to easily project the state to weight, filtered pars, + // filtered cov const FiltProjector proj{tmpStates.traj, tmpStates.weights}; // Boolean flag, to distinguish measurement and outlier states. This flag From c41bdd3feba4762750483aa5540899f62e1950cd Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Tue, 14 Nov 2023 16:48:31 +0100 Subject: [PATCH 08/16] do it --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 78 +++++++++++-------- Core/src/TrackFitting/BetheHeitlerApprox.cpp | 2 +- Examples/Python/src/TrackFitting.cpp | 3 +- 3 files changed, 47 insertions(+), 36 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 66579bb9905..e7741222bd7 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -98,6 +98,9 @@ struct BetheHeitlerApproxSingleCmp { /// mixture. To enable an approximation for continuous input variables, the /// weights, means and variances are internally parametrized as a Nth order /// polynomial. +/// @todo This class is rather unflexible: It forces two data representations, +/// making it a bit awkward to add a single parameterization. It would be good +/// to generalize this at some point. template class AtlasBetheHeitlerApprox { static_assert(NComponents > 0); @@ -112,32 +115,37 @@ class AtlasBetheHeitlerApprox { using Data = std::array; - constexpr static double noChangeLimit = 0.0001; - constexpr static double singleGaussianLimit = 0.002; - constexpr static double lowerLimit = 0.10; - constexpr static double higherLimit = 0.20; - private: - Data m_low_data; - Data m_high_data; - bool m_low_transform; - bool m_high_transform; + Data m_lowData; + Data m_highData; + bool m_lowTransform; + bool m_highTransform; + + constexpr static double m_noChangeLimit = 0.0001; + constexpr static double m_singleGaussianLimit = 0.002; + double m_lowLimit = 0.10; + double m_highLimit = 0.20; public: - /// Construct the Bethe-Heitler approximation description. Additional to the - /// coefficients of the polynomials, the information whether these values need - /// to be transformed beforehand must be given (see ATLAS code). + /// Construct the Bethe-Heitler approximation description with two + /// parameterizations, one for lower ranges, one for higher ranges. + /// Is it assumed that the lower limit of the high-x/x0 data is equal + /// to the upper limit of the low-x/x0 data. /// - /// @param low_data data for the lower x/x0 range - /// @param high_data data for the higher x/x0 range - /// @param low_transform whether the low data need to be transformed - /// @param high_transform whether the high data need to be transformed - constexpr AtlasBetheHeitlerApprox(const Data &low_data, const Data &high_data, - bool low_transform, bool high_transform) - : m_low_data(low_data), - m_high_data(high_data), - m_low_transform(low_transform), - m_high_transform(high_transform) {} + /// @param lowData data for the lower x/x0 range + /// @param highData data for the higher x/x0 range + /// @param lowTransform whether the low data need to be transformed + /// @param highTransform whether the high data need to be transformed + /// @param lowLimit the upper limit for the low data + /// @param highLimit the upper limit for the high data + constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData, + bool lowTransform, bool highTransform, double lowLimit=0.1, double highLimit=0.2) + : m_lowData(lowData), + m_highData(highData), + m_lowTransform(lowTransform), + m_highTransform(highTransform), + m_lowLimit(lowLimit), + m_highLimit(highLimit) {} /// Returns the number of components the returned mixture will have constexpr auto numComponents() const { return NComponents; } @@ -145,7 +153,7 @@ class AtlasBetheHeitlerApprox { /// Checks if an input is valid for the parameterization /// /// @param x pathlength in terms of the radiation length - constexpr bool validXOverX0(ActsScalar x) const { return x < higherLimit; } + constexpr bool validXOverX0(ActsScalar x) const { return x < m_highLimit; } /// Generates the mixture from the polynomials and reweights them, so /// that the sum of all weights is 1 @@ -194,7 +202,7 @@ class AtlasBetheHeitlerApprox { }; // Return no change - if (x < noChangeLimit) { + if (x < m_noChangeLimit) { Array ret(1); ret[0].weight = 1.0; @@ -204,19 +212,19 @@ class AtlasBetheHeitlerApprox { return ret; } // Return single gaussian approximation - if (x < singleGaussianLimit) { + if (x < m_singleGaussianLimit) { Array ret(1); ret[0] = BetheHeitlerApproxSingleCmp::mixture(x)[0]; return ret; } // Return a component representation for lower x0 - if (x < lowerLimit) { - return make_mixture(m_low_data, x, m_low_transform); + if (x < m_lowLimit) { + return make_mixture(m_lowData, x, m_lowTransform); } // Return a component representation for higher x0 // Cap the x because beyond the parameterization goes wild - const auto high_x = std::min(higherLimit, x); - return make_mixture(m_high_data, high_x, m_high_transform); + const auto high_x = std::min(m_highLimit, x); + return make_mixture(m_highData, high_x, m_highTransform); } /// Loads a parameterization from a file according to the Atlas file @@ -226,8 +234,10 @@ class AtlasBetheHeitlerApprox { /// the parameterization for low x/x0 /// @param high_parameters_path Path to the foo.par file that stores /// the parameterization for high x/x0 + /// @param lowLimit the upper limit for the low x/x0-data + /// @param highLimit the upper limit for the high x/x0-data static auto loadFromFiles(const std::string &low_parameters_path, - const std::string &high_parameters_path) { + const std::string &high_parameters_path, double lowLimit = 0.1, double highLimit = 0.2) { auto read_file = [](const std::string &filepath) { std::ifstream file(filepath); @@ -267,11 +277,11 @@ class AtlasBetheHeitlerApprox { return std::make_tuple(data, transform_code); }; - const auto [low_data, low_transform] = read_file(low_parameters_path); - const auto [high_data, high_transform] = read_file(high_parameters_path); + const auto [lowData, lowTransform] = read_file(low_parameters_path); + const auto [highData, highTransform] = read_file(high_parameters_path); - return AtlasBetheHeitlerApprox(low_data, high_data, low_transform, - high_transform); + return AtlasBetheHeitlerApprox(lowData, highData, lowTransform, + highTransform, lowLimit, highLimit); } }; diff --git a/Core/src/TrackFitting/BetheHeitlerApprox.cpp b/Core/src/TrackFitting/BetheHeitlerApprox.cpp index 8b862007050..06d4b802505 100644 --- a/Core/src/TrackFitting/BetheHeitlerApprox.cpp +++ b/Core/src/TrackFitting/BetheHeitlerApprox.cpp @@ -52,5 +52,5 @@ Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() { // clang-format on return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data, - cdf_cmps6_order5_data, true, true); + cdf_cmps6_order5_data, true, true, 0.2, 0.2); } diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index cdba4a4cb0c..a2f69a7d578 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -104,7 +104,8 @@ void addTrackFitting(Context& ctx) { py::class_(mex, "AtlasBetheHeitlerApprox") .def_static("loadFromFiles", &ActsExamples::BetheHeitlerApprox::loadFromFiles, - py::arg("lowParametersPath"), py::arg("lowParametersPath")) + py::arg("lowParametersPath"), py::arg("highParametersPath"), + py::arg("lowLimit") = 0.1, py::arg("highLimit") = 0.2) .def_static("makeDefault", []() { return Acts::makeDefaultBetheHeitlerApprox(); }); From 92946c96459b4f10013329cb84b189d94887d9b3 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 16 Nov 2023 10:17:46 +0100 Subject: [PATCH 09/16] improve comment --- .../Acts/TrackFitting/GaussianSumFitter.hpp | 14 +++++++---- .../Acts/TrackFitting/detail/GsfActor.hpp | 23 +++++++++++-------- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index 6e17181dfa3..0a5a0163c58 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -323,6 +323,7 @@ struct GaussianSumFitter { ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates); std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler; + double maxPathXOverX0 = fwdResult.maxPathXOverX0; ////////////////// // Backward pass @@ -409,13 +410,16 @@ struct GaussianSumFitter { } nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler; + maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0); if (nInvalidBetheHeitler > 0) { - ACTS_WARNING("Encountered " - << nInvalidBetheHeitler - << " cases where the material thickness exceeds the range " - "of the Bethe-Heitler-Approximation. Enable DEBUG output " - "for more information."); + ACTS_WARNING("Encountered " << nInvalidBetheHeitler + << " cases where x/X0 exceeds the range " + "of the Bethe-Heitler-Approximation. The " + "maximum x/X0 encountered was " + << maxPathXOverX0 + << ". Enable DEBUG output " + "for more information."); } //////////////////////////////////// diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index ea519b51398..055bb342b19 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -56,7 +56,9 @@ struct GsfResult { std::vector visitedSurfaces; std::vector surfacesVisitedBwdAgain; + /// Statistics about encounterings of invalid bethe heitler approximations std::size_t nInvalidBetheHeitler = 0; + double maxPathXOverX0 = 0.0; // Propagate potential errors to the outside Result result{Result::success()}; @@ -287,7 +289,7 @@ struct GsfActor { componentCache.clear(); convoluteComponents(state, stepper, navigator, tmpStates, componentCache, - result.nInvalidBetheHeitler); + result); if (componentCache.empty()) { ACTS_WARNING( @@ -328,7 +330,7 @@ struct GsfActor { const navigator_t& navigator, const TemporaryStates& tmpStates, std::vector& componentCache, - std::size_t& nInvalidBetheHeitler) const { + result_type &result) const { auto cmps = stepper.componentIterable(state.stepping); for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) { auto proxy = tmpStates.traj.getTrackState(idx); @@ -338,7 +340,7 @@ struct GsfActor { stepper.particleHypothesis(state.stepping)); applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx), - componentCache, nInvalidBetheHeitler); + componentCache, result); } } @@ -348,7 +350,7 @@ struct GsfActor { const BoundTrackParameters& old_bound, const double old_weight, std::vector& componentCaches, - std::size_t& nInvalidBetheHeitler) const { + result_type& result) const { const auto& surface = *navigator.currentSurface(state.navigation); const auto p_prev = old_bound.absoluteMomentum(); @@ -357,22 +359,25 @@ struct GsfActor { old_bound.position(state.stepping.geoContext), state.options.direction, MaterialUpdateStage::FullUpdate); - auto pathCorrection = surface.pathCorrection( + const auto pathCorrection = surface.pathCorrection( state.stepping.geoContext, old_bound.position(state.stepping.geoContext), old_bound.direction()); slab.scaleThickness(pathCorrection); + + const double pathXOverX0 = slab.thicknessInX0(); // Emit a warning if the approximation is not valid for this x/x0 - if (!m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0())) { - ++nInvalidBetheHeitler; + if (!m_cfg.bethe_heitler_approx->validXOverX0(pathXOverX0)) { + ++result.nInvalidBetheHeitler; + result.maxPathXOverX0 = std::max(result.maxPathXOverX0, pathXOverX0); ACTS_DEBUG( "Bethe-Heitler approximation encountered invalid value for x/x0=" - << slab.thicknessInX0() << " at surface " << surface.geometryId()); + << pathXOverX0 << " at surface " << surface.geometryId()); } // Get the mixture const auto mixture = - m_cfg.bethe_heitler_approx->mixture(slab.thicknessInX0()); + m_cfg.bethe_heitler_approx->mixture(pathXOverX0); // Create all possible new components for (const auto& gaussian : mixture) { From d3b6936fd5a1cd94755943dc6a3df5d9d0d378c8 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 16 Nov 2023 10:20:05 +0100 Subject: [PATCH 10/16] format --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 11 +++++++---- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 7 +++---- Core/src/TrackFitting/BetheHeitlerApprox.cpp | 4 ++-- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index e7741222bd7..9e2a0213424 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -98,7 +98,7 @@ struct BetheHeitlerApproxSingleCmp { /// mixture. To enable an approximation for continuous input variables, the /// weights, means and variances are internally parametrized as a Nth order /// polynomial. -/// @todo This class is rather unflexible: It forces two data representations, +/// @todo This class is rather unflexible: It forces two data representations, /// making it a bit awkward to add a single parameterization. It would be good /// to generalize this at some point. template @@ -120,7 +120,7 @@ class AtlasBetheHeitlerApprox { Data m_highData; bool m_lowTransform; bool m_highTransform; - + constexpr static double m_noChangeLimit = 0.0001; constexpr static double m_singleGaussianLimit = 0.002; double m_lowLimit = 0.10; @@ -139,7 +139,9 @@ class AtlasBetheHeitlerApprox { /// @param lowLimit the upper limit for the low data /// @param highLimit the upper limit for the high data constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData, - bool lowTransform, bool highTransform, double lowLimit=0.1, double highLimit=0.2) + bool lowTransform, bool highTransform, + double lowLimit = 0.1, + double highLimit = 0.2) : m_lowData(lowData), m_highData(highData), m_lowTransform(lowTransform), @@ -237,7 +239,8 @@ class AtlasBetheHeitlerApprox { /// @param lowLimit the upper limit for the low x/x0-data /// @param highLimit the upper limit for the high x/x0-data static auto loadFromFiles(const std::string &low_parameters_path, - const std::string &high_parameters_path, double lowLimit = 0.1, double highLimit = 0.2) { + const std::string &high_parameters_path, + double lowLimit = 0.1, double highLimit = 0.2) { auto read_file = [](const std::string &filepath) { std::ifstream file(filepath); diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 055bb342b19..f348451b660 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -330,7 +330,7 @@ struct GsfActor { const navigator_t& navigator, const TemporaryStates& tmpStates, std::vector& componentCache, - result_type &result) const { + result_type& result) const { auto cmps = stepper.componentIterable(state.stepping); for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) { auto proxy = tmpStates.traj.getTrackState(idx); @@ -363,7 +363,7 @@ struct GsfActor { state.stepping.geoContext, old_bound.position(state.stepping.geoContext), old_bound.direction()); slab.scaleThickness(pathCorrection); - + const double pathXOverX0 = slab.thicknessInX0(); // Emit a warning if the approximation is not valid for this x/x0 @@ -376,8 +376,7 @@ struct GsfActor { } // Get the mixture - const auto mixture = - m_cfg.bethe_heitler_approx->mixture(pathXOverX0); + const auto mixture = m_cfg.bethe_heitler_approx->mixture(pathXOverX0); // Create all possible new components for (const auto& gaussian : mixture) { diff --git a/Core/src/TrackFitting/BetheHeitlerApprox.cpp b/Core/src/TrackFitting/BetheHeitlerApprox.cpp index 06d4b802505..45eca1d0156 100644 --- a/Core/src/TrackFitting/BetheHeitlerApprox.cpp +++ b/Core/src/TrackFitting/BetheHeitlerApprox.cpp @@ -51,6 +51,6 @@ Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() { }}; // clang-format on - return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data, - cdf_cmps6_order5_data, true, true, 0.2, 0.2); + return AtlasBetheHeitlerApprox<6, 5>( + cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true, 0.2, 0.2); } From ea83c2cbaf461c77d347a023ff7b8d7a50965cd0 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 16 Nov 2023 10:29:26 +0100 Subject: [PATCH 11/16] fix --- Core/include/Acts/TrackFitting/GaussianSumFitter.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index eb3d7f088bd..adc918b6eac 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -322,7 +322,7 @@ struct GaussianSumFitter { ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates); std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler; - double maxPathXOverX0 = fwdResult.maxPathXOverX0; + double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0; ////////////////// // Backward pass From 89359613f55cac7d14b0f95f5a35a1458880081c Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 16 Nov 2023 10:40:38 +0100 Subject: [PATCH 12/16] update --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 9e2a0213424..b13f11c2d53 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -98,7 +98,7 @@ struct BetheHeitlerApproxSingleCmp { /// mixture. To enable an approximation for continuous input variables, the /// weights, means and variances are internally parametrized as a Nth order /// polynomial. -/// @todo This class is rather unflexible: It forces two data representations, +/// @todo This class is rather inflexible: It forces two data representations, /// making it a bit awkward to add a single parameterization. It would be good /// to generalize this at some point. template From a65868f07a252275c22ce2fc40d8877144caf564 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 16 Nov 2023 15:04:49 +0100 Subject: [PATCH 13/16] fix --- Core/include/Acts/TrackFitting/detail/GsfUtils.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp index 6bf15c5b0e5..8566e1c7e82 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp @@ -109,8 +109,8 @@ class ScopedGsfInfoPrinterAndChecker { assert(allFinite && "weights not finite at the start"); assert(allNormalized && "not normalized at the start"); } else { + assert((zeroComponents || allNormalized) && "not normalized at the end"); assert(allFinite && "weights not finite at the end"); - assert(allNormalized && "not normalized at the end"); } } From ee62aecad0d4497e5eea3068e123be0b788962d2 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Sat, 18 Nov 2023 11:57:58 +0100 Subject: [PATCH 14/16] update --- .../ActsExamples/TrackFitting/TrackFitterFunction.hpp | 2 +- .../Algorithms/TrackFitting/src/GsfFitterFunction.cpp | 5 +++-- Examples/Python/python/acts/examples/reconstruction.py | 1 + Examples/Python/src/TrackFitting.cpp | 10 ++++++---- 4 files changed, 11 insertions(+), 7 deletions(-) diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp index f272d217772..5d79fb3d5b4 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp @@ -93,7 +93,7 @@ std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, - double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + double weightCutoff, double momentumCutoff, Acts::ComponentMergeMethod componentMergeMethod, MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger); diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 6e68677d12e..7f732351690 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -79,7 +79,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction { std::size_t maxComponents = 0; double weightCutoff = 0; - const double momentumCutoff = 500_MeV; + double momentumCutoff = 0; bool abortOnError = false; bool disableAllMaterialHandling = false; MixtureReductionAlgorithm reductionAlg = @@ -180,7 +180,7 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, - double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + double weightCutoff, double momentumCutoff, Acts::ComponentMergeMethod componentMergeMethod, MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger) { // Standard fitter @@ -214,6 +214,7 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::move(trackFitter), std::move(directTrackFitter), geo); fitterFunction->maxComponents = maxComponents; fitterFunction->weightCutoff = weightCutoff; + fitterFunction->momentumCutoff = momentumCutoff; fitterFunction->mergeMethod = componentMergeMethod; fitterFunction->reductionAlg = mixtureReductionAlgorithm; diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 7ab4c4a072e..e68c178831a 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -982,6 +982,7 @@ def addTruthTrackingGsf( "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight, "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance, "weightCutoff": 1.0e-4, + "momentumCutoff": 250 * u.MeV, "level": customLogLevel(), } diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 97091012c3d..202ebfd94b1 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -118,19 +118,21 @@ void addTrackFitting(Context& ctx) { [](std::shared_ptr trackingGeometry, std::shared_ptr magneticField, BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, - double weightCutoff, Acts::ComponentMergeMethod componentMergeMethod, + double weightCutoff, double momentumCutoff, + Acts::ComponentMergeMethod componentMergeMethod, ActsExamples::MixtureReductionAlgorithm mixtureReductionAlgorithm, Logging::Level level) { return ActsExamples::makeGsfFitterFunction( trackingGeometry, magneticField, betheHeitlerApprox, - maxComponents, weightCutoff, componentMergeMethod, + maxComponents, weightCutoff, momentumCutoff, componentMergeMethod, mixtureReductionAlgorithm, *Acts::getDefaultLogger("GSFFunc", level)); }, py::arg("trackingGeometry"), py::arg("magneticField"), py::arg("betheHeitlerApprox"), py::arg("maxComponents"), - py::arg("weightCutoff"), py::arg("componentMergeMethod"), - py::arg("mixtureReductionAlgorithm"), py::arg("level")); + py::arg("weightCutoff"), py::arg("momentumCutoff"), + py::arg("componentMergeMethod"), py::arg("mixtureReductionAlgorithm"), + py::arg("level")); mex.def( "makeGlobalChiSquareFitterFunction", From b7aa0fe726c61c8e9b080f34c61778943e3f4418 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Sat, 18 Nov 2023 12:00:07 +0100 Subject: [PATCH 15/16] update --- .../include/ActsExamples/TrackFitting/TrackFitterFunction.hpp | 3 ++- Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp index 5d79fb3d5b4..cd7a061b450 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp @@ -93,7 +93,8 @@ std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, - double weightCutoff, double momentumCutoff, Acts::ComponentMergeMethod componentMergeMethod, + double weightCutoff, double momentumCutoff, + Acts::ComponentMergeMethod componentMergeMethod, MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger); diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 7f732351690..0fb5b4ad94e 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -180,7 +180,8 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, - double weightCutoff, double momentumCutoff, Acts::ComponentMergeMethod componentMergeMethod, + double weightCutoff, double momentumCutoff, + Acts::ComponentMergeMethod componentMergeMethod, MixtureReductionAlgorithm mixtureReductionAlgorithm, const Acts::Logger& logger) { // Standard fitter From 213347c040b6b959883ac9886e678a61a14375e1 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Sat, 18 Nov 2023 19:53:55 +0100 Subject: [PATCH 16/16] update --- Core/include/Acts/TrackFitting/GsfOptions.hpp | 2 +- Examples/Python/python/acts/examples/reconstruction.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GsfOptions.hpp b/Core/include/Acts/TrackFitting/GsfOptions.hpp index 9ccdc6baa3f..fa55a4fff71 100644 --- a/Core/include/Acts/TrackFitting/GsfOptions.hpp +++ b/Core/include/Acts/TrackFitting/GsfOptions.hpp @@ -104,7 +104,7 @@ struct GsfOptions { double weightCutoff = 1.e-4; - double momentumCutoff = 500_MeV; + double momentumCutoff = 0; bool abortOnError = false; diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index e68c178831a..f76c7edd56a 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -982,7 +982,7 @@ def addTruthTrackingGsf( "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight, "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance, "weightCutoff": 1.0e-4, - "momentumCutoff": 250 * u.MeV, + "momentumCutoff": 100 * u.MeV, "level": customLogLevel(), }