UnivTimeKGLogNormal.cc

Go to the documentation of this file.

#include "UnivTimeKGConfig.h"

#include <tls/UnivTimeKGLogNormal.h>

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

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

using namespace std;


namespace UnivTimeKG {

  /*
    Implementation of the time model based on the log-normal distribution (GAP 2013-022)
  */

  LogNormalTimeModel::LogNormalTimeModel()
    : m(0), s(0), moff(0), soff(0) { }


  LogNormalTimeModel::LogNormalTimeModel(const int i) :
    m(0),
    s(0),
    moff(0),
    soff(0)
  {
    nParams = 2;
    icomp = i;

    for (int ilnpar = 0; ilnpar <= 1; ++ilnpar) {
      for (int ipar = 0; ipar <= 1; ++ipar) {
        std::ostringstream stringStream;
        std::string baseDir = DATA_DIR;  // from settings.h
        stringStream << baseDir << "/LogNormal/parameters/pars_"
                     << icomp << "_" << ilnpar << "_" << ipar;
        addInterpolationTable(stringStream.str());
      }
    }
  }


  double
  LogNormalTimeModel::getShapeParameter(const unsigned int ipar, const vector<double>& pars, const double DX)
  {
    if (ipar == 0) {
      const double g1 = pars[0];
      const double g2 = pars[1];
      return 1e6 * pow((1 + pow((0.888889 + 0.000222222 * DX), 2.0)), (-0.5 * g2)) * pow((4000.0 + DX), (-g1));
    } else if (ipar == 1) {
      const double a = pars[2];
      const double b = pars[3];
      return a + b * DX;
    } else {
      cout << "bad parameter number!" << endl;
      throw;
    }
  }


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


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


  // for testing
  void
  LogNormalTimeModel::setShapeParametersDirectly(const double mm, const double ss)
  {
    m = mm;
    s = ss;
  }


  double
  LogNormalTimeModel::pdf(const double t)
  {
    if (t > 0)
      return ROOT::Math::lognormal_pdf(t, m, s, 0);
    else
      return 0;
  }


  double
  LogNormalTimeModel::cdf(const double t)
  {
    if (t > 0)
      return ROOT::Math::lognormal_cdf(t, m, s, 0);
    else
      return 0;
  }


  double
  LogNormalTimeModel::invcdf(const double quantile)
  {
    return ROOT::Math::lognormal_quantile(quantile, m, s);
  }

}