diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index 201ab564bdb..1f3dc8fdd4e 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -508,7 +508,7 @@ struct GaussianSumFitter { const auto [finalPars, finalCov] = detail::Gsf::mergeGaussianMixture( params.components(), [](const auto& cmp) { - const auto& [weight_l, pars_l, opt_cov_l] = cmp; + auto&& [weight_l, pars_l, opt_cov_l] = cmp; return std::tie(weight_l, pars_l, *opt_cov_l); }, params.referenceSurface(), options.componentMergeMethod); diff --git a/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp index 59a7a11e3ad..d542fcd4a33 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfComponentMerging.hpp @@ -20,7 +20,9 @@ namespace Acts::detail::Gsf { -/// Reduce Gaussian mixture +/// Reduce Gaussian mixture into a single parameter vector and covariance, using +/// the specified method to reduce the mixture. +/// /// @param mixture The mixture iterable /// @param surface The surface, on which the bound state is /// @param method How to reduce the mixture @@ -29,7 +31,8 @@ std::tuple mergeGaussianMixture( std::span mixture, const Surface &surface, ComponentMergeMethod method); -/// Merge two GSF components into one +/// Merge two Gaussian mixture components into one +/// /// @param a First component /// @param b Second component /// @param surface The surface, on which the bound state is @@ -93,15 +96,83 @@ auto angleDescriptionSwitch(const Surface &surface, Callable &&callable) { } } +/// Combine multiple components into one representative parameter object. The +/// function takes a range and a projector to allow for arbitrary input to be +/// combined. +/// +/// @param cmps The component range to merge +/// @param projector A projector to extract the weight and parameters from the components +/// @param angleDesc The angle description object +/// @return parameters +template +BoundVector mergeGaussianMixtureMean(const component_range_t &cmps, + const projector_t &projector, + const angle_desc_t &angleDesc) { + // Early return in case of range with length 1 + if (cmps.size() == 1) { + const auto [weight_l, pars_l] = projector(cmps.front()); + return pars_l; + } + + // Do the (circular) mean with complex arithmetic. + // For normal values, just keep the real values. For angles, use the complex + // phase. Weighting then transparently happens by multiplying a real-valued + // weight. + using CVec = Eigen::Matrix, eBoundSize, 1>; + CVec cMean = CVec::Zero(); + double sumOfWeights = 0; + + for (const auto &cmp : cmps) { + const auto [weight_l, pars_l] = projector(cmp); + CVec cPars_l = pars_l; + + const auto setPolar = [&](const auto &desc) { + cPars_l[desc.idx] = std::polar(1.0, pars_l[desc.idx] / desc.constant); + }; + std::apply([&](auto... dsc) { (setPolar(dsc), ...); }, angleDesc); + + sumOfWeights += weight_l; + cMean += weight_l * cPars_l; + } + + cMean /= sumOfWeights; + + BoundVector mean = cMean.real(); + + const auto getArg = [&](const auto &desc) { + mean[desc.idx] = desc.constant * std::arg(cMean[desc.idx]); + }; + std::apply([&](auto... dsc) { (getArg(dsc), ...); }, angleDesc); + + return mean; +} + +/// Combine multiple components into one representative parameter object. The +/// function takes a range and a projector to allow for arbitrary input to be +/// combined. +/// +/// @param cmps The component range to merge +/// @param projector A projector to extract the weight and parameters from the components +/// @return parameters +template +BoundVector mergeGaussianMixtureMaxWeight(const component_range_t &cmps, + const projector_t &projector) { + const auto maxWeightIt = + std::ranges::max_element(cmps, {}, [&](const auto &cmp) { + const auto [weight_l, pars_l] = projector(cmp); + return weight_l; + }); + const auto [weight_max, pars_max] = projector(*maxWeightIt); + return pars_max; +} + /// Compute the covariance of a Gaussian mixture given the mean. The function -/// takes iterators to allow for arbitrary ranges to be combined. The dimension -/// of the vectors is inferred from the inputs. +/// takes a range and a projector to allow for arbitrary input to be combined. /// /// @param cmps The component range to merge -/// @param projector A projector to extract the weight, parameters and covariance -/// from the components +/// @param projector A projector to extract the weight and parameters from the components /// @param mean The precomputed mean of the mixture -/// @param sumOfWeights The precomputed sum of weights of the mixture /// @param angleDesc The angle description object /// @return covariance template mergeGaussianMixtureMeanCov( const component_range_t &cmps, const projector_t &projector, const angle_desc_t &angleDesc) { - // Early return in case of range with length 1 if (cmps.size() == 1) { - const auto &[beginWeight, beginPars, beginCov] = projector(cmps.front()); - return {beginPars / beginWeight, beginCov / beginWeight}; + const auto [weight_l, pars_l, cov_l] = projector(cmps.front()); + return {pars_l, cov_l}; } - // Do the (circular) mean with complex arithmetic. - // For normal values, just keep the real values. For angles, use the complex - // phase. Weighting then transparently happens by multiplying a real-valued - // weight. - using CVec = Eigen::Matrix, eBoundSize, 1>; - CVec cMean = CVec::Zero(); - double sumOfWeights{0.0}; - - for (const auto &cmp : cmps) { - const auto &[weight_l, pars_l, cov_l] = projector(cmp); - - CVec pars_l_c = pars_l; - - const auto setPolar = [&](const auto &desc) { - pars_l_c[desc.idx] = std::polar(1.0, pars_l[desc.idx] / desc.constant); - }; - std::apply([&](auto... dsc) { (setPolar(dsc), ...); }, angleDesc); - - sumOfWeights += weight_l; - cMean += weight_l * pars_l_c; - } - - cMean /= sumOfWeights; - - BoundVector mean = cMean.real(); - - const auto getArg = [&](const auto &desc) { - mean[desc.idx] = desc.constant * std::arg(cMean[desc.idx]); - }; - std::apply([&](auto... dsc) { (getArg(dsc), ...); }, angleDesc); + const BoundVector mean = mergeGaussianMixtureMean( + cmps, + [&](const auto &c) { + const auto [weight_l, pars_l, cov_l] = projector(c); + return std::tuple(weight_l, pars_l); + }, + angleDesc); // MARK: fpeMaskBegin(FLTUND, 1, #2347) const BoundMatrix cov = - mergeGaussianMixtureCov(cmps, projector, mean, sumOfWeights, angleDesc); + mergeGaussianMixtureCov(cmps, projector, mean, angleDesc); // MARK: fpeMaskEnd(FLTUND) return {mean, cov}; } -/// Reduce Gaussian mixture +/// Reduce Gaussian mixture into a single parameter vector, using the specified +/// method to reduce the mixture. +/// +/// @param cmps The component range to merge +/// @param projector A projector to extract the weight and parameters from the components +/// @param angleDesc The angle description object +/// @param method How to reduce the mixture +/// @return parameters +template +BoundVector mergeGaussianMixtureParams(const component_range_t &cmps, + const projector_t &projector, + const angle_desc_t &angleDesc, + ComponentMergeMethod method) { + if (method == ComponentMergeMethod::eMean) { + return mergeGaussianMixtureMean(cmps, projector, angleDesc); + } else if (method == ComponentMergeMethod::eMaxWeight) { + return mergeGaussianMixtureMaxWeight(cmps, projector); + } else { + throw std::logic_error("Invalid component merge method"); + } +} + +/// Reduce Gaussian mixture into a single parameter vector and covariance, using +/// the specified method to reduce the mixture. /// /// @param cmps The component range to merge /// @param projector A projector to extract the weight, parameters and covariance @@ -223,20 +293,37 @@ std::tuple mergeGaussianMixture( if (method == ComponentMergeMethod::eMean) { return {mean, cov}; } else if (method == ComponentMergeMethod::eMaxWeight) { - const auto maxWeightIt = - std::ranges::max_element(cmps, {}, [&](const auto &cmp) { - const auto &[weight_l, pars_l, cov_l] = projector(cmp); - return weight_l; + const BoundVector mode = + mergeGaussianMixtureMaxWeight(cmps, [&](const auto &c) { + const auto [weight_l, pars_l, cov_l] = projector(c); + return std::tuple(weight_l, pars_l); }); - const auto &[weight_l, pars_l, cov_l] = projector(*maxWeightIt); - - return {pars_l, cov}; + return {mode, cov}; } else { throw std::logic_error("Invalid component merge method"); } } -/// Reduce Gaussian mixture +/// Reduce Gaussian mixture into a single parameter vector, using the specified +/// method to reduce the mixture. +/// +/// @param cmps The component range to merge +/// @param projector A projector to extract the weight and parameters from the components +/// @param surface The surface, which the bound state is on +/// @param method How to reduce the mixture +/// @return parameters +template +BoundVector mergeGaussianMixtureParams(const component_range_t &cmps, + const projector_t &projector, + const Surface &surface, + ComponentMergeMethod method) { + return angleDescriptionSwitch(surface, [&](const auto &desc) { + return mergeGaussianMixtureParams(cmps, projector, desc, method); + }); +} + +/// Reduce Gaussian mixture into a single parameter vector and covariance, using +/// the specified method to reduce the mixture. /// /// @param cmps The component range to merge /// @param projector A projector to extract the weight, parameters and covariance diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp index 44135fc34e8..0b86035859e 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp @@ -226,14 +226,14 @@ struct MultiTrajectoryProjector { const auto proxy = mt.getTrackState(idx); switch (type) { case StatesType::ePredicted: - return std::make_tuple(weights.at(idx), proxy.predicted(), - proxy.predictedCovariance()); + return std::tuple(weights.at(idx), proxy.predicted(), + proxy.predictedCovariance()); case StatesType::eFiltered: - return std::make_tuple(weights.at(idx), proxy.filtered(), - proxy.filteredCovariance()); + return std::tuple(weights.at(idx), proxy.filtered(), + proxy.filteredCovariance()); case StatesType::eSmoothed: - return std::make_tuple(weights.at(idx), proxy.smoothed(), - proxy.smoothedCovariance()); + return std::tuple(weights.at(idx), proxy.smoothed(), + proxy.smoothedCovariance()); default: throw std::invalid_argument( "Incorrect StatesType, should be ePredicted" diff --git a/Core/src/TrackFitting/detail/GsfComponentMerging.cpp b/Core/src/TrackFitting/detail/GsfComponentMerging.cpp index 831caf3bbca..2b864be2414 100644 --- a/Core/src/TrackFitting/detail/GsfComponentMerging.cpp +++ b/Core/src/TrackFitting/detail/GsfComponentMerging.cpp @@ -15,7 +15,12 @@ namespace Acts { std::tuple detail::Gsf::mergeGaussianMixture( std::span mixture, const Surface &surface, ComponentMergeMethod method) { - return mergeGaussianMixture(mixture, std::identity{}, surface, method); + return mergeGaussianMixture( + mixture, + [](const GsfComponent &c) { + return std::tie(c.weight, c.boundPars, c.boundCov); + }, + surface, method); } GsfComponent detail::Gsf::mergeTwoComponents(const GsfComponent &a, diff --git a/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp b/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp index 41fb0b1ffa7..7e3d258bc2f 100644 --- a/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp @@ -301,7 +301,11 @@ BOOST_AUTO_TEST_CASE(test_with_data_circular) { using detail::Gsf::CyclicAngle; const auto d = std::tuple, CyclicAngle>{}; const auto [mean_test, boundCov_test] = detail::Gsf::mergeGaussianMixture( - cmps, std::identity{}, d, ComponentMergeMethod::eMean); + cmps, + [](const GsfComponent &c) { + return std::tie(c.weight, c.boundPars, c.boundCov); + }, + d, ComponentMergeMethod::eMean); BOOST_CHECK(std::abs(detail::difference_periodic( mean_data[0], mean_test[0], 2 * std::numbers::pi)) < 1.e-1);