From 9ab0de717209622f1e338a6f778d14bd7ea04c72 Mon Sep 17 00:00:00 2001
From: christos <christos@cern.ch>
Date: Sun, 25 Oct 2020 04:11:05 +0100
Subject: [PATCH] GSF PosteriorWeightsCalculator : Reduce number of allocations
---
.../src/GsfBetheHeitlerEffects.cxx | 2 +
.../src/KLGaussianMixtureReduction.cxx | 2 +-
.../src/PosteriorWeightsCalculator.cxx | 52 ++++++++++++-------
3 files changed, 35 insertions(+), 21 deletions(-)
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
index 928bebe17ac..832e00f48d0 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
@@ -477,6 +477,8 @@ Trk::GsfBetheHeitlerEffects::compute(
varianceInverseMomentum =
mixture[componentIndex].variance / (momentum * momentum);
} // end backwards propagation if clause
+
+ //set in the cache and increase the elements
cache.elements[cache.numElements] = { weight,
deltaP,
varianceInverseMomentum };
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/KLGaussianMixtureReduction.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/KLGaussianMixtureReduction.cxx
index 5755e9564bc..b02026d63ad 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/KLGaussianMixtureReduction.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/KLGaussianMixtureReduction.cxx
@@ -221,7 +221,7 @@ findMerges(Component1D* componentsIn,
if (inputSize < 0 ||
inputSize > GSFConstants::maxComponentsAfterConvolution ||
reducedSize > inputSize) {
- throw std::runtime_error("Invalid InputSize or reducedSize");
+ throw std::runtime_error("findMerges :Invalid InputSize or reducedSize");
}
// We need just one for the full duration of a job
const static std::vector<triangularToIJ> convert = createToIJMaxRowCols();
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/PosteriorWeightsCalculator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/PosteriorWeightsCalculator.cxx
index 66bf47aa1dc..17d69368843 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/PosteriorWeightsCalculator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/PosteriorWeightsCalculator.cxx
@@ -13,8 +13,10 @@
#include "TrkEventPrimitives/FitQuality.h"
#include "TrkEventPrimitives/LocalParameters.h"
#include "TrkEventPrimitives/ProjectionMatricesSet.h"
+#include "TrkGaussianSumFilter/GsfConstants.h"
#include "TrkParameters/TrackParameters.h"
-
+#include <array>
+#include <stdexcept>
namespace {
using namespace Trk;
@@ -100,6 +102,17 @@ calculateWeight_2D_3(const TrackParameters* componentTrackParameters,
det, 0.5 * ((r.transpose() * R.inverse() * r)(0, 0)));
}
+struct componentsCache
+{
+ struct element
+ {
+ double determinantR;
+ double chi2;
+ };
+ std::array<element, GSFConstants::maxComponentsAfterConvolution> elements;
+ size_t numElements = 0;
+};
+
} // end of anonymous namespace
std::vector<Trk::ComponentParameters>
@@ -111,6 +124,10 @@ Trk::PosteriorWeightsCalculator::weights(MultiComponentState&& predictedState,
if (predictedStateSize == 0) {
return {};
}
+ if (predictedStateSize > GSFConstants::maxComponentsAfterConvolution) {
+ throw std::runtime_error(
+ "PosteriorWeightsCalculator :Invalid predictedState size");
+ }
const Trk::LocalParameters& measurementLocalParameters =
measurement.localParameters();
int nLocCoord = measurement.localParameters().dimension();
@@ -121,11 +138,7 @@ Trk::PosteriorWeightsCalculator::weights(MultiComponentState&& predictedState,
std::vector<Trk::ComponentParameters> returnMultiComponentState{};
returnMultiComponentState.reserve(predictedStateSize);
- std::vector<double> componentDeterminantR;
- componentDeterminantR.reserve(predictedStateSize);
- std::vector<double> componentChi2;
- componentChi2.reserve(predictedStateSize);
-
+ componentsCache determinantRandChi2{};
// Calculate chi2 and determinant of each component.
double minimumChi2(10.e10); // Initalise high
@@ -135,17 +148,15 @@ Trk::PosteriorWeightsCalculator::weights(MultiComponentState&& predictedState,
const Trk::TrackParameters* componentTrackParameters =
(*component).first.get();
-
if (!componentTrackParameters) {
continue;
}
-
const AmgSymMatrix(5)* predictedCov =
componentTrackParameters->covariance();
-
if (!predictedCov) {
continue;
}
+
std::pair<double, double> result(0, 0);
switch (nLocCoord) {
case 1: {
@@ -202,20 +213,19 @@ Trk::PosteriorWeightsCalculator::weights(MultiComponentState&& predictedState,
if (result.first == 0) {
continue;
}
- // Compute Chi2
- componentDeterminantR.push_back(result.first);
- componentChi2.push_back(result.second);
-
+ // Cache R and Chi2
+ determinantRandChi2.elements[determinantRandChi2.numElements] = {
+ result.first, result.second
+ };
+ ++determinantRandChi2.numElements;
if (result.second < minimumChi2) {
minimumChi2 = result.second;
}
} // end loop over components
- if (componentDeterminantR.size() != predictedState.size() ||
- componentChi2.size() != predictedState.size()) {
+ if (determinantRandChi2.numElements != predictedState.size()) {
return {};
}
-
// Calculate posterior weights.
unsigned int index(0);
double sumWeights(0.);
@@ -227,14 +237,16 @@ Trk::PosteriorWeightsCalculator::weights(MultiComponentState&& predictedState,
double priorWeight = (*component).second;
// Extract common factor to avoid numerical problems during exponentiation
- double chi2 = componentChi2[index] - minimumChi2;
+ double chi2 = determinantRandChi2.elements[index].chi2 - minimumChi2;
double updatedWeight(0.);
// Determinant can not be belowe 1e-19 in CLHEP .... rather ugly but protect
// against 0 determinants Normally occur when the component is a poor fit
- if (componentDeterminantR[index] > 1e-20) {
- updatedWeight = priorWeight * sqrt(1. / componentDeterminantR[index]) *
- exp(-0.5 * chi2);
+ if (determinantRandChi2.elements[index].determinantR > 1e-20) {
+ updatedWeight =
+ priorWeight *
+ sqrt(1. / determinantRandChi2.elements[index].determinantR) *
+ exp(-0.5 * chi2);
} else {
updatedWeight = 1e-10;
}
--
GitLab