#include <GsfBetheHeitlerUpdator.h>

Inheritance diagram for GsfBetheHeitlerUpdator:

Classes
struct	GSContainer

class	Polynomial

Public Types
enum	CorrectionFlag { NoCorrection = 0, MeanCorrection = 1, FullCorrection = 2 }

Public Types inherited from GsfMaterialEffectsUpdator
typedef materialEffect::Covariance	Covariance

typedef materialEffect::CovIndex	CovIndex

typedef materialEffect::Effect	Effect

Public Member Functions
GsfBetheHeitlerUpdator *	clone () const override

	GsfBetheHeitlerUpdator (const std::string fileName, const int correctionFlag)
	constructor with explicit filename and correction flag More...

Public Member Functions inherited from GsfMaterialEffectsUpdator
	GsfMaterialEffectsUpdator (float mass, uint32_t is)

float	mass () const

size_t	size () const

virtual TrajectoryStateOnSurface	updateState (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir) const

virtual	~GsfMaterialEffectsUpdator ()

Private Member Functions
void	compute (const TrajectoryStateOnSurface &, const PropagationDirection, Effect[]) const override
	Computation: generates vectors of weights, means and standard deviations. More...

float	correctedFirstMean (const float, const GSContainer &) const
	Correction for mean of component 1. More...

float	correctedFirstVar (const float, const GSContainer &) const
	Correction for variance of component 1. More...

void	correctWeights (GSContainer &) const
	Correction for weight of component 1. More...

void	getMixtureParameters (const float, GSContainer &) const
	Filling of mixture (in terms of z=E/E0) More...

void	readParameters (const std::string)
	Read parametrization from file. More...

Polynomial	readPolynomial (std::ifstream &, const unsigned int)
	Read coefficients of one polynomial from file. More...

Private Attributes
int	theCorrectionFlag
	values to be transformed by logistic / exp. function? More...

int	theNrComponents

Polynomial	thePolyMeans [MaxSize]
	parametrisation of weight for each component More...

Polynomial	thePolyVars [MaxSize]
	parametrisation of mean for each componentP More...

Polynomial	thePolyWeights [MaxSize]
	correction of 1st or 1st&2nd moments More...

int	theTransformationCode
	number of components used for parameterisation More...

Static Private Attributes
static int	MaxOrder = 6

static int	MaxSize = 6

Additional Inherited Members
Protected Member Functions inherited from GsfMaterialEffectsUpdator
void	resize (size_t is)

Detailed Description

Description of electron energy loss according to Bethe-Heitler as a sum of Gaussian components. The weights and parameters of the Gaussians as a function of x/X0 are parametrized as polynomials. The coefficients of these polynomials are read from a file at construction time.

Definition at line 20 of file GsfBetheHeitlerUpdator.h.

Member Enumeration Documentation

enum GsfBetheHeitlerUpdator::CorrectionFlag

Enumerator
NoCorrection
MeanCorrection
FullCorrection

Definition at line 52 of file GsfBetheHeitlerUpdator.h.

52 { NoCorrection = 0, MeanCorrection = 1, FullCorrection = 2 };

GsfBetheHeitlerUpdator::FullCorrection

Definition: GsfBetheHeitlerUpdator.h:52

GsfBetheHeitlerUpdator::MeanCorrection

Definition: GsfBetheHeitlerUpdator.h:52

GsfBetheHeitlerUpdator::NoCorrection

Definition: GsfBetheHeitlerUpdator.h:52

Constructor & Destructor Documentation

GsfBetheHeitlerUpdator::GsfBetheHeitlerUpdator	(	const std::string	fileName,
		const int	correctionFlag
	)

constructor with explicit filename and correction flag

Definition at line 53 of file GsfBetheHeitlerUpdator.cc.

References readParameters(), GsfMaterialEffectsUpdator::resize(), theCorrectionFlag, and theNrComponents.

Referenced by clone().

     : GsfMaterialEffectsUpdator(0.000511, 6), theNrComponents(0), theCorrectionFlag(correctionFlag) {
   if (theCorrectionFlag == 1)
     edm::LogInfo("GsfBetheHeitlerUpdator") << "1st moment of mixture will be corrected";
   if (theCorrectionFlag >= 2)
     edm::LogInfo("GsfBetheHeitlerUpdator") << "1st and 2nd moments of mixture will be corrected";
 
   readParameters(fileName);
   assert(theNrComponents <= 6);
   resize(theNrComponents);
 }

Member Function Documentation

GsfBetheHeitlerUpdator* GsfBetheHeitlerUpdator::clone ( void ) const

inlineoverridevirtual

Implements GsfMaterialEffectsUpdator.

Definition at line 55 of file GsfBetheHeitlerUpdator.h.

References MillePedeFileConverter_cfg::fileName, GsfBetheHeitlerUpdator(), and AlCaHLTBitMon_QueryRunRegistry::string.

55 { return new GsfBetheHeitlerUpdator(*this); }

GsfBetheHeitlerUpdator::GsfBetheHeitlerUpdator

GsfBetheHeitlerUpdator(const std::string fileName, const int correctionFlag)

constructor with explicit filename and correction flag

Definition: GsfBetheHeitlerUpdator.cc:53

void GsfBetheHeitlerUpdator::compute	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir,
		Effect	effects[]
	)		const

overrideprivatevirtual

Computation: generates vectors of weights, means and standard deviations.

Implements GsfMaterialEffectsUpdator.

Definition at line 96 of file GsfBetheHeitlerUpdator.cc.

References alongMomentum, correctedFirstMean(), correctedFirstVar(), correctWeights(), materialEffect::Effect::deltaCov, materialEffect::Effect::deltaP, materialEffect::elos, f, getMixtureParameters(), mps_fire::i, MediumProperties::isValid(), TrajectoryStateOnSurface::localMomentum(), PV3DBase< T, PVType, FrameType >::mag(), Surface::mediumProperties(), AlCaHLTBitMon_ParallelJobs::p, MediumProperties::radLen(), TrajectoryStateOnSurface::surface(), theCorrectionFlag, theNrComponents, materialEffect::Effect::weight, and PV3DBase< T, PVType, FrameType >::z().

                                                              {
   //
   // Get surface and check presence of medium properties
   //
   const Surface& surface = TSoS.surface();
   //
   // calculate components: first check associated material constants
   //
   float rl(0.f);
   float p(0.f);
   if (surface.mediumProperties().isValid()) {
     LocalVector pvec = TSoS.localMomentum();
     p = pvec.mag();
     rl = surface.mediumProperties().radLen() / fabs(pvec.z()) * p;
   }
   //
   // produce multi-state only in case of x/X0>0
   //
   if (rl > 0.0001f) {
     //
     // limit x/x0 to valid range for parametrisation
     // should be done in a more elegant way ...
     //
     if (rl < 0.01f)
       rl = 0.01f;
     if (rl > 0.20f)
       rl = 0.20f;
 
     float mixtureData[3][theNrComponents];
     GSContainer mixture{mixtureData[0], mixtureData[1], mixtureData[2]};
 
     getMixtureParameters(rl, mixture);
     correctWeights(mixture);
     if (theCorrectionFlag >= 1)
       mixture.second[0] = correctedFirstMean(rl, mixture);
     if (theCorrectionFlag >= 2)
       mixture.third[0] = correctedFirstVar(rl, mixture);
 
     for (int i = 0; i < theNrComponents; i++) {
       float varPinv;
       effects[i].weight *= mixture.first[i];
       if (propDir == alongMomentum) {
         //
         // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
         // then convert sig(p) to sig(1/p).
         //
         effects[i].deltaP += p * (mixture.second[i] - 1.f);
         //    float f = 1./p/mixture.second[i]/mixture.second[i];
         // patch to ensure consistency between for- and backward propagation
         float f = 1.f / (p * mixture.second[i]);
         varPinv = f * f * mixture.third[i];
       } else {
         //
         // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside
         // convert to obtain equivalent delta(p)
         //
         effects[i].deltaP += p * (1.f / mixture.second[i] - 1.f);
         varPinv = mixture.third[i] / (p * p);
       }
       using namespace materialEffect;
       effects[i].deltaCov[elos] += varPinv;
     }
   }
 }

float GsfBetheHeitlerUpdator::correctedFirstMean	(	const float	rl,
		const GSContainer &	mixture
	)		const

private

Correction for mean of component 1.

Definition at line 207 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, mps_fire::i, SiStripPI::max, SiStripPI::mean, min(), GsfBetheHeitlerUpdator::GSContainer::second, and theNrComponents.

Referenced by compute().

                                                                                                  {
   //
   // calculate difference true mean - weighted sum
   //
   float mean = BetheHeitlerMean(rl);
   for (int i = 1; i < theNrComponents; i++)
     mean -= mixture.first[i] * mixture.second[i];
   //
   // return corrected mean for first component
   //
   return std::max(std::min(mean / mixture.first[0], 1.f), 0.f);
 }

float GsfBetheHeitlerUpdator::correctedFirstVar	(	const float	rl,
		const GSContainer &	mixture
	)		const

private

Correction for variance of component 1.

Definition at line 222 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, mps_fire::i, SiStripPI::max, GsfBetheHeitlerUpdator::GSContainer::second, theNrComponents, GsfBetheHeitlerUpdator::GSContainer::third, and trigObjTnPSource_cfi::var.

Referenced by compute().

                                                                                                 {
   //
   // calculate difference true variance - weighted sum
   //
   float var = BetheHeitlerVariance(rl) + BetheHeitlerMean(rl) * BetheHeitlerMean(rl) -
               mixture.first[0] * mixture.second[0] * mixture.second[0];
   for (int i = 1; i < theNrComponents; i++)
     var -= mixture.first[i] * (mixture.second[i] * mixture.second[i] + mixture.third[i]);
   //
   // return corrected variance for first component
   //
   return std::max(var / mixture.first[0], 0.f);
 }

void GsfBetheHeitlerUpdator::correctWeights ( GSContainer & mixture ) const

private

Correction for weight of component 1.

Definition at line 190 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, mps_fire::i, and theNrComponents.

Referenced by compute().

                                                                       {
   //
   // get sum of weights
   //
   float wsum(0);
   for (int i = 0; i < theNrComponents; i++)
     wsum += mixture.first[i];
   //
   // rescale to obtain 1
   //
   wsum = 1.f / wsum;
   for (int i = 0; i < theNrComponents; i++)
     mixture.first[i] *= wsum;
 }

void GsfBetheHeitlerUpdator::getMixtureParameters	(	const float	rl,
		GSContainer &	mixture
	)		const

private

Filling of mixture (in terms of z=E/E0)

Definition at line 165 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, mps_fire::i, GsfBetheHeitlerUpdator::GSContainer::second, theNrComponents, thePolyMeans, thePolyVars, thePolyWeights, theTransformationCode, GsfBetheHeitlerUpdator::GSContainer::third, and z.

Referenced by compute().

                                                                                             {
   float weight[theNrComponents], z[theNrComponents], vz[theNrComponents];
   for (int i = 0; i < theNrComponents; i++) {
     weight[i] = thePolyWeights[i](rl);
     z[i] = thePolyMeans[i](rl);
     vz[i] = thePolyVars[i](rl);
   }
   if (theTransformationCode)
     for (int i = 0; i < theNrComponents; i++) {
       mixture.first[i] = logisticFunction(weight[i]);
       mixture.second[i] = logisticFunction(z[i]);
       mixture.third[i] = unsafe_expf<4>(vz[i]);
       ;
     }
   else  // theTransformationCode
     for (int i = 0; i < theNrComponents; i++) {
       mixture.first[i] = weight[i];
       mixture.second[i] = z[i];
       mixture.third[i] = vz[i] * vz[i];
     }
 }

void GsfBetheHeitlerUpdator::readParameters ( const std::string fileName )

private

Read parametrization from file.

Definition at line 65 of file GsfBetheHeitlerUpdator.cc.

References MillePedeFileConverter_cfg::fileName, edm::FileInPath::fullPath(), MaxOrder, Skims_PA_cff::name, readPolynomial(), AlCaHLTBitMon_QueryRunRegistry::string, theNrComponents, thePolyMeans, thePolyVars, thePolyWeights, and theTransformationCode.

Referenced by GsfBetheHeitlerUpdator().

                                                                     {
   std::string name = "TrackingTools/GsfTracking/data/";
   name += fileName;
 
   edm::FileInPath parFile(name);
   edm::LogInfo("GsfBetheHeitlerUpdator") << "Reading GSF parameterization "
                                          << "of Bethe-Heitler energy loss from " << parFile.fullPath();
   std::ifstream ifs(parFile.fullPath().c_str());
 
   ifs >> theNrComponents;
   int orderP;
   ifs >> orderP;
   ifs >> theTransformationCode;
 
   assert(orderP < MaxOrder);
 
   for (int ic = 0; ic != theNrComponents; ++ic) {
     thePolyWeights[ic] = readPolynomial(ifs, orderP);
     thePolyMeans[ic] = readPolynomial(ifs, orderP);
     thePolyVars[ic] = readPolynomial(ifs, orderP);
   }
 }

GsfBetheHeitlerUpdator::Polynomial GsfBetheHeitlerUpdator::readPolynomial	(	std::ifstream &	aStream,
		const unsigned int	order
	)

private

Read coefficients of one polynomial from file.

Definition at line 88 of file GsfBetheHeitlerUpdator.cc.

References mps_fire::i.

Referenced by readParameters().

                                                                                                     {
   float coeffs[order + 1];
   for (unsigned int i = 0; i < (order + 1); ++i)
     aStream >> coeffs[i];
   return Polynomial(coeffs, order + 1);
 }