GsfBetheHeitlerUpdator.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h"

#include "DataFormats/GeometrySurface/interface/MediumProperties.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <string>
#include <fstream>
#include <cmath>

#include<cassert>

namespace {
  inline float logisticFunction (const float x) {return 1./(1.+unsafe_expf<4>(-x));}
  inline float BetheHeitlerMean (const float rl) {
    return unsafe_expf<4>(-rl);
  }
  inline float BetheHeitlerVariance (const float rl)
  {
#if defined(__clang__) || defined(__INTEL_COMPILER)
    const
#else
    constexpr
#endif
    float l3ol2 = std::log(3.)/std::log(2.);
    float mean = BetheHeitlerMean(rl);
    return unsafe_expf<4>(-rl*l3ol2) -  mean*mean;
  }
}

/*
namespace {
  inline float logisticFunction (const float x) {return 1.f/(1.f+std::exp(-x));}
  inline float BetheHeitlerMean (const float rl) {
    return std::exp(-rl);
  }
  inline float BetheHeitlerVariance (const float rl)
  {
#if __clang__
    const
#else    
    constexpr
#endif
    float l3ol2 = std::log(3.)/std::log(2.);
    return std::exp(-rl*l3ol2) -  std::exp(-2*rl);
  }
}
*/

GsfBetheHeitlerUpdator::GsfBetheHeitlerUpdator(const std::string fileName,
                           const int correctionFlag) :
  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);
}

void GsfBetheHeitlerUpdator::readParameters (const std::string fileName)
{  
  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 unsigned int order) {
  float coeffs[order+1];
  for ( unsigned int i=0; i<(order+1); ++i ) aStream >> coeffs[i];
  return Polynomial(coeffs,order+1);
}

void
GsfBetheHeitlerUpdator::compute (const TrajectoryStateOnSurface& TSoS,
                 const PropagationDirection propDir, Effect effects[]) const 
{
   //
  // 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;
    }
  }
}
//
// Mixture parameters (in z)
//
void 
GsfBetheHeitlerUpdator::getMixtureParameters (const float rl,
                          GSContainer & mixture) const
{
 
  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];
  }
}

//
// Correct weights
//
void
GsfBetheHeitlerUpdator::correctWeights (GSContainer & mixture) const
{
  //
  // 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;
}
//
// Correct means
//
float
GsfBetheHeitlerUpdator::correctedFirstMean (const float rl,
                        const GSContainer & mixture) const
{
  //
  // 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);
}
//
// Correct variances
//
float
GsfBetheHeitlerUpdator::correctedFirstVar (const float rl,
                       const GSContainer & mixture) const
{
  //
  // 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);
}