diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfConstants.h index 134b4fd3c833e2df984a90fc1192d6938d241d1e..200dbf99f144a6b630e9e42eb45dff6559ef35e1 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 8af03a16a9e7037318df80089e49e789bd1bcfb6..d83706261f5b2f9680d6b806b0b1bd686f3ae702 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 c4fc4676c8e27a537bc83b088a9568ea6a3a42c1..928bebe17ac527b0dfd5e0783a6c13a87a340eb9 100644 --- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx +++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx @@ -352,8 +352,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(); @@ -363,9 +362,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; } @@ -386,9 +384,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; } @@ -463,29 +460,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 22a07e1aa0686a35e9ad3ffe9452100901fcbc39..54db3b3bd1184d30a70cf4527ed13ac88886ae31 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; }