SdReconstruction/LDFFinderKG/EnergyConversion.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_EnergyConversion_h_
#define _LDFFinderKG_EnergyConversion_h_

#include <utl/Math.h>
#include <vector>
#include <cmath>


namespace LDFFinderKG {

  struct EnergyConversion {

    static const unsigned int fgOrderInTheta = 3; // cubic: up to and including x^3
    static const unsigned int fgOrderInS38 = 2; // quadratic: up to and including log10(s38)^2

    double fCicReferenceAngle = 0;
    double fCicReferenceS38 = 0;
    std::pair<double, double> fCicS38Range = { 0, 0 };
    // [ (a0, a1, a2), (b0, b1, b2), (c0, c1, c2) ]
    std::vector<double> fCicParameters[fgOrderInTheta];

    double fEnergyConstant = 0;
    double fEnergySlope = 0;

    //

    void SetCICParameters(const std::vector<double>& par0,
                          const std::vector<double>& par1,
                          const std::vector<double>& par2);

    double GetS38(const double s1000, const double cosTheta) const;

    void GetEnergy(const double cosTheta, const double s1000, 
                   const double s1000Err, const double s1000Sys,
                   double& energy, double& energyErr, double& energySys) const;

  private:
    double GetX(const double cosTheta) const
    { return utl::Sqr(cosTheta) - utl::Sqr(cos(fCicReferenceAngle)); }

    double
    AttenuationFunction(const double x, const double s38)
      const
    {
      // limit S38 to valid range
      const double& s1 = fCicS38Range.first;
      const double& s2 = fCicS38Range.second;
      const double s = (s38 < s1) ?  s1 : ((s38 > s2) ? s2 : s38);
      const auto& c = fCicParameters;
      const double y = std::log10(s / fCicReferenceS38);
      using utl::EvalPoly;
      return EvalPoly({1., EvalPoly(c[0], y), EvalPoly(c[1], y), EvalPoly(c[2], y)}, x);
    }


  };

}


#endif