diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
index 6f71521e4d2d2ad8a2f7e556881f3004f62fab58..4de6f11955e499e48a75f615358f1fc20acdadef 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GaussianSumFitter.h
@@ -11,7 +11,6 @@
#include "TrkGaussianSumFilter/IMultiStateExtrapolator.h"
#include "TrkGaussianSumFilterUtils/GsfMeasurementUpdator.h"
-#include "TrkGaussianSumFilterUtils/QuickCloseComponentsMultiStateMerger.h"
//
#include "TrkMultiComponentStateOnSurface/MultiComponentStateOnSurface.h"
//
@@ -119,7 +118,7 @@ private:
const ParticleHypothesis particleHypothesis = nonInteracting) const;
/** Gsf smoothe trajectory*/
- MultiTrajectory fit(
+ MultiTrajectory smootherFit(
const EventContext& ctx,
Trk::IMultiStateExtrapolator::Cache&,
const MultiTrajectory&,
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
index 617b80bdb6620f30a2f248cd3603e71072c2b9ef..907cbade8d93e654275cf8907589b7130b91e00a 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GaussianSumFitter.cxx
@@ -60,95 +60,6 @@ buildFitQuality(const Trk::SmoothedTrajectory& smoothedTrajectory)
return std::make_unique<Trk::FitQuality>(chiSquared, numberDoF);
}
-/*
- * Helper to combine forward with
- * smoother MultiComponentStates
- */
-Trk::MultiComponentState
-combine(const Trk::MultiComponentState& forwardsMultiState,
- const Trk::MultiComponentState& smootherMultiState,
- unsigned int maximumNumberOfComponents)
-{
-
- std::unique_ptr<Trk::MultiComponentState> combinedMultiState =
- std::make_unique<Trk::MultiComponentState>();
-
- // Loop over all components in forwards multi-state
- for (const auto& forwardsComponent : forwardsMultiState) {
- // Need to check that all components have associated weight matricies
- const AmgSymMatrix(5)* forwardMeasuredCov =
- forwardsComponent.first->covariance();
- // Loop over all components in the smoother multi-state
- for (const auto& smootherComponent : smootherMultiState) {
- // Need to check that all components have associated weight matricies
- const AmgSymMatrix(5)* smootherMeasuredCov =
- smootherComponent.first->covariance();
- if (!smootherMeasuredCov && !forwardMeasuredCov) {
- return {};
- }
-
- if (!forwardMeasuredCov) {
- Trk::ComponentParameters smootherComponentOnly(
- smootherComponent.first->clone(), smootherComponent.second);
- combinedMultiState->push_back(std::move(smootherComponentOnly));
- continue;
- }
-
- if (!smootherMeasuredCov) {
- Trk::ComponentParameters forwardComponentOnly(
- forwardsComponent.first->clone(), forwardsComponent.second);
- combinedMultiState->push_back(std::move(forwardComponentOnly));
- continue;
- }
-
- const AmgSymMatrix(5) summedCovariance =
- *forwardMeasuredCov + *smootherMeasuredCov;
- const AmgSymMatrix(5) K =
- *forwardMeasuredCov * summedCovariance.inverse();
- const AmgVector(5) newParameters =
- forwardsComponent.first->parameters() +
- K * (smootherComponent.first->parameters() -
- forwardsComponent.first->parameters());
- const AmgVector(5) parametersDiff =
- forwardsComponent.first->parameters() -
- smootherComponent.first->parameters();
-
- AmgSymMatrix(5) covarianceOfNewParameters =
- AmgSymMatrix(5)(K * *smootherMeasuredCov);
-
- std::unique_ptr<Trk::TrackParameters> combinedTrackParameters =
- (forwardsComponent.first)
- ->associatedSurface()
- .createUniqueTrackParameters(newParameters[Trk::loc1],
- newParameters[Trk::loc2],
- newParameters[Trk::phi],
- newParameters[Trk::theta],
- newParameters[Trk::qOverP],
- std::move(covarianceOfNewParameters));
- const AmgSymMatrix(5) invertedSummedCovariance =
- summedCovariance.inverse();
- // Determine the scaling factor for the new weighting. Determined from the
- // PDF of the many-dimensional gaussian
- double exponent =
- parametersDiff.transpose() * invertedSummedCovariance * parametersDiff;
- double weightScalingFactor = exp(-0.5 * exponent);
- double combinedWeight = smootherComponent.second *
- forwardsComponent.second * weightScalingFactor;
- Trk::ComponentParameters combinedComponent(
- std::move(combinedTrackParameters), combinedWeight);
- combinedMultiState->push_back(std::move(combinedComponent));
- }
- }
- // Component reduction on the combined state
- Trk::MultiComponentState mergedState =
- Trk::QuickCloseComponentsMultiStateMerger::merge(
- std::move(*combinedMultiState), maximumNumberOfComponents);
- // Before return the weights of the states need to be renormalised to one.
- Trk::MultiComponentStateHelpers::renormaliseState(mergedState);
-
- return mergedState;
-}
-
/*
* Helper to return the MultiComponent TSOS
* that we will push back to the Trajectory DataVector
@@ -156,8 +67,8 @@ combine(const Trk::MultiComponentState& forwardsMultiState,
*/
Trk::MultiComponentStateOnSurface*
smootherHelper(Trk::MultiComponentState&& updatedState,
- std::unique_ptr<const Trk::MeasurementBase> measurement,
- std::unique_ptr<Trk::FitQualityOnSurface> fitQuality,
+ std::unique_ptr<const Trk::MeasurementBase>&& measurement,
+ std::unique_ptr<Trk::FitQualityOnSurface>&& fitQuality,
bool islast)
{
@@ -396,7 +307,7 @@ Trk::GaussianSumFitter::fit(
}
// Perform GSF smoother operation
MultiTrajectory smoothedTrajectory =
- fit(ctx, extrapolatorCache, forwardTrajectory, particleHypothesis);
+ smootherFit(ctx, extrapolatorCache, forwardTrajectory, particleHypothesis);
if (smoothedTrajectory.empty()) {
m_SmootherFailure.fetch_add(1, std::memory_order_relaxed);
return nullptr;
@@ -506,7 +417,7 @@ Trk::GaussianSumFitter::fit(
// Perform GSF smoother operation
MultiTrajectory smoothedTrajectory =
- fit(ctx, extrapolatorCache, forwardTrajectory, particleHypothesis, ccot);
+ smootherFit(ctx, extrapolatorCache, forwardTrajectory, particleHypothesis, ccot);
if (smoothedTrajectory.empty()) {
m_SmootherFailure.fetch_add(1, std::memory_order_relaxed);
return nullptr;
@@ -960,7 +871,7 @@ Trk::GaussianSumFitter::stepForwardFit(
* Implementation of the smoothing of the trajectory.
*/
Trk::GaussianSumFitter::MultiTrajectory
-Trk::GaussianSumFitter::fit(
+Trk::GaussianSumFitter::smootherFit(
const EventContext& ctx,
Trk::IMultiStateExtrapolator::Cache& extrapolatorCache,
const MultiTrajectory& forwardTrajectory,
@@ -1109,7 +1020,7 @@ Trk::GaussianSumFitter::fit(
continue;
}
// Update with the measurement
- auto updatedState = Trk::GsfMeasurementUpdator::update(
+ updatedState = Trk::GsfMeasurementUpdator::update(
std::move(extrapolatedState), *measurement, fitQuality);
if (updatedState.empty()) {
ATH_MSG_WARNING(
@@ -1125,7 +1036,8 @@ Trk::GaussianSumFitter::fit(
const Trk::MultiComponentState& forwardsMultiState =
(*trackStateOnSurface)->components();
Trk::MultiComponentState combinedfitterState =
- combine(forwardsMultiState, updatedState, m_maximumNumberOfComponents);
+ Trk::MultiComponentStateCombiner::combine(
+ forwardsMultiState, updatedState, m_maximumNumberOfComponents);
if (combinedfitterState.empty()) {
ATH_MSG_WARNING("Could not combine state from forward fit with "
"smoother state");
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/TrkGaussianSumFilterUtils/MultiComponentStateCombiner.h b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/TrkGaussianSumFilterUtils/MultiComponentStateCombiner.h
index 33709d5744794ab136d4f103e9438f30aad3df65..5a053e76d12b11daec6ca1e9b57cab96299e1d14 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/TrkGaussianSumFilterUtils/MultiComponentStateCombiner.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/TrkGaussianSumFilterUtils/MultiComponentStateCombiner.h
@@ -47,6 +47,13 @@ combineCovWithWeight(const AmgVector(5) & firstParameters,
const AmgSymMatrix(5) & secondMeasuredCov,
const double secondWeight);
-}
+/** @brief Helper to combine forward with smoother MultiComponentStates
+ */
+Trk::MultiComponentState
+combine(const Trk::MultiComponentState& forwardsMultiState,
+ const Trk::MultiComponentState& smootherMultiState,
+ unsigned int maximumNumberOfComponents);
+
+}//end of MultiComponentStateCombiner namespace
} // end Trk namespace
#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/MultiComponentStateCombiner.cxx b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/MultiComponentStateCombiner.cxx
index 71af9e32c1b032446c5af2e29fb8a3519ab7db06..5e8ba142084bfc3ac71d9f2d232d2d4ac3008cc4 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/MultiComponentStateCombiner.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/MultiComponentStateCombiner.cxx
@@ -12,6 +12,7 @@
#include "TrkGaussianSumFilterUtils/MultiComponentStateCombiner.h"
#include "TrkGaussianSumFilterUtils/MultiComponentStateModeCalculator.h"
+#include "TrkGaussianSumFilterUtils/QuickCloseComponentsMultiStateMerger.h"
//
#include "CxxUtils/phihelper.h"
#include "TrkParameters/TrackParameters.h"
@@ -248,7 +249,8 @@ Trk::MultiComponentStateCombiner::combineWithWeight(
finalParameters, finalWeight, secondParameters, secondWeight);
const AmgSymMatrix(5)* firstMeasuredCov = firstTrackParameters->covariance();
- const AmgSymMatrix(5)* secondMeasuredCov = secondTrackParameters->covariance();
+ const AmgSymMatrix(5)* secondMeasuredCov =
+ secondTrackParameters->covariance();
// Check to see if first track parameters are measured or not
if (firstMeasuredCov && secondMeasuredCov) {
AmgSymMatrix(5) finalMeasuredCov(*firstMeasuredCov);
@@ -317,3 +319,95 @@ Trk::MultiComponentStateCombiner::combineCovWithWeight(
parameterDifference.transpose();
}
+// The following does heave use of Eigen
+// for covariance. Avoid out-of-line calls
+// to Eigen
+#if defined(__GNUC__)
+[[gnu::flatten]]
+#endif
+Trk::MultiComponentState
+Trk::MultiComponentStateCombiner::combine(
+ const Trk::MultiComponentState& forwardsMultiState,
+ const Trk::MultiComponentState& smootherMultiState,
+ unsigned int maximumNumberOfComponents)
+{
+
+ std::unique_ptr<Trk::MultiComponentState> combinedMultiState =
+ std::make_unique<Trk::MultiComponentState>();
+
+ // Loop over all components in forwards multi-state
+ for (const auto& forwardsComponent : forwardsMultiState) {
+ // Need to check that all components have associated weight matricies
+ const AmgSymMatrix(5)* forwardMeasuredCov =
+ forwardsComponent.first->covariance();
+ // Loop over all components in the smoother multi-state
+ for (const auto& smootherComponent : smootherMultiState) {
+ // Need to check that all components have associated weight matricies
+ const AmgSymMatrix(5)* smootherMeasuredCov =
+ smootherComponent.first->covariance();
+ if (!smootherMeasuredCov && !forwardMeasuredCov) {
+ return {};
+ }
+
+ if (!forwardMeasuredCov) {
+ Trk::ComponentParameters smootherComponentOnly(
+ smootherComponent.first->clone(), smootherComponent.second);
+ combinedMultiState->push_back(std::move(smootherComponentOnly));
+ continue;
+ }
+
+ if (!smootherMeasuredCov) {
+ Trk::ComponentParameters forwardComponentOnly(
+ forwardsComponent.first->clone(), forwardsComponent.second);
+ combinedMultiState->push_back(std::move(forwardComponentOnly));
+ continue;
+ }
+
+ const AmgSymMatrix(5) summedCovariance =
+ *forwardMeasuredCov + *smootherMeasuredCov;
+ const AmgSymMatrix(5) K =
+ *forwardMeasuredCov * summedCovariance.inverse();
+ const AmgVector(5) newParameters =
+ forwardsComponent.first->parameters() +
+ K * (smootherComponent.first->parameters() -
+ forwardsComponent.first->parameters());
+ const AmgVector(5) parametersDiff =
+ forwardsComponent.first->parameters() -
+ smootherComponent.first->parameters();
+
+ AmgSymMatrix(5) covarianceOfNewParameters =
+ AmgSymMatrix(5)(K * *smootherMeasuredCov);
+
+ std::unique_ptr<Trk::TrackParameters> combinedTrackParameters =
+ (forwardsComponent.first)
+ ->associatedSurface()
+ .createUniqueTrackParameters(newParameters[Trk::loc1],
+ newParameters[Trk::loc2],
+ newParameters[Trk::phi],
+ newParameters[Trk::theta],
+ newParameters[Trk::qOverP],
+ std::move(covarianceOfNewParameters));
+ const AmgSymMatrix(5) invertedSummedCovariance =
+ summedCovariance.inverse();
+ // Determine the scaling factor for the new weighting. Determined from the
+ // PDF of the many-dimensional gaussian
+ double exponent =
+ parametersDiff.transpose() * invertedSummedCovariance * parametersDiff;
+ double weightScalingFactor = exp(-0.5 * exponent);
+ double combinedWeight = smootherComponent.second *
+ forwardsComponent.second * weightScalingFactor;
+ Trk::ComponentParameters combinedComponent(
+ std::move(combinedTrackParameters), combinedWeight);
+ combinedMultiState->push_back(std::move(combinedComponent));
+ }
+ }
+ // Component reduction on the combined state
+ Trk::MultiComponentState mergedState =
+ Trk::QuickCloseComponentsMultiStateMerger::merge(
+ std::move(*combinedMultiState), maximumNumberOfComponents);
+ // Before return the weights of the states need to be renormalised to one.
+ Trk::MultiComponentStateHelpers::renormaliseState(mergedState);
+
+ return mergedState;
+}
+