GSF : simplify the energy loss cache

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

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;

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
 }
 

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;
 }