#include <GsfBetheHeitlerUpdator.h>

Inheritance diagram for GsfBetheHeitlerUpdator:

Classes
class	Polynomial

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

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
virtual std::vector < AlgebraicSymMatrix55 >	deltaLocalErrors (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir) const

virtual std::vector< double >	deltaPs (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir) const

	GsfMaterialEffectsUpdator (float mass)

float	mass () const

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

virtual std::vector< double >	weights (const TrajectoryStateOnSurface &TSoS, const PropagationDirection propDir) const

virtual	~GsfMaterialEffectsUpdator ()

Protected Member Functions
virtual bool	newArguments (const TrajectoryStateOnSurface &, const PropagationDirection) const
	check of arguments for use with cached values More...

virtual void	storeArguments (const TrajectoryStateOnSurface &, const PropagationDirection) const
	storage of arguments for later use of More...

Private Types
typedef std::vector< Triplet < double, double, double > >	GSContainer

Private Member Functions
double	BetheHeitlerMean (const double rl) const
	First moment of the Bethe-Heitler distribution (in z=E/E0) More...

double	BetheHeitlerVariance (const double rl) const
	Second moment of the Bethe-Heitler distribution (in z=E/E0) More...

virtual void	compute (const TrajectoryStateOnSurface &, const PropagationDirection) const
	Computation: generates vectors of weights, means and standard deviations. More...

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

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

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

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

double	logisticFunction (const double x) const
	Logistic function (needed for transformation of weight and mean) 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
	parametrisation of variance for each component More...

float	theLastDz
	correction of 1st or 1st&2nd moments More...

float	theLastP

PropagationDirection	theLastPropDir

float	theLastRadLength

int	theNrComponents

std::vector< Polynomial >	thePolyMeans
	parametrisation of weight for each component More...

std::vector< Polynomial >	thePolyVars
	parametrisation of mean for each component More...

std::vector< Polynomial >	thePolyWeights
	values to be transformed by logistic / exp. function? More...

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

Additional Inherited Members
Protected Attributes inherited from GsfMaterialEffectsUpdator
std::vector< AlgebraicSymMatrix55 >	theDeltaCovs

std::vector< double >	theDeltaPs

std::vector< double >	theWeights

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

Member Typedef Documentation

typedef std::vector< Triplet<double,double,double> > GsfBetheHeitlerUpdator::GSContainer

private

Definition at line 59 of file GsfBetheHeitlerUpdator.h.

Member Enumeration Documentation

enum GsfBetheHeitlerUpdator::CorrectionFlag

Enumerator
NoCorrection
MeanCorrection
FullCorrection

Definition at line 45 of file GsfBetheHeitlerUpdator.h.

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

GsfBetheHeitlerUpdator::FullCorrection

Definition: GsfBetheHeitlerUpdator.h:45

GsfBetheHeitlerUpdator::MeanCorrection

Definition: GsfBetheHeitlerUpdator.h:45

GsfBetheHeitlerUpdator::NoCorrection

Definition: GsfBetheHeitlerUpdator.h:45

Constructor & Destructor Documentation

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

constructor with explicit filename and correction flag

Definition at line 11 of file GsfBetheHeitlerUpdator.cc.

References readParameters(), and theCorrectionFlag.

Referenced by clone().

                                                      :
   //                           const CorrectionFlag correctionFlag) :
   GsfMaterialEffectsUpdator(0.000511),
   theNrComponents(0),
   theCorrectionFlag(correctionFlag),
   theLastDz(0.),
   theLastP(-1.),
   theLastPropDir(alongMomentum),
   theLastRadLength(-1.)
 {
   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);
 }

Member Function Documentation

double GsfBetheHeitlerUpdator::BetheHeitlerMean ( const double rl ) const

inlineprivate

First moment of the Bethe-Heitler distribution (in z=E/E0)

Definition at line 71 of file GsfBetheHeitlerUpdator.h.

References funct::exp().

Referenced by correctedFirstMean(), and correctedFirstVar().

   {
     return exp(-rl);
   }

double GsfBetheHeitlerUpdator::BetheHeitlerVariance ( const double rl ) const

inlineprivate

Second moment of the Bethe-Heitler distribution (in z=E/E0)

Definition at line 76 of file GsfBetheHeitlerUpdator.h.

References funct::exp(), and funct::log().

Referenced by correctedFirstVar().

   {
     return exp(-rl*log(3.)/log(2.)) - exp(-2*rl);
   }

virtual GsfBetheHeitlerUpdator* GsfBetheHeitlerUpdator::clone ( void ) const

inlinevirtual

Implements GsfMaterialEffectsUpdator.

Definition at line 48 of file GsfBetheHeitlerUpdator.h.

References GsfBetheHeitlerUpdator().

   {
     return new GsfBetheHeitlerUpdator(*this);
   }

void GsfBetheHeitlerUpdator::compute	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir
	)		const

privatevirtual

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

Implements GsfMaterialEffectsUpdator.

Definition at line 62 of file GsfBetheHeitlerUpdator.cc.

References alongMomentum, correctedFirstMean(), correctedFirstVar(), correctWeights(), benchmark_cfg::errors, f, first, getMixtureParameters(), i, TrajectoryStateOnSurface::localMomentum(), PV3DBase< T, PVType, FrameType >::mag(), Surface::mediumProperties(), L1TEmulatorMonitor_cff::p, MediumProperties::radLen(), edm::second(), storeArguments(), TrajectoryStateOnSurface::surface(), theCorrectionFlag, GsfMaterialEffectsUpdator::theDeltaCovs, GsfMaterialEffectsUpdator::theDeltaPs, theNrComponents, GsfMaterialEffectsUpdator::theWeights, and PV3DBase< T, PVType, FrameType >::z().

 {
   //
   // clear cache
   //
   theWeights.clear();
   theDeltaPs.clear();
   theDeltaCovs.clear();
   //
   // Get surface and check presence of medium properties
   //
   const Surface& surface = TSoS.surface();
   //
   // calculate components: first check associated material constants
   //
   double rl(0.);
   double p(0.);
   if ( surface.mediumProperties() ) {
     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.0001 ) {
     //
     // limit x/x0 to valid range for parametrisation
     // should be done in a more elegant way ...
     //
     if ( rl<0.01 )  rl = 0.01;
     if ( rl>0.20 )  rl = 0.20;
 
     GSContainer mixture;
     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++ ) {
       double varPinv;
       theWeights.push_back(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). 
     //
     theDeltaPs.push_back(p*(mixture[i].second-1));
     //    double f = 1./p/mixture[i].second/mixture[i].second;
     // patch to ensure consistency between for- and backward propagation
     double 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)
     //
     theDeltaPs.push_back(p*(1/mixture[i].second-1));
     varPinv = mixture[i].third/p/p;
       }
       AlgebraicSymMatrix55 errors;
       errors(0,0) = varPinv;
       theDeltaCovs.push_back(errors);
     }
   }
   else {
     theWeights.push_back(1.);
     theDeltaPs.push_back(0.);
     theDeltaCovs.push_back(AlgebraicSymMatrix55());
   }
   //
   // Save arguments to avoid duplication of computation
   //
   storeArguments(TSoS,propDir); 
 }

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

private

Correction for mean of component 1.

Definition at line 190 of file GsfBetheHeitlerUpdator.cc.

References BetheHeitlerMean(), i, max(), plotscripts::mean(), and min.

Referenced by compute().

 {
   if ( mixture.empty() )  return 0.;
   //
   // calculate difference true mean - weighted sum
   //
   double mean = BetheHeitlerMean(rl);
   for ( GSContainer::const_iterator i=mixture.begin()+1;
     i!=mixture.end(); i++ )  mean -= (*i).first*(*i).second;
   //
   // return corrected mean for first component
   //
   return std::max(std::min(mean/mixture[0].first,1.),0.);
 }

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

private

Correction for variance of component 1.

Definition at line 209 of file GsfBetheHeitlerUpdator.cc.

References BetheHeitlerMean(), BetheHeitlerVariance(), first, i, and max().

Referenced by compute().

 {
   if ( mixture.empty() )  return 0.;
   //
   // calculate difference true variance - weighted sum
   //
   double var = BetheHeitlerVariance(rl) +
     BetheHeitlerMean(rl)*BetheHeitlerMean(rl) -
     mixture[0].first*mixture[0].second*mixture[0].second;
   for ( GSContainer::const_iterator i=mixture.begin()+1;
     i!=mixture.end(); i++ )
     var -= (*i).first*((*i).second*(*i).second+(*i).third);
   //
   // return corrected variance for first component
   //
   return std::max(var/mixture[0].first,0.);
 }

void GsfBetheHeitlerUpdator::correctWeights ( GSContainer & mixture ) const

private

Correction for weight of component 1.

Definition at line 171 of file GsfBetheHeitlerUpdator.cc.

References i.

Referenced by compute().

 {
   if ( mixture.empty() )  return;
   //
   // get sum of weights
   //
   double wsum(0);
   for ( GSContainer::const_iterator i=mixture.begin();
     i!=mixture.end(); i++ )  wsum += (*i).first;
   //
   // rescale to obtain 1
   //
   for ( GSContainer::iterator i=mixture.begin();
     i!=mixture.end(); i++ )  (*i).first /= wsum;
 }

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

private

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

Definition at line 145 of file GsfBetheHeitlerUpdator.cc.

References funct::exp(), i, logisticFunction(), theNrComponents, thePolyMeans, thePolyVars, thePolyWeights, theTransformationCode, CommonMethods::weight(), and detailsBasic3DVector::z.

Referenced by compute().

 {
   mixture.clear();
   mixture.reserve(theNrComponents);
 
   for ( int i=0; i<theNrComponents; i++ ) {
 
     double weight = thePolyWeights[i](rl);
     if ( theTransformationCode )  weight = logisticFunction(weight);
 
     double z = thePolyMeans[i](rl);
     if ( theTransformationCode )  z = logisticFunction(z);
 
     double vz = thePolyVars[i](rl);
     if ( theTransformationCode )  vz = exp(vz);
     else                          vz = vz*vz;
 
     mixture.push_back(Triplet<double,double,double>(weight,z,vz));
   }
 }

double GsfBetheHeitlerUpdator::logisticFunction ( const double x ) const

inlineprivate

Logistic function (needed for transformation of weight and mean)

Definition at line 69 of file GsfBetheHeitlerUpdator.h.

References funct::exp().

Referenced by getMixtureParameters().

69 {return 1./(1.+exp(-x));}

funct::exp

Exp< T >::type exp(const T &t)

Definition: Exp.h:22

x

Definition: DDAxes.h:10

bool GsfBetheHeitlerUpdator::newArguments	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir
	)		const

protectedvirtual

check of arguments for use with cached values

Reimplemented from GsfMaterialEffectsUpdator.

Definition at line 232 of file GsfBetheHeitlerUpdator.cc.

References TrajectoryStateOnSurface::localMomentum(), PV3DBase< T, PVType, FrameType >::mag(), Surface::mediumProperties(), MediumProperties::radLen(), TrajectoryStateOnSurface::surface(), theLastDz, theLastP, theLastPropDir, theLastRadLength, Vector3DBase< T, FrameTag >::unit(), and PV3DBase< T, PVType, FrameType >::z().

                                                                 {
   return TSoS.localMomentum().unit().z()!=theLastDz ||
     TSoS.localMomentum().mag()!=theLastP || propDir!=theLastPropDir ||
     TSoS.surface().mediumProperties()->radLen()!=theLastRadLength;
 }

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

private

Read parametrization from file.

Definition at line 31 of file GsfBetheHeitlerUpdator.cc.

References convertXMLtoSQLite_cfg::fileName, edm::FileInPath::fullPath(), mergeVDriftHistosByStation::name, readPolynomial(), 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;
   for ( int ic=0; ic<theNrComponents; ic++ ) {
     thePolyWeights.push_back(readPolynomial(ifs,orderP));
     thePolyMeans.push_back(readPolynomial(ifs,orderP));
     thePolyVars.push_back(readPolynomial(ifs,orderP));
   }
 }

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

private

Read coefficients of one polynomial from file.

Definition at line 54 of file GsfBetheHeitlerUpdator.cc.

References i.

Referenced by readParameters().

                                      {
   std::vector<double> coeffs(order+1);
   for ( int i=0; i<(order+1); i++ ) aStream >> coeffs[i];
   return Polynomial(coeffs);
 }

void GsfBetheHeitlerUpdator::storeArguments	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir
	)		const

protectedvirtual

storage of arguments for later use of

Reimplemented from GsfMaterialEffectsUpdator.

Definition at line 241 of file GsfBetheHeitlerUpdator.cc.

References TrajectoryStateOnSurface::localMomentum(), PV3DBase< T, PVType, FrameType >::mag(), Surface::mediumProperties(), MediumProperties::radLen(), TrajectoryStateOnSurface::surface(), theLastDz, theLastP, theLastPropDir, theLastRadLength, Vector3DBase< T, FrameTag >::unit(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by compute().

                                                                    {
   theLastDz = TSoS.localMomentum().unit().z();
   theLastP = TSoS.localMomentum().mag();
   theLastPropDir = propDir;
   theLastRadLength = TSoS.surface().mediumProperties()->radLen();
 }