#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)

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_BHnumberOfComponentsHighX0 {}

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 36 of file ElectronCombinedMaterialEffects.h.

Constructor & Destructor Documentation

◆ ElectronCombinedMaterialEffects()

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

Definition at line 221 of file ElectronCombinedMaterialEffects.cxx.

                                                      {
   // Read the std  polynomials
   {
     const 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());
     }
  
     readPolys readin = fillFromFile(resolvedFileName);
     m_BHnumberOfComponents = readin.numberOfComponents;
     m_BHpolynomialWeights = readin.weights;
     m_BHpolynomialMeans = readin.means;
     m_BHpolynomialVariances = readin.variances;
   }
   // Read the high X0 polynomials
   {
     const 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());
     }
     readPolys readin = fillFromFile(resolvedFileName);
     m_BHnumberOfComponentsHighX0 = readin.numberOfComponents;
     m_BHpolynomialWeightsHighX0 = readin.weights;
     m_BHpolynomialMeansHighX0 = readin.means;
     m_BHpolynomialVariancesHighX0 = readin.variances;
   }
   if (m_BHnumberOfComponentsHighX0 != m_BHnumberOfComponents) {
     std::ostringstream ss;
     ss << " numberOfComponentsHighX0 != numberOfComponents";
     throw std::logic_error(ss.str());
   }
 }

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 325 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 in x0 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;
   }
   //clip to upper range
   if (pathlengthInX0 > s_upperRange) {
     pathlengthInX0 = s_upperRange;
   }
  
   // Get proper mixture parameters
   MixtureParameters mixture;
   if (pathlengthInX0 > s_xOverRange) {
     mixture = getTransformedMixtureParameters(m_BHpolynomialWeightsHighX0,
                                               m_BHpolynomialMeansHighX0,
                                               m_BHpolynomialVariancesHighX0,
                                               pathlengthInX0,
                                               m_BHnumberOfComponents);
   } else {
     mixture = getTransformedMixtureParameters(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;
     } else { // For backwards propagation
       deltaP = momentum * (1. / mixture[componentIndex].mean - 1.);
       varianceInverseMomentum = mixture[componentIndex].variance / (momentum * momentum);
     }
  
     // 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 265 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
    * This for electrons comes from the Bethe-Heitler
    */
   GsfMaterial::EnergyLoss cache_energyLoss;
   this->BetheHeitler(cache_energyLoss,
                      componentParameters,
                      materialProperties,
                      pathLength,
                      direction);
   // Protect if there are no new energy loss
   // components.
   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.numEntries = 0; //cahce to be filled
   for (int i = 0; i < cache_energyLoss.numElements; ++i) {
     cache.weights[i] = cache_energyLoss.elements[i].weight;
     cache.deltaPs[i] = cache_energyLoss.elements[i].deltaP;
     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;
       //Here set the "delta" for the covariance
       //due to material effects
       cache.covariances[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.covariances[i].setZero();
     }
     ++cache.numEntries;
   }  // end for loop over energy loss components
 }