GsfTracking/interface/GsfBetheHeitlerUpdator.h

0001 #ifndef GsfBetheHeitlerUpdator_h_
0002 #define GsfBetheHeitlerUpdator_h_
0003
0004 #include "TrackingTools/GsfTracking/interface/GsfMaterialEffectsUpdator.h"
0005
0006 #include "DataFormats/Math/interface/approx_exp.h"
0007 #include "DataFormats/Math/interface/approx_log.h"
0008
0009 #include <iosfwd>
0010 #include <string>
0011
0012 /** \class GsfBetheHeitlerUpdator
0013  *  Description of electron energy loss according to Bethe-Heitler
0014  *  as a sum of Gaussian components. The weights and parameters of the
0015  *  Gaussians as a function of x/X0 are parametrized as polynomials.
0016  *  The coefficients of these polynomials are read from a file at
0017  * construction time.
0018  */
0019
0020 class GsfBetheHeitlerUpdator final : public GsfMaterialEffectsUpdator {
0021 private:
0022   static constexpr int MaxSize = 6;
0023   static constexpr int MaxOrder = 6;
0024
0025   /** Helper class for construction & evaluation of a polynomial
0026    */
0027   class Polynomial {
0028   public:
0029     /// Default constructor (needed for construction of a vector)
0030     Polynomial() {}
0031     /** Constructor from a vector of coefficients
0032      *  (in decreasing order of powers of x)
0033      */
0034     Polynomial(float coefficients[], int is) : m_size(is) {
0035       for (int i = 0; i != m_size; ++i)
0036         theCoeffs[i] = coefficients[i];
0037     }
0038     /// Evaluation of the polynomial
0039     float operator()(float x) const {
0040       float sum = theCoeffs[0];
0041       for (int i = 1; i != m_size; ++i)
0042         sum = x * sum + theCoeffs[i];
0043       return sum;
0044     }
0045
0046   private:
0047     float theCoeffs[MaxOrder] = {0};
0048     int m_size = 0;
0049   };
0050
0051 public:
0052   enum CorrectionFlag { NoCorrection = 0, MeanCorrection = 1, FullCorrection = 2 };
0053
0054 public:
0055   GsfBetheHeitlerUpdator *clone() const override { return new GsfBetheHeitlerUpdator(*this); }
0056
0057 public:
0058   /// constructor with explicit filename and correction flag
0059   GsfBetheHeitlerUpdator(const std::string fileName, const int correctionFlag);
0060
0061 private:
0062   struct GSContainer {
0063     float *first, *second, *third;
0064   };
0065
0066   /// Computation: generates vectors of weights, means and standard deviations
0067   void compute(const TrajectoryStateOnSurface &, const PropagationDirection, Effect[]) const override;
0068
0069 private:
0070   /// Read parametrization from file
0071   void readParameters(const std::string);
0072   /// Read coefficients of one polynomial from file
0073   Polynomial readPolynomial(std::ifstream &, const unsigned int);
0074
0075   /// Filling of mixture (in terms of z=E/E0)
0076   void getMixtureParameters(const float, GSContainer &) const;
0077   /// Correction for weight of component 1
0078   void correctWeights(GSContainer &) const;
0079   /// Correction for mean of component 1
0080   float correctedFirstMean(const float, const GSContainer &) const;
0081   /// Correction for variance of component 1
0082   float correctedFirstVar(const float, const GSContainer &) const;
0083
0084 private:
0085   int theNrComponents;        /// number of components used for parameterisation
0086   int theTransformationCode;  /// values to be transformed by logistic / exp. function?
0087   int theCorrectionFlag;      /// correction of 1st or 1st&2nd moments
0088
0089   Polynomial thePolyWeights[MaxSize];  /// parametrisation of weight for each component
0090   Polynomial thePolyMeans[MaxSize];    /// parametrisation of mean for each componentP
0091   Polynomial thePolyVars[MaxSize];     /// parametrisation of variance for each component
0092 };
0093
0094 #endif