Trk::ElectronCombinedMaterialEffects Class Reference

#include <ElectronCombinedMaterialEffects.h>

Collaboration diagram for Trk::ElectronCombinedMaterialEffects:

Classes
struct	ComponentValues
	Helper Struct for multiple Gaussian components. More...

Public Types
using	MixtureParameters = std::array< ComponentValues, GSFConstants::maxNumberofMatComponents >
using	Polynomial = std::array< double, GSFConstants::polynomialCoefficients >

Public Member Functions
	ElectronCombinedMaterialEffects (const std::string &parameterisationFileName, const std::string &parameterisationFileNameHighX0)
	ElectronCombinedMaterialEffects ()=default
	ElectronCombinedMaterialEffects (const ElectronCombinedMaterialEffects &)=default
	ElectronCombinedMaterialEffects (ElectronCombinedMaterialEffects &&)=default
ElectronCombinedMaterialEffects &	operator= (const ElectronCombinedMaterialEffects &)=default
ElectronCombinedMaterialEffects &	operator= (ElectronCombinedMaterialEffects &&)=default
	~ElectronCombinedMaterialEffects ()=default
void	compute (GsfMaterial::Combined &, const Trk::ComponentParameters &, const Trk::MaterialProperties &, double, Trk::PropDirection=anyDirection) const

Private Member Functions
void	BetheHeitler (GsfMaterial::EnergyLoss &cache, const ComponentParameters &componentParameters, const MaterialProperties &materialProperties, double pathLenght, PropDirection direction=anyDirection) const

Private Attributes
int	m_BHnumberOfComponents {}
int	m_BHtransformationCode {}
int	m_BHnumberOfComponentsHighX0 {}
int	m_BHtransformationCodeHighX0 {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialWeights {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialMeans {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialVariances {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialWeightsHighX0 {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialMeansHighX0 {}
std::array< Polynomial, GSFConstants::maxNumberofMatComponents >	m_BHpolynomialVariancesHighX0 {}

Detailed Description

Definition at line 25 of file ElectronCombinedMaterialEffects.h.

Member Typedef Documentation

MixtureParameters

using Trk::ElectronCombinedMaterialEffects::MixtureParameters = std::array<ComponentValues, GSFConstants::maxNumberofMatComponents>

Definition at line 35 of file ElectronCombinedMaterialEffects.h.

Polynomial

using Trk::ElectronCombinedMaterialEffects::Polynomial = std::array<double, GSFConstants::polynomialCoefficients>

Definition at line 38 of file ElectronCombinedMaterialEffects.h.

Constructor & Destructor Documentation

ElectronCombinedMaterialEffects() [1/4]

Trk::ElectronCombinedMaterialEffects::ElectronCombinedMaterialEffects	(	const std::string &	parameterisationFileName,
		const std::string &	parameterisationFileNameHighX0
	)

Definition at line 198 of file ElectronCombinedMaterialEffects.cxx.

{
  // The following is a bit repetitive code
  // we could consider refactoring
  // The low X0 polynomials
  {
    std::string resolvedFileName =
      PathResolver::find_file(parameterisationFileName, "DATAPATH");
    if (resolvedFileName.empty()) {
      std::ostringstream ss;
      ss << "Parameterisation file : " << parameterisationFileName
         << " not found";
      throw std::logic_error(ss.str());
    }
 
    const char* filename = resolvedFileName.c_str();
    std::ifstream fin(filename);
    if (fin.bad()) {
      std::ostringstream ss;
      ss << "Error opening file: " << resolvedFileName;
      throw std::logic_error(ss.str());
    }
 
    fin >> m_BHnumberOfComponents;
    int orderPolynomial = 0;
    fin >> orderPolynomial;
    fin >> m_BHtransformationCode;
    if (not inRange(
          m_BHnumberOfComponents, 0, GSFConstants::maxNumberofMatComponents)) {
      std::ostringstream ss;
      ss << "numberOfComponents Parameter out of range 0- "
         << GSFConstants::maxNumberofMatComponents << " : "
         << m_BHnumberOfComponents;
      throw std::logic_error(ss.str());
    }
    if (orderPolynomial != (GSFConstants::polynomialCoefficients - 1)) {
      std::ostringstream ss;
      ss << "orderPolynomial  order !=  "
         << (GSFConstants::polynomialCoefficients - 1);
      throw std::logic_error(ss.str());
    }
    if (not inRange(m_BHtransformationCode, 0, 1)) {
      std::ostringstream ss;
      ss << "transformationCode Parameter out of range 0-1: "
         << m_BHtransformationCode;
      throw std::logic_error(ss.str());
    }
    if (!fin) {
      std::ostringstream ss;
      ss << "Error while reading file : " << resolvedFileName;
      throw std::logic_error(ss.str());
    }
    // Fill the polynomials
    int componentIndex = 0;
    for (; componentIndex < m_BHnumberOfComponents; ++componentIndex) {
      m_BHpolynomialWeights[componentIndex] = readPolynomial(fin);
      m_BHpolynomialMeans[componentIndex] = readPolynomial(fin);
      m_BHpolynomialVariances[componentIndex] = readPolynomial(fin);
    }
  }
  // Read the high X0 polynomials
  {
    std::string resolvedFileName =
      PathResolver::find_file(parameterisationFileNameHighX0, "DATAPATH");
    if (resolvedFileName.empty()) {
      std::ostringstream ss;
      ss << "Parameterisation file : " << parameterisationFileNameHighX0
         << " not found";
      throw std::logic_error(ss.str());
    }
 
    const char* filename = resolvedFileName.c_str();
    std::ifstream fin(filename);
    if (fin.bad()) {
      std::ostringstream ss;
      ss << "Error opening file: " << resolvedFileName;
      throw std::logic_error(ss.str());
    }
    fin >> m_BHnumberOfComponentsHighX0;
    int orderPolynomial = 0;
    fin >> orderPolynomial;
    fin >> m_BHtransformationCodeHighX0;
    //
    if (not inRange(m_BHnumberOfComponentsHighX0,
                    0,
                    GSFConstants::maxNumberofMatComponents)) {
      std::ostringstream ss;
      ss << "numberOfComponentsHighX0 Parameter out of range 0- "
         << GSFConstants::maxNumberofMatComponents << " : "
         << m_BHnumberOfComponentsHighX0;
      throw std::logic_error(ss.str());
    }
    if (m_BHnumberOfComponentsHighX0 != m_BHnumberOfComponents) {
      std::ostringstream ss;
      ss << " numberOfComponentsHighX0 != numberOfComponents";
      throw std::logic_error(ss.str());
    }
    if (orderPolynomial != (GSFConstants::polynomialCoefficients - 1)) {
      std::ostringstream ss;
      ss << "orderPolynomial  order !=  "
         << (GSFConstants::polynomialCoefficients - 1);
      throw std::logic_error(ss.str());
    }
    if (not inRange(m_BHtransformationCodeHighX0, 0, 1)) {
      std::ostringstream ss;
      ss << "transformationCode Parameter out of range "
            "0-1: "
         << m_BHtransformationCodeHighX0;
      throw std::logic_error(ss.str());
    }
 
    // Fill the polynomials
    int componentIndex = 0;
    for (; componentIndex < m_BHnumberOfComponentsHighX0; ++componentIndex) {
      m_BHpolynomialWeightsHighX0[componentIndex] = readPolynomial(fin);
      m_BHpolynomialMeansHighX0[componentIndex] = readPolynomial(fin);
      m_BHpolynomialVariancesHighX0[componentIndex] = readPolynomial(fin);
    }
  }
}

ElectronCombinedMaterialEffects() [2/4]

Trk::ElectronCombinedMaterialEffects::ElectronCombinedMaterialEffects ( )

default

ElectronCombinedMaterialEffects() [3/4]

Trk::ElectronCombinedMaterialEffects::ElectronCombinedMaterialEffects ( const ElectronCombinedMaterialEffects &  )

default

ElectronCombinedMaterialEffects() [4/4]

Trk::ElectronCombinedMaterialEffects::ElectronCombinedMaterialEffects ( ElectronCombinedMaterialEffects &&  )

default

~ElectronCombinedMaterialEffects()

Trk::ElectronCombinedMaterialEffects::~ElectronCombinedMaterialEffects ( )

default

Member Function Documentation

BetheHeitler()

void Trk::ElectronCombinedMaterialEffects::BetheHeitler ( GsfMaterial::EnergyLoss &  cache, const ComponentParameters &  componentParameters, const MaterialProperties &  materialProperties, double  pathLenght, Trk::PropDirection  direction = anyDirection  ) const

private

Definition at line 385 of file ElectronCombinedMaterialEffects.cxx.

{
  cache.numElements = 0;
 
  const Trk::TrackParameters* trackParameters = componentParameters.params.get();
  const Amg::Vector3D& globalMomentum = trackParameters->momentum();
 
  const double radiationLength = materialProperties.x0();
  const double momentum = globalMomentum.mag();
  double pathlengthInX0 = pathLength / radiationLength;
 
  if (pathlengthInX0 < s_singleGaussianRange) {
    cache.elements[0] = { 1., 0., 0. };
    cache.numElements = 1;
    return;
  }
 
  // If the amount of material is between 0.0001 and 0.01 return the gaussian
  // approximation to the Bethe-Heitler distribution
  if (pathlengthInX0 < s_lowerRange) {
    const double meanZ = std::exp(-1. * pathlengthInX0);
    const double sign = (direction == Trk::oppositeMomentum) ? 1. : -1.;
    const double varZ =
      std::exp(-1. * pathlengthInX0 * std::log(3.) / std::log(2.)) -
      std::exp(-2. * pathlengthInX0);
    double deltaP(0.);
    double varQoverP(0.);
    if (direction == Trk::alongMomentum) {
      deltaP = sign * momentum * (1. - meanZ);
      varQoverP = 1. / (meanZ * meanZ * momentum * momentum) * varZ;
    } else {
      deltaP = sign * momentum * (1. / meanZ - 1.);
      varQoverP = varZ / (momentum * momentum);
    }
    cache.elements[0] = { 1., deltaP, varQoverP };
    cache.numElements = 1;
    return;
  }
 
  // Now we do the full calculation
  if (pathlengthInX0 > s_upperRange) {
    pathlengthInX0 = s_upperRange;
  }
 
  // Get proper mixture parameters
  MixtureParameters mixture;
  if (pathlengthInX0 > s_xOverRange) {
    if (m_BHtransformationCodeHighX0) {
      mixture = getTransformedMixtureParameters(m_BHpolynomialWeightsHighX0,
                                               m_BHpolynomialMeansHighX0,
                                               m_BHpolynomialVariancesHighX0,
                                               pathlengthInX0,
                                               m_BHnumberOfComponents);
    } else {
      mixture = getMixtureParameters(m_BHpolynomialWeightsHighX0,
                                     m_BHpolynomialMeansHighX0,
                                     m_BHpolynomialVariancesHighX0,
                                     pathlengthInX0,
                                     m_BHnumberOfComponents);
    }
  } else {
    if (m_BHtransformationCode) {
      mixture = getTransformedMixtureParameters(m_BHpolynomialWeights,
                                               m_BHpolynomialMeans,
                                               m_BHpolynomialVariances,
                                               pathlengthInX0,
                                               m_BHnumberOfComponents);
    } else {
      mixture = getMixtureParameters(m_BHpolynomialWeights,
                                     m_BHpolynomialMeans,
                                     m_BHpolynomialVariances,
                                     pathlengthInX0,
                                     m_BHnumberOfComponents);
    }
  }
  // Correct the mixture
  correctWeights(mixture, m_BHnumberOfComponents);
  int componentIndex = 0;
  double weightToBeRemoved(0.);
  int componentWithHighestMean(0);
  for (; componentIndex < m_BHnumberOfComponents; ++componentIndex) {
    if (mixture[componentIndex].mean > mixture[componentWithHighestMean].mean) {
      componentWithHighestMean = componentIndex;
    }
    if (mixture[componentIndex].mean >= s_componentMeanCut) {
      continue;
    }
    weightToBeRemoved += mixture[componentIndex].weight;
  }
  // Fill the cache to be returned
  componentIndex = 0;
  for (; componentIndex < m_BHnumberOfComponents; ++componentIndex) {
    double varianceInverseMomentum = 0;
    // This is not mathematically correct but it does stabilize the GSF
    if (mixture[componentIndex].mean < s_componentMeanCut) {
      continue;
    }
    double weight = mixture[componentIndex].weight;
    if (componentIndex == componentWithHighestMean) {
      weight += weightToBeRemoved;
    }
    double deltaP(0.);
    if (direction == alongMomentum) {
      // For forward propagation
      deltaP = momentum * (mixture[componentIndex].mean - 1.);
      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.);
      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 };
    ++cache.numElements;
  } // end for loop over all components
}

compute()

void Trk::ElectronCombinedMaterialEffects::compute	(	GsfMaterial::Combined &	cache,
		const Trk::ComponentParameters &	componentParameters,
		const Trk::MaterialProperties &	materialProperties,
		double	pathLength,
		Trk::PropDirection	direction = anyDirection
	)		const

Definition at line 322 of file ElectronCombinedMaterialEffects.cxx.

{
  const AmgSymMatrix(5)* measuredCov = componentParameters.params->covariance();
  /*
   * 1.  Retrieve multiple scattering corrections
   */
  GsfMaterial::Scattering cache_multipleScatter;
  scattering(
    cache_multipleScatter, componentParameters, materialProperties, pathLength);
  /*
   * 2. Retrieve energy loss corrections
   */
  GsfMaterial::EnergyLoss cache_energyLoss;
  this->BetheHeitler(cache_energyLoss,
                     componentParameters,
                     materialProperties,
                     pathLength,
                     direction);
  // Protect if there are no new energy loss
  // components
  // we want at least on dummy to "combine"
  // with scattering
  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
   */
  // Cache is to be filled so 0 entries here
  cache.numEntries = 0;
  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[i] = combinedWeight;
    cache.deltaPs[i] = combinedDeltaP;
    if (measuredCov) {
      // Create the covariance
      const double covPhi = cache_multipleScatter.deltaPhiCov;
      const double covTheta = cache_multipleScatter.deltaThetaCov;
      const double covQoverP = cache_energyLoss.elements[i].deltaQOvePCov;
      cache.deltaCovariances[i] << 0, 0, 0, 0, 0, // 5
        0, 0, 0, 0, 0,                            // 10
        0, 0, covPhi, 0, 0,                       // 15
        0, 0, 0, covTheta, 0,                     // 20
        0, 0, 0, 0, covQoverP;
    } else {
      cache.deltaCovariances[i].setZero();
    }
    ++cache.numEntries;
  } // end for loop over energy loss components
}

operator=() [1/2]

ElectronCombinedMaterialEffects& Trk::ElectronCombinedMaterialEffects::operator= ( const ElectronCombinedMaterialEffects &  )

default

operator=() [2/2]

ElectronCombinedMaterialEffects& Trk::ElectronCombinedMaterialEffects::operator= ( ElectronCombinedMaterialEffects &&  )

default

Member Data Documentation

m_BHnumberOfComponents

int Trk::ElectronCombinedMaterialEffects::m_BHnumberOfComponents {}

private

Definition at line 70 of file ElectronCombinedMaterialEffects.h.

m_BHnumberOfComponentsHighX0

int Trk::ElectronCombinedMaterialEffects::m_BHnumberOfComponentsHighX0 {}

private

Definition at line 72 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialMeans

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialMeans {}

private

Definition at line 78 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialMeansHighX0

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialMeansHighX0 {}

private

Definition at line 84 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialVariances

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialVariances {}

private

Definition at line 80 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialVariancesHighX0

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialVariancesHighX0 {}

private

Definition at line 86 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialWeights

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialWeights {}

private

Definition at line 76 of file ElectronCombinedMaterialEffects.h.

m_BHpolynomialWeightsHighX0

std::array<Polynomial, GSFConstants::maxNumberofMatComponents> Trk::ElectronCombinedMaterialEffects::m_BHpolynomialWeightsHighX0 {}

private

Definition at line 82 of file ElectronCombinedMaterialEffects.h.

m_BHtransformationCode

int Trk::ElectronCombinedMaterialEffects::m_BHtransformationCode {}

private

Definition at line 71 of file ElectronCombinedMaterialEffects.h.

m_BHtransformationCodeHighX0

int Trk::ElectronCombinedMaterialEffects::m_BHtransformationCodeHighX0 {}

private

Definition at line 73 of file ElectronCombinedMaterialEffects.h.

---

The documentation for this class was generated from the following files:

• ElectronCombinedMaterialEffects.h
• ElectronCombinedMaterialEffects.cxx