GsfBetheHeitlerUpdator Class Reference

#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 constexpr int	MaxOrder = 6
static constexpr 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

CorrectionFlag

enum GsfBetheHeitlerUpdator::CorrectionFlag

Enumerator
NoCorrection
MeanCorrection
FullCorrection

Definition at line 52 of file GsfBetheHeitlerUpdator.h.

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

Constructor & Destructor Documentation

GsfBetheHeitlerUpdator()

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 cms::cuda::assert(), MillePedeFileConverter_cfg::fileName, 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

clone()

GsfBetheHeitlerUpdator* GsfBetheHeitlerUpdator::clone ( void  ) const

inlineoverridevirtual

Implements GsfMaterialEffectsUpdator.

Definition at line 55 of file GsfBetheHeitlerUpdator.h.

References GsfBetheHeitlerUpdator().

{ return new GsfBetheHeitlerUpdator(*this); }

compute()

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

correctedFirstMean()

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, WZElectronSkims53X_cff::max, SiStripPI::mean, SiStripPI::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);
}

correctedFirstVar()

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, WZElectronSkims53X_cff::max, GsfBetheHeitlerUpdator::GSContainer::second, theNrComponents, GsfBetheHeitlerUpdator::GSContainer::third, and ALCARECOEcalPhiSym_cff::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);
}

correctWeights()

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, cms::cuda::for(), 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;
}

getMixtureParameters()

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, btvMC_cff::vz, 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];
    }
}

readParameters()

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

private

Read parametrization from file.

Definition at line 65 of file GsfBetheHeitlerUpdator.cc.

References cms::cuda::assert(), MillePedeFileConverter_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);
  }
}

readPolynomial()

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, and eventshapeDQM_cfi::order.

Referenced by readParameters().

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

Member Data Documentation

MaxOrder

constexpr int GsfBetheHeitlerUpdator::MaxOrder = 6

staticprivate

Definition at line 23 of file GsfBetheHeitlerUpdator.h.

Referenced by readParameters().

MaxSize

constexpr int GsfBetheHeitlerUpdator::MaxSize = 6

staticprivate

Definition at line 22 of file GsfBetheHeitlerUpdator.h.

theCorrectionFlag

int GsfBetheHeitlerUpdator::theCorrectionFlag

private

values to be transformed by logistic / exp. function?

Definition at line 87 of file GsfBetheHeitlerUpdator.h.

Referenced by compute(), and GsfBetheHeitlerUpdator().

theNrComponents

int GsfBetheHeitlerUpdator::theNrComponents

private

Definition at line 85 of file GsfBetheHeitlerUpdator.h.

Referenced by compute(), correctedFirstMean(), correctedFirstVar(), correctWeights(), getMixtureParameters(), GsfBetheHeitlerUpdator(), and readParameters().

thePolyMeans

Polynomial GsfBetheHeitlerUpdator::thePolyMeans[MaxSize]

private

parametrisation of weight for each component

Definition at line 90 of file GsfBetheHeitlerUpdator.h.

Referenced by getMixtureParameters(), and readParameters().

thePolyVars

Polynomial GsfBetheHeitlerUpdator::thePolyVars[MaxSize]

private

parametrisation of mean for each componentP

Definition at line 91 of file GsfBetheHeitlerUpdator.h.

Referenced by getMixtureParameters(), and readParameters().

thePolyWeights

Polynomial GsfBetheHeitlerUpdator::thePolyWeights[MaxSize]

private

correction of 1st or 1st&2nd moments

Definition at line 89 of file GsfBetheHeitlerUpdator.h.

Referenced by getMixtureParameters(), and readParameters().

theTransformationCode

int GsfBetheHeitlerUpdator::theTransformationCode

private

number of components used for parameterisation

Definition at line 86 of file GsfBetheHeitlerUpdator.h.

Referenced by getMixtureParameters(), and readParameters().