#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
virtual GsfBetheHeitlerUpdator *	clone () const

	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
virtual void	compute (const TrajectoryStateOnSurface &, const PropagationDirection, Effect[]) const
	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 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 22 of file GsfBetheHeitlerUpdator.h.

Member Enumeration Documentation

enum GsfBetheHeitlerUpdator::CorrectionFlag

Enumerator
NoCorrection
MeanCorrection
FullCorrection

Definition at line 55 of file GsfBetheHeitlerUpdator.h.

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

GsfBetheHeitlerUpdator::FullCorrection

Definition: GsfBetheHeitlerUpdator.h:55

GsfBetheHeitlerUpdator::MeanCorrection

Definition: GsfBetheHeitlerUpdator.h:55

GsfBetheHeitlerUpdator::NoCorrection

Definition: GsfBetheHeitlerUpdator.h:55

Constructor & Destructor Documentation

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

constructor with explicit filename and correction flag

Definition at line 56 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

virtual GsfBetheHeitlerUpdator* GsfBetheHeitlerUpdator::clone ( void ) const

inlinevirtual

Implements GsfMaterialEffectsUpdator.

Definition at line 58 of file GsfBetheHeitlerUpdator.h.

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

   {
     return new GsfBetheHeitlerUpdator(*this);
   }

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

privatevirtual

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

Implements GsfMaterialEffectsUpdator.

Definition at line 107 of file GsfBetheHeitlerUpdator.cc.

References alongMomentum, correctedFirstMean(), correctedFirstVar(), correctWeights(), materialEffect::Effect::deltaCov, materialEffect::Effect::deltaP, materialEffect::elos, f, getMixtureParameters(), 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 223 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, i, hpstanc_transforms::max, RecoTauDiscriminantConfiguration::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 241 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, i, hpstanc_transforms::max, GsfBetheHeitlerUpdator::GSContainer::second, theNrComponents, GsfBetheHeitlerUpdator::GSContainer::third, and JetChargeProducer_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 204 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, 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 176 of file GsfBetheHeitlerUpdator.cc.

References GsfBetheHeitlerUpdator::GSContainer::first, 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 73 of file GsfBetheHeitlerUpdator.cc.

References MillePedeFileConverter_cfg::fileName, edm::FileInPath::fullPath(), MaxOrder, dataset::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 int	order
	)

private

Read coefficients of one polynomial from file.

Definition at line 99 of file GsfBetheHeitlerUpdator.cc.

References i.

Referenced by readParameters().

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