UnivTimeKGGamma.cc

Go to the documentation of this file.

#include "UnivTimeKGConfig.h"

#include <tls/UnivTimeKGGamma.h>

#include <utl/GeometryUtilities.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>

#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>
#include <Math/QuantFuncMathCore.h>

#include <TMath.h>
#include <TF1.h>

#include <string>
#include <cmath>
#include <iostream>
#include <sstream>
#include <algorithm>

using namespace std;
using namespace utl;


namespace UnivTimeKG {

  /*
    Implementation of the time model based on the gumbel distribution
  */

  GammaTimeModel::GammaTimeModel()
    : m(0), s(0), l(0), moff(0), soff(0) { }


  GammaTimeModel::GammaTimeModel(const int i) :
    m(0),
    s(0),
    l(0),
    moff(0),
    soff(0)
  {
    nParams = 2;  // m, s, l of generalized gamma distribution
    icomp = i;

    for (unsigned int ipar = 0; ipar < nParams; ++ipar) {

      ostringstream stringStream;
      string baseDir = DATA_DIR;  // from settings.h
      stringStream << baseDir << "/Gamma/parameters/"
                   << icomp << "_" << ipar << ".xml";

      Reader XMLReader(stringStream.str(), Reader::eNONE);

      Branch topB = XMLReader.GetTopBranch();

      vector<double> rps;
      topB.GetChild("rpoints").GetData(rps);

      Branch parsB = topB.GetChild("parameters");

      fRpoints.push_back(rps);

      vector<Polynomial> ipolys;

      const string ptypes[8] = {"DX0", "DX1", "DX2", "DX3", "Zenith0", "Zenith1", "Energy0", "Energy1"};

      for (unsigned int j = 0; j < sizeof(ptypes) / sizeof(ptypes[0]); ++j) {

        vector<double> coeff;
        parsB.GetChild(ptypes[j]).GetData(coeff);
        ipolys.push_back(Polynomial(coeff));

      }

      fPolys.push_back(ipolys);

    }
  }


  double
  GammaTimeModel::getShapeParameter(const unsigned int /*ipar*/, const vector<double>& /*pars*/, const double /*DX*/)
  {
    return 0;
  }


  void
  GammaTimeModel::CalculateModel(const double DX, const double r, const double ppsi,
                                 const double theta, const double lgE, vector<double>& output)
  {
    // left-right symmetry in showers
    const double psi = fabs(utl::NormalizeAngleMinusPiPi(ppsi));

    const double r0 = 1e3;
    const double rn = r / r0;
    const double DXn = DX / 750;
    const double lgEn = lgE - 19;

    for (unsigned int ipar = 0; ipar < nParams; ++ipar) {

      vector<Polynomial> ipoly = fPolys[ipar];
      vector<double> pars;

      // const double rmin = *min_element(fRpoints[ipar].begin(), fRpoints[ipar].end());
      // const double rmax = *max_element(fRpoints[ipar].begin(), fRpoints[ipar].end());
      // const double rstep = 50;

      for (unsigned int i = 0; i < ipoly.size(); ++i) {

        double rpolyval;

        // Extrapolate before/after first/last point used to fit polynomial
        // if (r < rmin) {
        //   const double fder = (ipoly[i]((rmin + rstep) / r0) - ipoly[i](rmin / r0)) / rstep;
        //   rpolyval = ipoly[i](rmin / r0) + (r - rmin) / r0 * fder;
        // } else if (r > rmax) {
        //   const double fder = (ipoly[i](rmax / r0) - ipoly[i]((rmax - rstep) / r0)) / rstep;
        //   rpolyval = ipoly[i](rmax / r0) + (r - rmax) / r0 * fder;
        // } else
        //   rpolyval = ipoly[i](rn);

        rpolyval = ipoly[i](rn);

        // debug output
        // cout << "Parameter: " << ipar << ", rmin: " << rmin << ", rmax: " << rmax << ", value: " << rpolyval << endl;

        pars.push_back(rpolyval);

      }

      const double pari = pars[0] + DXn * (pars[1] + DXn * (pars[2] + DXn * pars[3]))
                          + sin(theta) * (pars[4] * cos(psi) + pars[5] * DXn)
                          + lgEn * (pars[6] + DXn * pars[7]);

      output.push_back(pari < 0 ? 0 : pari);

    }
  }


  void
  GammaTimeModel::setShapeParameters(const double DX, const double r, const double psi,
                                     const double theta, const double lgE)
  {
    vector<double> params;
    CalculateModel(DX, r, psi, theta, lgE, params);
    m = params[0] + moff;
    s = params[1] + soff;
    // l = params[2];
    l = 0;
  }


  void
  GammaTimeModel::setParameterOffsets(const double m, const double s)
  {
    moff = m;
    soff = s;
  }


  void
  GammaTimeModel::setShapeParametersDirectly(const double mm, const double ss, const double ll)
  {
    m = mm;
    s = ss;
    l = ll;
  }


  // Implementation of the generalized gamma distribution (log to be numerically more robust)
  // reference: http://reliawiki.org/index.php/The_Generalized_Gamma_Distribution
  double
  GammaTimeModel::pdf(const double t)
  {
    if (t <= 0)
      return 0;

    if (l == 0)
      return ROOT::Math::lognormal_pdf(t, m, s, 0);

    const double lnres = (-TMath::Exp(l * (log(t) - m) / s) / l / l
                          - m / l / s + log(t) * (1 / s / l - 1)
                          + log(l) * (1 - 2 / l / l) - log(s)
                          - TMath::LnGamma(1 / l / l));

    return TMath::Exp(lnres);
  }


  // This is analytically correct due to defintion of TMath::Gamma := 1/Gamma(a) * usual defintion
  // of lower incomplete gamma function
  double
  GammaTimeModel::cdf_analytical(const double* const x, const double* const /*p*/)
  {
    const double t = x[0];

    if (t <= 0)
      return 0;

    if (l == 0)
      return ROOT::Math::lognormal_cdf(t, m, s, 0);

    const double z = (log(t) - m) / s;
    return TMath::Gamma(1 / l / l, 1 / l / l * TMath::Exp(l * z));
  }


  double
  GammaTimeModel::cdf(const double t)
  {
    const double xs[] = { t };
    return cdf_analytical(xs, 0);
  }


  double
  GammaTimeModel::invcdf(const double quantile)
  {
    if (l == 0)
      return ROOT::Math::lognormal_quantile(quantile, m, s);

    const TF1 cdfgamma("gamma_cdf", this, &GammaTimeModel::cdf_analytical, 0., 100000., 0);
    return cdfgamma.GetX(quantile);
  }

}