GsfTracking/src/GsfBetheHeitlerUpdator.cc

0001 #include "TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0004 #include "FWCore/ParameterSet/interface/FileInPath.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include <string>
0008 #include <fstream>
0009 #include <cmath>
0010
0011 #include <cassert>
0012
0013 namespace {
0014   /// Logistic function (needed for transformation of weight and mean)
0015   inline float logisticFunction(const float x) { return 1. / (1. + unsafe_expf<4>(-x)); }
0016   /// First moment of the Bethe-Heitler distribution (in z=E/E0)
0017   inline float BetheHeitlerMean(const float rl) { return unsafe_expf<4>(-rl); }
0018   /// Second moment of the Bethe-Heitler distribution (in z=E/E0)
0019   inline float BetheHeitlerVariance(const float rl) {
0020 #if defined(__clang__) || defined(__INTEL_COMPILER)
0021     const
0022 #else
0023     constexpr
0024 #endif
0025         float l3ol2 = std::log(3.) / std::log(2.);
0026     float mean = BetheHeitlerMean(rl);
0027     return unsafe_expf<4>(-rl * l3ol2) - mean * mean;
0028   }
0029 }  // namespace
0030
0031 /*
0032 namespace {
0033  /// Logistic function (needed for transformation of weight and mean)
0034   inline float logisticFunction (const float x) {return 1.f/(1.f+std::exp(-x));}
0035   /// First moment of the Bethe-Heitler distribution (in z=E/E0)
0036   inline float BetheHeitlerMean (const float rl) {
0037     return std::exp(-rl);
0038   }
0039   /// Second moment of the Bethe-Heitler distribution (in z=E/E0)
0040   inline float BetheHeitlerVariance (const float rl)
0041   {
0042 #if __clang__
0043     const
0044 #else
0045     constexpr
0046 #endif
0047     float l3ol2 = std::log(3.)/std::log(2.);
0048     return std::exp(-rl*l3ol2) -  std::exp(-2*rl);
0049   }
0050 }
0051 */
0052
0053 GsfBetheHeitlerUpdator::GsfBetheHeitlerUpdator(const std::string fileName, const int correctionFlag)
0054     : GsfMaterialEffectsUpdator(0.000511, 6), theNrComponents(0), theCorrectionFlag(correctionFlag) {
0055   if (theCorrectionFlag == 1)
0056     edm::LogInfo("GsfBetheHeitlerUpdator") << "1st moment of mixture will be corrected";
0057   if (theCorrectionFlag >= 2)
0058     edm::LogInfo("GsfBetheHeitlerUpdator") << "1st and 2nd moments of mixture will be corrected";
0059
0060   readParameters(fileName);
0061   assert(theNrComponents <= 6);
0062   resize(theNrComponents);
0063 }
0064
0065 void GsfBetheHeitlerUpdator::readParameters(const std::string fileName) {
0066   std::string name = "TrackingTools/GsfTracking/data/";
0067   name += fileName;
0068
0069   edm::FileInPath parFile(name);
0070   edm::LogInfo("GsfBetheHeitlerUpdator") << "Reading GSF parameterization "
0071                                          << "of Bethe-Heitler energy loss from " << parFile.fullPath();
0072   std::ifstream ifs(parFile.fullPath().c_str());
0073
0074   ifs >> theNrComponents;
0075   int orderP;
0076   ifs >> orderP;
0077   ifs >> theTransformationCode;
0078
0079   assert(orderP < MaxOrder);
0080
0081   for (int ic = 0; ic != theNrComponents; ++ic) {
0082     thePolyWeights[ic] = readPolynomial(ifs, orderP);
0083     thePolyMeans[ic] = readPolynomial(ifs, orderP);
0084     thePolyVars[ic] = readPolynomial(ifs, orderP);
0085   }
0086 }
0087
0088 GsfBetheHeitlerUpdator::Polynomial GsfBetheHeitlerUpdator::readPolynomial(std::ifstream& aStream,
0089                                                                           const unsigned int order) {
0090   float coeffs[order + 1];
0091   for (unsigned int i = 0; i < (order + 1); ++i)
0092     aStream >> coeffs[i];
0093   return Polynomial(coeffs, order + 1);
0094 }
0095
0096 void GsfBetheHeitlerUpdator::compute(const TrajectoryStateOnSurface& TSoS,
0097                                      const PropagationDirection propDir,
0098                                      Effect effects[]) const {
0099   //
0100   // Get surface and check presence of medium properties
0101   //
0102   const Surface& surface = TSoS.surface();
0103   //
0104   // calculate components: first check associated material constants
0105   //
0106   float rl(0.f);
0107   float p(0.f);
0108   if (surface.mediumProperties().isValid()) {
0109     LocalVector pvec = TSoS.localMomentum();
0110     p = pvec.mag();
0111     rl = surface.mediumProperties().radLen() / fabs(pvec.z()) * p;
0112   }
0113   //
0114   // produce multi-state only in case of x/X0>0
0115   //
0116   if (rl > 0.0001f) {
0117     //
0118     // limit x/x0 to valid range for parametrisation
0119     // should be done in a more elegant way ...
0120     //
0121     if (rl < 0.01f)
0122       rl = 0.01f;
0123     if (rl > 0.20f)
0124       rl = 0.20f;
0125
0126     float mixtureData[3][theNrComponents];
0127     GSContainer mixture{mixtureData[0], mixtureData[1], mixtureData[2]};
0128
0129     getMixtureParameters(rl, mixture);
0130     correctWeights(mixture);
0131     if (theCorrectionFlag >= 1)
0132       mixture.second[0] = correctedFirstMean(rl, mixture);
0133     if (theCorrectionFlag >= 2)
0134       mixture.third[0] = correctedFirstVar(rl, mixture);
0135
0136     for (int i = 0; i < theNrComponents; i++) {
0137       float varPinv;
0138       effects[i].weight *= mixture.first[i];
0139       if (propDir == alongMomentum) {
0140         //
0141         // for forward propagation: calculate in p (linear in 1/z=p_inside/p_outside),
0142         // then convert sig(p) to sig(1/p).
0143         //
0144         effects[i].deltaP += p * (mixture.second[i] - 1.f);
0145         //    float f = 1./p/mixture.second[i]/mixture.second[i];
0146         // patch to ensure consistency between for- and backward propagation
0147         float f = 1.f / (p * mixture.second[i]);
0148         varPinv = f * f * mixture.third[i];
0149       } else {
0150         //
0151         // for backward propagation: delta(1/p) is linear in z=p_outside/p_inside
0152         // convert to obtain equivalent delta(p)
0153         //
0154         effects[i].deltaP += p * (1.f / mixture.second[i] - 1.f);
0155         varPinv = mixture.third[i] / (p * p);
0156       }
0157       using namespace materialEffect;
0158       effects[i].deltaCov[elos] += varPinv;
0159     }
0160   }
0161 }
0162 //
0163 // Mixture parameters (in z)
0164 //
0165 void GsfBetheHeitlerUpdator::getMixtureParameters(const float rl, GSContainer& mixture) const {
0166   float weight[theNrComponents], z[theNrComponents], vz[theNrComponents];
0167   for (int i = 0; i < theNrComponents; i++) {
0168     weight[i] = thePolyWeights[i](rl);
0169     z[i] = thePolyMeans[i](rl);
0170     vz[i] = thePolyVars[i](rl);
0171   }
0172   if (theTransformationCode)
0173     for (int i = 0; i < theNrComponents; i++) {
0174       mixture.first[i] = logisticFunction(weight[i]);
0175       mixture.second[i] = logisticFunction(z[i]);
0176       mixture.third[i] = unsafe_expf<4>(vz[i]);
0177       ;
0178     }
0179   else  // theTransformationCode
0180     for (int i = 0; i < theNrComponents; i++) {
0181       mixture.first[i] = weight[i];
0182       mixture.second[i] = z[i];
0183       mixture.third[i] = vz[i] * vz[i];
0184     }
0185 }
0186
0187 //
0188 // Correct weights
0189 //
0190 void GsfBetheHeitlerUpdator::correctWeights(GSContainer& mixture) const {
0191   //
0192   // get sum of weights
0193   //
0194   float wsum(0);
0195   for (int i = 0; i < theNrComponents; i++)
0196     wsum += mixture.first[i];
0197   //
0198   // rescale to obtain 1
0199   //
0200   wsum = 1.f / wsum;
0201   for (int i = 0; i < theNrComponents; i++)
0202     mixture.first[i] *= wsum;
0203 }
0204 //
0205 // Correct means
0206 //
0207 float GsfBetheHeitlerUpdator::correctedFirstMean(const float rl, const GSContainer& mixture) const {
0208   //
0209   // calculate difference true mean - weighted sum
0210   //
0211   float mean = BetheHeitlerMean(rl);
0212   for (int i = 1; i < theNrComponents; i++)
0213     mean -= mixture.first[i] * mixture.second[i];
0214   //
0215   // return corrected mean for first component
0216   //
0217   return std::max(std::min(mean / mixture.first[0], 1.f), 0.f);
0218 }
0219 //
0220 // Correct variances
0221 //
0222 float GsfBetheHeitlerUpdator::correctedFirstVar(const float rl, const GSContainer& mixture) const {
0223   //
0224   // calculate difference true variance - weighted sum
0225   //
0226   float var = BetheHeitlerVariance(rl) + BetheHeitlerMean(rl) * BetheHeitlerMean(rl) -
0227               mixture.first[0] * mixture.second[0] * mixture.second[0];
0228   for (int i = 1; i < theNrComponents; i++)
0229     var -= mixture.first[i] * (mixture.second[i] * mixture.second[i] + mixture.third[i]);
0230   //
0231   // return corrected variance for first component
0232   //
0233   return std::max(var / mixture.first[0], 0.f);
0234 }