diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfBetheHeitlerEffects.h b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfBetheHeitlerEffects.h
index 591ad30217889a8f39e992dd4ff593c09791be76..4dd080aa15b834d7bb5d7fa20a35a746ce0a48f9 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfBetheHeitlerEffects.h
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/GsfBetheHeitlerEffects.h
@@ -27,15 +27,43 @@ class GsfBetheHeitlerEffects
{
public:
+ static constexpr int maxNumberofComponents = 8;
+ static constexpr int maxOrderPolynomial = 10;
+ struct ComponentValues
+ {
+ // Default ctors/dtor/assignment operators
+ ComponentValues() = default;
+ ~ComponentValues() = default;
+ ComponentValues(const ComponentValues&) = default;
+ ComponentValues& operator=(const ComponentValues&) = default;
+ ComponentValues(ComponentValues&&) = default;
+ ComponentValues& operator=(ComponentValues&&) = default;
+ // Constructor with arguments
+ ComponentValues(double aWeight, double aMean, double aVariance)
+ : weight(aWeight)
+ , mean(aMean)
+ , variance(aVariance)
+ {}
+ double weight;
+ double mean;
+ double variance;
+ };
+
+ using MixtureParameters = std::array<ComponentValues, maxNumberofComponents>;
+
/** Helper class for construction and evaluation of polynomial */
class Polynomial
{
public:
- // Default constructor
Polynomial() = default;
-
- /** Constructor from a vector of coefficients (in decreasing order of powers
- * of x */
+ ~Polynomial() = default;
+ Polynomial(const Polynomial&) = default;
+ Polynomial& operator=(const Polynomial&) = default;
+ Polynomial(Polynomial&&) = default;
+ Polynomial& operator=(Polynomial&&) = default;
+
+ /** Constructor from a vector of coefficients
+ * (in decreasing order of powers of x) */
Polynomial(const std::vector<double>& coefficients)
: m_coefficients(coefficients){};
@@ -56,26 +84,6 @@ public:
std::vector<double> m_coefficients;
};
- struct ComponentValues
- {
- // Default ctors/dtor/assignment operators
- ComponentValues() = default;
- ~ComponentValues() = default;
- ComponentValues(const ComponentValues&) = default;
- ComponentValues& operator=(const ComponentValues&) = default;
- ComponentValues(ComponentValues&&) = default;
- ComponentValues& operator=(ComponentValues&&) = default;
- // Constructor with arguments
- ComponentValues(double aWeight, double aMean, double aVariance)
- : weight(aWeight)
- , mean(aMean)
- , variance(aVariance)
- {}
- double weight;
- double mean;
- double variance;
- };
-
GsfBetheHeitlerEffects(const std::string&,
const std::string&,
const IInterface*);
@@ -93,16 +101,13 @@ public:
ParticleHypothesis particleHypothesis =
nonInteracting) const override final;
- typedef std::vector<ComponentValues> MixtureParameters;
private:
-
// Read polynomial fit parameters from a specified file
bool readParameters();
// Read coeffients for a single polynomial fit
Polynomial readPolynomial(std::ifstream&, const int);
-
std::vector<Polynomial> m_polynomialWeights;
std::vector<Polynomial> m_polynomialMeans;
std::vector<Polynomial> m_polynomialVariances;
@@ -110,7 +115,6 @@ private:
std::vector<Polynomial> m_polynomialMeansHighX0;
std::vector<Polynomial> m_polynomialVariancesHighX0;
-
int m_numberOfComponents;
int m_transformationCode;
int m_correctionFlag;
@@ -122,11 +126,10 @@ private:
double m_xOverRange;
double m_upperRange;
double m_componentMeanCut;
-
+
bool m_useHighX0;
std::string m_parameterisationFileName;
std::string m_parameterisationFileNameHighX0;
-
};
}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
index 3c6063f1c849aa314188a699efce428357dcf5ff..4607881ddba3f9e457b4e6ea0da367f84eb89afa 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilter/src/GsfBetheHeitlerEffects.cxx
@@ -49,47 +49,39 @@ betheHeitlerVariance(const double r)
}
void
-correctWeights(Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture)
+correctWeights(Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture,
+ const int numberOfComponents)
{
- if (mixture.empty()) {
+ if (numberOfComponents < 1) {
return;
}
-
// Obtain the sum of weights
double weightSum(0.);
-
- Trk::GsfBetheHeitlerEffects::MixtureParameters::const_iterator component =
- mixture.begin();
- for (; component != mixture.end(); ++component) {
- weightSum += (*component).weight;
+ for (int i = 0; i < numberOfComponents; ++i) {
+ weightSum += mixture[i].weight;
}
-
// Rescale so that total weighting is 1
- Trk::GsfBetheHeitlerEffects::MixtureParameters::iterator modifiableComponent =
- mixture.begin();
-
- for (; modifiableComponent != mixture.end(); ++modifiableComponent) {
- (*modifiableComponent).weight /= weightSum;
+ for (int i = 0; i < numberOfComponents; ++i) {
+ mixture[i].weight /= weightSum;
}
}
double
correctedFirstMean(
const double pathlengthInX0,
- const Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture)
+ const Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture,
+ const int numberOfComponents)
{
-
- if (mixture.empty()) {
- return 0.;
+ if (numberOfComponents < 1) {
+ return 0;
}
// Obtain the difference between the true and weighted sum
double meanBH = betheHeitlerMean(pathlengthInX0);
- Trk::GsfBetheHeitlerEffects::MixtureParameters::const_iterator component =
- mixture.begin() + 1;
- for (; component != mixture.end(); ++component) {
- meanBH -= (*component).weight * (*component).mean;
+ for (int i = 1; i < numberOfComponents; ++i) {
+ meanBH -= mixture[i].weight * mixture[i].mean;
}
+
// return the corrected mean for the first component
return std::max(std::min(meanBH / mixture[0].weight, 1.), 0.);
}
@@ -97,25 +89,25 @@ correctedFirstMean(
double
correctedFirstVariance(
const double pathlengthInX0,
- const Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture)
+ const Trk::GsfBetheHeitlerEffects::MixtureParameters& mixture,
+ const int numberOfComponents)
{
-
- if (mixture.empty()) {
- return 0.;
+ if (numberOfComponents < 1) {
+ return 0;
}
// Obtain the difference between the true and weighed sum variances
double varianceBH =
betheHeitlerVariance(pathlengthInX0) +
(betheHeitlerMean(pathlengthInX0) * betheHeitlerMean(pathlengthInX0));
+
varianceBH -= mixture[0].weight * mixture[0].mean * mixture[0].mean;
- Trk::GsfBetheHeitlerEffects::MixtureParameters::const_iterator component =
- mixture.begin() + 1;
- for (; component != mixture.end(); ++component) {
- varianceBH -= (*component).weight * ((*component).mean * (*component).mean +
- (*component).variance);
+ for (int i = 1; i < numberOfComponents; ++i) {
+ varianceBH -= mixture[i].weight *
+ (mixture[i].mean * mixture[i].mean + mixture[i].variance);
}
+
return std::max(varianceBH / mixture[0].weight, 0.);
}
@@ -123,20 +115,20 @@ Trk::GsfBetheHeitlerEffects::MixtureParameters
getTranformedMixtureParameters(
const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialWeights,
const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialMeans,
- const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialVariances,
+ const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>&
+ polynomialVariances,
const double pathlengthInX0,
const int numberOfComponents)
{
- Trk::GsfBetheHeitlerEffects::MixtureParameters mixture;
- mixture.reserve(numberOfComponents);
+ Trk::GsfBetheHeitlerEffects::MixtureParameters mixture{};
for (int i = 0; i < numberOfComponents; ++i) {
const double updatedWeight = polynomialWeights[i](pathlengthInX0);
const double updatedMean = polynomialMeans[i](pathlengthInX0);
const double updatedVariance = polynomialVariances[i](pathlengthInX0);
- mixture.emplace_back(logisticFunction(updatedWeight),
- logisticFunction(updatedMean),
- exp(updatedVariance));
+ mixture[i] = { logisticFunction(updatedWeight),
+ logisticFunction(updatedMean),
+ exp(updatedVariance) };
}
return mixture;
}
@@ -145,22 +137,24 @@ Trk::GsfBetheHeitlerEffects::MixtureParameters
getMixtureParameters(
const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialWeights,
const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialMeans,
- const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>& polynomialVariances,
+ const std::vector<Trk::GsfBetheHeitlerEffects::Polynomial>&
+ polynomialVariances,
const double pathlengthInX0,
const int numberOfComponents)
{
- Trk::GsfBetheHeitlerEffects::MixtureParameters mixture;
- mixture.reserve(numberOfComponents);
+ Trk::GsfBetheHeitlerEffects::MixtureParameters mixture{};
for (int i = 0; i < numberOfComponents; ++i) {
const double updatedWeight = polynomialWeights[i](pathlengthInX0);
const double updatedMean = polynomialMeans[i](pathlengthInX0);
const double updatedVariance = polynomialVariances[i](pathlengthInX0);
- mixture.emplace_back(
- updatedWeight, updatedMean, updatedVariance * updatedVariance);
+ mixture[i] = { updatedWeight,
+ updatedMean,
+ updatedVariance * updatedVariance };
}
return mixture;
}
+
} // end of Anonymous namespace for Helper methods
Trk::GsfBetheHeitlerEffects::GsfBetheHeitlerEffects(const std::string& type,
@@ -217,6 +211,8 @@ Trk::GsfBetheHeitlerEffects::initialize()
bool
Trk::GsfBetheHeitlerEffects::readParameters()
{
+
+ // Read std polynomial
std::string resolvedFileName =
PathResolver::find_file(m_parameterisationFileName, "DATAPATH");
if (!resolvedFileName.empty()) {
@@ -238,14 +234,14 @@ Trk::GsfBetheHeitlerEffects::readParameters()
fin >> orderPolynomial;
fin >> m_transformationCode;
//
- if (not inRange(m_numberOfComponents, 0, 8)) {
- ATH_MSG_ERROR("numberOfComponents Parameter out of range 0-8: "
- << m_numberOfComponents);
+ if (not inRange(m_numberOfComponents, 0, maxNumberofComponents)) {
+ ATH_MSG_ERROR("numberOfComponents Parameter out of range 0- "
+ << maxNumberofComponents << " : " << m_numberOfComponents);
return false;
}
- if (not inRange(orderPolynomial, 0, 10)) {
- ATH_MSG_ERROR(
- "orderPolynomial Parameter out of range 0-10: " << orderPolynomial);
+ if (not inRange(orderPolynomial, 0, maxOrderPolynomial)) {
+ ATH_MSG_ERROR("orderPolynomial Parameter out of range 0- "
+ << maxOrderPolynomial << " : " << orderPolynomial);
return false;
}
if (not inRange(m_transformationCode, 0, 10)) {
@@ -265,6 +261,7 @@ Trk::GsfBetheHeitlerEffects::readParameters()
m_polynomialVariances.push_back(readPolynomial(fin, orderPolynomial));
}
+ // Read the high X0 polynomial
if (m_useHighX0) {
resolvedFileName =
PathResolver::find_file(m_parameterisationFileNameHighX0, "DATAPATH");
@@ -288,14 +285,20 @@ Trk::GsfBetheHeitlerEffects::readParameters()
fin >> orderPolynomial;
fin >> m_transformationCodeHighX0;
//
- if (not inRange(m_numberOfComponentsHighX0, 0, 8)) {
- ATH_MSG_ERROR("numberOfComponents Parameter out of range 0-8: "
+ if (not inRange(m_numberOfComponentsHighX0, 0, maxNumberofComponents)) {
+ ATH_MSG_ERROR("numberOfComponentsHighX0 Parameter out of range 0- "
+ << maxNumberofComponents << " : "
<< m_numberOfComponentsHighX0);
return false;
}
- if (not inRange(orderPolynomial, 0, 10)) {
- ATH_MSG_ERROR(
- "orderPolynomial Parameter out of range 0-10: " << orderPolynomial);
+ if (m_numberOfComponentsHighX0 != m_numberOfComponents) {
+ ATH_MSG_ERROR(" numberOfComponentsHighX0 != numberOfComponents");
+ return false;
+ }
+
+ if (not inRange(orderPolynomial, 0, maxOrderPolynomial)) {
+ ATH_MSG_ERROR("orderPolynomial Parameter out of range 0- "
+ << maxOrderPolynomial << " : " << orderPolynomial);
return false;
}
if (not inRange(m_transformationCodeHighX0, 0, 10)) {
@@ -354,9 +357,6 @@ Trk::GsfBetheHeitlerEffects::compute(
double pathlengthInX0 = pathLength / radiationLength;
if (pathlengthInX0 < m_singleGaussianRange) {
- ATH_MSG_DEBUG("Trying to apply energy loss to "
- << pathlengthInX0
- << " x/x0. No Bethe-Heitler effects applied");
cache.weights.push_back(1.);
cache.deltaPs.push_back(0.);
cache.deltaQOvePCov.push_back(0.);
@@ -423,17 +423,19 @@ Trk::GsfBetheHeitlerEffects::compute(
}
// Correct the mixture
- correctWeights(mixture);
+ correctWeights(mixture, m_numberOfComponents);
if (m_correctionFlag == 1) {
- mixture[0].mean = correctedFirstMean(pathlengthInX0, mixture);
+ mixture[0].mean =
+ correctedFirstMean(pathlengthInX0, mixture, m_numberOfComponents);
}
if (m_correctionFlag == 2) {
- mixture[0].mean = correctedFirstMean(pathlengthInX0, mixture);
- mixture[0].variance = correctedFirstVariance(pathlengthInX0, mixture);
+ mixture[0].mean =
+ correctedFirstMean(pathlengthInX0, mixture, m_numberOfComponents);
+ mixture[0].variance =
+ correctedFirstVariance(pathlengthInX0, mixture, m_numberOfComponents);
}
- //
int componentIndex = 0;
double weightToBeRemoved(0.);
int componentWithHighestMean(0);
@@ -476,10 +478,6 @@ Trk::GsfBetheHeitlerEffects::compute(
varianceInverseMomentum =
mixture[componentIndex].variance / (momentum * momentum);
} // end backwards propagation if clause
-
- AmgSymMatrix(5) newCovarianceMatrix;
- newCovarianceMatrix.setZero();
- newCovarianceMatrix(Trk::qOverP, Trk::qOverP) = varianceInverseMomentum;
cache.deltaQOvePCov.push_back(varianceInverseMomentum);
} // end for loop over all components
}