#include <GsfBetheHeitlerUpdator.h>

Inheritance diagram for GsfBetheHeitlerUpdator:

Classes
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 Types
typedef Triplet< float, float, float >	GSContainer

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 26 of file GsfBetheHeitlerUpdator.h.

Member Typedef Documentation

typedef Triplet<float,float,float> GsfBetheHeitlerUpdator::GSContainer

private

Definition at line 72 of file GsfBetheHeitlerUpdator.h.

Member Enumeration Documentation

enum GsfBetheHeitlerUpdator::CorrectionFlag

Enumerator
NoCorrection
MeanCorrection
FullCorrection

Definition at line 59 of file GsfBetheHeitlerUpdator.h.

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

GsfBetheHeitlerUpdator::FullCorrection

Definition: GsfBetheHeitlerUpdator.h:59

GsfBetheHeitlerUpdator::MeanCorrection

Definition: GsfBetheHeitlerUpdator.h:59

GsfBetheHeitlerUpdator::NoCorrection

Definition: GsfBetheHeitlerUpdator.h:59

Constructor & Destructor Documentation

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

constructor with explicit filename and correction flag

Definition at line 46 of file GsfBetheHeitlerUpdator.cc.

References assert(), 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 62 of file GsfBetheHeitlerUpdator.h.

References GsfBetheHeitlerUpdator().

   {
     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 97 of file GsfBetheHeitlerUpdator.cc.

References alongMomentum, correctedFirstMean(), correctedFirstVar(), correctWeights(), materialEffect::Effect::deltaCov, materialEffect::Effect::deltaP, materialEffect::elos, f, Triplet< T1, T2, T3 >::first, getMixtureParameters(), i, MediumProperties::isValid(), TrajectoryStateOnSurface::localMomentum(), PV3DBase< T, PVType, FrameType >::mag(), Surface::mediumProperties(), AlCaHLTBitMon_ParallelJobs::p, MediumProperties::radLen(), Triplet< T1, T2, T3 >::second, TrajectoryStateOnSurface::surface(), theCorrectionFlag, theNrComponents, Triplet< T1, T2, T3 >::third, 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;
 
     #if __clang__
     std::vector<GSContainer> mixtureHolder(theNrComponents);
     GSContainer *mixture = mixtureHolder.data();
     #else
     GSContainer mixture[theNrComponents];
     #endif
     getMixtureParameters(rl,mixture);
     correctWeights(mixture);
     if ( theCorrectionFlag>=1 )
       mixture[0].second = correctedFirstMean(rl,mixture);
     if ( theCorrectionFlag>=2 )
       mixture[0].third = correctedFirstVar(rl,mixture);
 
     for ( int i=0; i<theNrComponents; i++ ) {
       float varPinv;
       effects[i].weight*=mixture[i].first;
       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[i].second-1);
     //    float f = 1./p/mixture[i].second/mixture[i].second;
     // patch to ensure consistency between for- and backward propagation
     float f = 1./p/mixture[i].second;
     varPinv = f*f*mixture[i].third;
       }
       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/mixture[i].second-1);
     varPinv = mixture[i].third/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 212 of file GsfBetheHeitlerUpdator.cc.

References f, plotBeamSpotDB::first, i, bookConverter::max, timingPdfMaker::mean, min(), edm::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[i].first*mixture[i].second;
   //
   // return corrected mean for first component
   //
   return std::max(std::min(mean/mixture[0].first,1.f),0.f);
 }

float GsfBetheHeitlerUpdator::correctedFirstVar	(	const float	rl,
		const GSContainer	mixture[]
	)		const

private

Correction for variance of component 1.

Definition at line 230 of file GsfBetheHeitlerUpdator.cc.

References f, Triplet< T1, T2, T3 >::first, plotBeamSpotDB::first, i, bookConverter::max, Triplet< T1, T2, T3 >::second, edm::second(), theNrComponents, and MetTreeProducer::var().

Referenced by compute().

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

void GsfBetheHeitlerUpdator::correctWeights ( GSContainer mixture[] ) const

private

Correction for weight of component 1.

Definition at line 194 of file GsfBetheHeitlerUpdator.cc.

References plotBeamSpotDB::first, i, and theNrComponents.

Referenced by compute().

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

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

private

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

Definition at line 169 of file GsfBetheHeitlerUpdator.cc.

References i, theNrComponents, thePolyMeans, thePolyVars, thePolyWeights, theTransformationCode, histoStyle::weight, and detailsBasic3DVector::z.

Referenced by compute().

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

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

private

Read parametrization from file.

Definition at line 63 of file GsfBetheHeitlerUpdator.cc.

References assert(), convertXMLtoSQLite_cfg::fileName, edm::FileInPath::fullPath(), MaxOrder, mergeVDriftHistosByStation::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 89 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);
 }