#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
Protected Member Functions
virtual bool	newArguments (const TrajectoryStateOnSurface &, const PropagationDirection) const
	check of arguments for use with cached values
virtual void	storeArguments (const TrajectoryStateOnSurface &, const PropagationDirection) const
	storage of arguments for later use of
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)
double	BetheHeitlerVariance (const double rl) const
	Second moment of the Bethe-Heitler distribution (in z=E/E0)
virtual void	compute (const TrajectoryStateOnSurface &, const PropagationDirection) const
	Computation: generates vectors of weights, means and standard deviations.
double	correctedFirstMean (const double, const GSContainer &) const
	Correction for mean of component 1.
double	correctedFirstVar (const double, const GSContainer &) const
	Correction for variance of component 1.
void	correctWeights (GSContainer &) const
	Correction for weight of component 1.
void	getMixtureParameters (const double, GSContainer &) const
	Filling of mixture (in terms of z=E/E0)
double	logisticFunction (const double x) const
	Logistic function (needed for transformation of weight and mean)
void	readParameters (const std::string)
	Read parametrization from file.
Polynomial	readPolynomial (std::ifstream &, const int)
	Read coefficients of one polynomial from file.
Private Attributes
int	theCorrectionFlag
	parametrisation of variance for each component
float	theLastDz
	correction of 1st or 1st&2nd moments
float	theLastP
PropagationDirection	theLastPropDir
float	theLastRadLength
int	theNrComponents
std::vector< Polynomial >	thePolyMeans
	parametrisation of weight for each component
std::vector< Polynomial >	thePolyVars
	parametrisation of mean for each component
std::vector< Polynomial >	thePolyWeights
	values to be transformed by logistic / exp. function?
int	theTransformationCode
	number of components used for parameterisation

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.

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

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 [inline, private]

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 [inline, private]

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 [inline, virtual]

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 `[private, virtual]`

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(), AlCaHLTBitMon_ParallelJobs::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(), timingPdfMaker::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 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 [inline, private]

Logistic function (needed for transformation of weight and mean)

Definition at line 69 of file GsfBetheHeitlerUpdator.h.

References funct::exp().

Referenced by getMixtureParameters().

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

bool GsfBetheHeitlerUpdator::newArguments	(	const TrajectoryStateOnSurface &	TSoS,
		const PropagationDirection	propDir
	)		const `[protected, virtual]`

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 `[protected, virtual]`

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();
}