SdReconstruction/LDFFinderKG/EnergyConversion.cc

Go to the documentation of this file.

#include "EnergyConversion.h"
#include <utl/Math.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <cmath>

using namespace std;
using namespace utl;


namespace LDFFinderKG {

  void
  EnergyConversion::SetCICParameters(const vector<double>& par0,
                                     const vector<double>& par1,
                                     const vector<double>& par2)
  {
    fCicParameters[0] = par0;
    fCicParameters[1] = par1;
    fCicParameters[2] = par2;

    const unsigned int n = fgOrderInS38 + 1;
    for (unsigned int i = 0; i < fgOrderInTheta; ++i)
      if (fCicParameters[i].size() != n)
        throw InvalidConfigurationException("Number of CIC parameters is out of bounds!");
  }


  double
  EnergyConversion::GetS38(const double s1000, const double cosTheta)
    const
  {
    const double eps = 1e-4;
    const double x = GetX(cosTheta);

    double s38 = fCicReferenceS38;
    unsigned int i = 100;
    for (; i; --i) {
      const double s38new = s1000 / AttenuationFunction(x, s38);
      const double diff = std::abs(s38 - s38new);
      s38 = s38new;
      if (s38 <= 0 || diff / s38 < eps)
        break;
    }

    if (!i) {
      ostringstream err;
      err << "S38 assigment did not converge, last value S38 = " << s38;
      ERROR(err);
      return 0;
    }

    if (s38 <= 0) {
      ostringstream warn;
      warn << "Got negative approximation for S38 = " << s38 << " with input values "
              "S1000 = " << s1000 << ", "
              "x = " << x << " (theta = " << acos(cosTheta)/utl::degree << " deg)";
      WARNING(warn);
      return 0;
    }

    return s38;
  }


  void
  EnergyConversion::GetEnergy(const double cosTheta, const double s1000, 
                              const double s1000Err, const double s1000Sys,
                              double& energy, double& energyErr, double& energySys)
    const
  {
    const double s38 = GetS38(s1000, cosTheta);
    if (!s38) {
      energy = energyErr = energySys = 0;
      return;
    }
    energy = fEnergyConstant * std::pow(s38, fEnergySlope);
    const double err = fEnergySlope * energy / s1000;
    energyErr = err * s1000Err;
    energySys = err * s1000Sys;
  }

}