From f268d6241aa5bfb5f49caea8fb27257749439701 Mon Sep 17 00:00:00 2001
From: Christos Anastopoulos <christos.anastopoulos@cern.ch>
Date: Fri, 23 Oct 2020 10:24:20 +0000
Subject: [PATCH] GSF : simplify the energy loss cache
---
.../TrkGaussianSumFilter/GsfConstants.h | 4 ++-
.../IBetheHeitlerEffects.h | 33 ++++++-----------
.../src/GsfBetheHeitlerEffects.cxx | 29 +++++++--------
.../src/GsfCombinedMaterialEffects.cxx | 36 +++++++------------
4 files changed, 37 insertions(+), 65 deletions(-)
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h
index 134b4fd3c833..200dbf99f144 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h
@@ -7,6 +7,8 @@
* @brief Collect constants we use
* in GSF and their meaning in one place
*/
+#ifndef GSFCONSTANTS_H
+#define GSFCONSTANTS_H
#include <cstdint>
namespace GSFConstants {
@@ -58,4 +60,4 @@ constexpr int8_t maxComponentsAfterConvolution =
maxNumberofBHComponents * maxNumberofStateComponents;
}
-
+#endif
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBetheHeitlerEffects.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBetheHeitlerEffects.h
index 8af03a16a9e7..d83706261f5b 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBetheHeitlerEffects.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/IBetheHeitlerEffects.h
@@ -16,36 +16,23 @@
#include "GaudiKernel/IAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
-#include "TrkMultiComponentStateOnSurface/ComponentParameters.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkEventPrimitives/PropDirection.h"
+#include "TrkGaussianSumFilter/GsfConstants.h"
+#include "TrkMultiComponentStateOnSurface/ComponentParameters.h"
+#include <array>
namespace Trk {
struct GSFEnergyLossCache
{
- GSFEnergyLossCache()
- {
- weights.reserve(6);
- deltaPs.reserve(6);
- deltaQOvePCov.reserve(6);
- }
-
- GSFEnergyLossCache(GSFEnergyLossCache&&) noexcept = default;
- GSFEnergyLossCache& operator=(GSFEnergyLossCache&&) noexcept = default;
- GSFEnergyLossCache(const GSFEnergyLossCache&) noexcept = default;
- GSFEnergyLossCache& operator=(const GSFEnergyLossCache&) noexcept = default;
- ~GSFEnergyLossCache() noexcept = default;
-
- std::vector<double> weights;
- std::vector<double> deltaPs;
- std::vector<double> deltaQOvePCov;
-
- void reset()
+ struct element
{
- weights.clear();
- deltaPs.clear();
- deltaQOvePCov.clear();
- }
+ double weight = 0;
+ double deltaP = 0;
+ double deltaQOvePCov = 0;
+ };
+ std::array<element, GSFConstants::maxNumberofBHComponents> elements = {};
+ int numElements = 0;
};
class MaterialProperties;
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
index ca42169aa614..0f69367e559d 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
@@ -24,7 +24,7 @@ namespace {
using BH = Trk::GsfBetheHeitlerEffects;
inline bool
-inRange(int var, int lo, int hi)
+inRange(int var, int lo, int hi)
{
return ((var <= hi) and (var >= lo));
}
@@ -356,8 +356,7 @@ Trk::GsfBetheHeitlerEffects::compute(
Trk::PropDirection direction,
Trk::ParticleHypothesis) const
{
- // Clear cache
- cache.reset();
+ cache.numElements = 0;
const Trk::TrackParameters* trackParameters = componentParameters.first.get();
const Amg::Vector3D& globalMomentum = trackParameters->momentum();
@@ -367,9 +366,8 @@ Trk::GsfBetheHeitlerEffects::compute(
double pathlengthInX0 = pathLength / radiationLength;
if (pathlengthInX0 < m_singleGaussianRange) {
- cache.weights.push_back(1.);
- cache.deltaPs.push_back(0.);
- cache.deltaQOvePCov.push_back(0.);
+ cache.elements[0] = { 1., 0., 0. };
+ cache.numElements = 1;
return;
}
@@ -390,9 +388,8 @@ Trk::GsfBetheHeitlerEffects::compute(
deltaP = sign * momentum * (1. / meanZ - 1.);
varQoverP = varZ / (momentum * momentum);
}
- cache.deltaPs.push_back(deltaP);
- cache.weights.push_back(1.);
- cache.deltaQOvePCov.push_back(varQoverP);
+ cache.elements[0] = { 1., deltaP, varQoverP };
+ cache.numElements = 1;
return;
}
@@ -467,29 +464,27 @@ Trk::GsfBetheHeitlerEffects::compute(
if (mixture[componentIndex].mean < m_componentMeanCut) {
continue;
}
+ double weight = mixture[componentIndex].weight;
if (componentIndex == componentWithHighestMean) {
- cache.weights.push_back(mixture[componentIndex].weight +
- weightToBeRemoved);
- } else {
- cache.weights.push_back(mixture[componentIndex].weight);
+ weight += weightToBeRemoved;
}
-
double deltaP(0.);
if (direction == alongMomentum) {
// For forward propagation
deltaP = momentum * (mixture[componentIndex].mean - 1.);
- cache.deltaPs.push_back(deltaP);
const double f = 1. / (momentum * mixture[componentIndex].mean);
varianceInverseMomentum = f * f * mixture[componentIndex].variance;
} // end forward propagation if clause
else {
// For backwards propagation
deltaP = momentum * (1. / mixture[componentIndex].mean - 1.);
- cache.deltaPs.push_back(deltaP);
varianceInverseMomentum =
mixture[componentIndex].variance / (momentum * momentum);
} // end backwards propagation if clause
- cache.deltaQOvePCov.push_back(varianceInverseMomentum);
+ cache.elements[cache.numElements] = { weight,
+ deltaP,
+ varianceInverseMomentum };
+ ++cache.numElements;
} // end for loop over all components
}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
index 22a07e1aa068..54db3b3bd118 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx
@@ -84,34 +84,24 @@ Trk::GsfCombinedMaterialEffects::compute(
// components
// we want at least on dummy to "combine"
// with scattering
- if (cache_energyLoss.weights.empty()) {
- cache_energyLoss.weights.push_back(1.);
- cache_energyLoss.deltaPs.push_back(0.);
- cache_energyLoss.deltaQOvePCov.push_back(0.);
+ if (cache_energyLoss.numElements == 0) {
+ cache_energyLoss.elements[0] = { 1, 0, 0 };
+ cache_energyLoss.numElements = 1;
}
/*
* 3. Combine the multiple scattering with each of the energy loss components
*/
- auto energyLoss_weightsIterator = cache_energyLoss.weights.begin();
- auto energyLoss_deltaPsIterator = cache_energyLoss.deltaPs.begin();
- auto energyLoss_deltaQOvePCovIterator =
- cache_energyLoss.deltaQOvePCov.begin();
-
- for (; energyLoss_weightsIterator != cache_energyLoss.weights.end();
- ++energyLoss_weightsIterator,
- ++energyLoss_deltaPsIterator,
- ++energyLoss_deltaQOvePCovIterator) {
-
- double combinedWeight = (*energyLoss_weightsIterator);
- double combinedDeltaP = (*energyLoss_deltaPsIterator);
+ for (int i = 0; i < cache_energyLoss.numElements; ++i) {
+
+ double combinedWeight = cache_energyLoss.elements[i].weight;
+ double combinedDeltaP = cache_energyLoss.elements[i].deltaP;
cache.weights.push_back(combinedWeight);
cache.deltaPs.push_back(combinedDeltaP);
-
if (measuredCov) {
// Create the covariance
const double covPhi = cache_multipleScatter.deltaPhiCov;
const double covTheta = cache_multipleScatter.deltaThetaCov;
- const double covQoverP = (*energyLoss_deltaQOvePCovIterator);
+ const double covQoverP = cache_energyLoss.elements[i].deltaQOvePCov;
AmgSymMatrix(5) cov;
cov << 0, 0, 0, 0, 0, // 5
0, 0, 0, 0, 0, // 10
@@ -173,8 +163,8 @@ Trk::GsfCombinedMaterialEffects::energyLoss(
PropDirection direction,
ParticleHypothesis particleHypothesis) const
{
- // Reset the cache
- cache.reset();
+
+ cache.numElements = 0;
// Request track parameters from component parameters
const Trk::TrackParameters* trackParameters = componentParameters.first.get();
@@ -198,7 +188,6 @@ Trk::GsfCombinedMaterialEffects::energyLoss(
// update for mean energy loss
const double deltaE = energyLoss ? energyLoss->deltaE() : 0;
const double sigmaDeltaE = energyLoss ? energyLoss->sigmaDeltaE() : 0;
-
// Calculate the pathlength encountered by the track
const double p = globalMomentum.mag();
const double m = s_particleMasses.mass[particleHypothesis];
@@ -208,7 +197,6 @@ Trk::GsfCombinedMaterialEffects::energyLoss(
// Calculate energy loss values uncertainty
const double sigmaQoverP = sigmaDeltaE / pow(beta * p, 2);
- cache.weights.push_back(1.);
- cache.deltaPs.push_back(deltaE);
- cache.deltaQOvePCov.push_back(sigmaQoverP * sigmaQoverP);
+ cache.elements[0] = { 1., deltaE, sigmaQoverP * sigmaQoverP };
+ cache.numElements = 1;
}
--
GitLab