UnivParamTime.cc

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <fstream>

#include "UnivParamTime.h"
#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>
#include <Math/QuantFuncMathCore.h>
#include "UnivParamTimeConfig.h"
#include <tls/Atmosphere.h>

using namespace UnivParamTimeNS;

UnivParamTime::UnivParamTime(int DetectorType)
{

  if (DetectorType < 0 || DetectorType > 2) {
    printf("UnivParamTime: Wrong Detector Type %d \n", DetectorType);
    exit(1);
  }
  fDetectorType = DetectorType;
  fOffsetM_Mu = 0.;
  fUseThetaInterpolation = false;

  // Read (m,s) table
  {
    std::string DAT = DATA_DIR; // from UnivParamTimeConfig.h
    DAT += "/";
    if (DetectorType == 0)           DAT += "TimeModel-WCD.dat";
    else if (DetectorType == 1)     DAT += "TimeModel-Scin.dat";
    else if (DetectorType == 2)     DAT += "TimeModel-MD.dat";

    FILE* fp = fopen(DAT.c_str(), "r");

    int ok = fscanf(fp, "%d %d  %d %d   %lf %lf ", &ith_min, &ith_max, &ir_min, &ir_max, &logE_ref, &Xmax_ref);
    if (ok != 6) {
      printf("ERROR Reading UnivParamTime parameters\n");
      exit(1);
    }
    if (ir_min < 0 || ir_min >= nParDist || ir_min > ir_max || ir_max < 0 || ir_max >= nParDist ||
        ith_min < 0 || ith_min >= nParTheta || ith_min > ith_max || ith_max < 0 || ith_max >= nParTheta) {
      printf("ERROR Reading UnivParamTime parameters\n");
      exit(1);
    }


    for (int icomp = 0; icomp < 4; icomp++)
      for (int iMS = 0; iMS < 2; iMS++)
        for (int ir = 0; ir < nParDist; ir++) {
          if (icomp > 0 && DetectorType == 2) continue;

          int nread = fscanf(fp, "%lf %lf %lf %lf ", &TimeModelPar_Stage1[icomp][ir][iMS][0], &TimeModelPar_Stage1[icomp][ir][iMS][1],
                             &TimeModelPar_Stage1[icomp][ir][iMS][2], &TimeModelPar_Stage1[icomp][ir][iMS][3]);
          if (nread != nPar_Stage1)  {
            printf("Error reading parameters 1 %d \n", nread);
            exit(1);
          }


          nread = fscanf(fp, "%lf %lf %lf %lf %lf %lf ",
                         &TimeModelPar_Stage2_Theta[icomp][ir][iMS][0][0], &TimeModelPar_Stage2_Theta[icomp][ir][iMS][0][1],
                         &TimeModelPar_Stage2_Theta[icomp][ir][iMS][1][0], &TimeModelPar_Stage2_Theta[icomp][ir][iMS][1][1],
                         &TimeModelPar_Stage2_Theta[icomp][ir][iMS][2][0], &TimeModelPar_Stage2_Theta[icomp][ir][iMS][2][1]);
          if (nread != (nPar_Stage2 * nPar_Stage2_Theta))  {
            printf("Error reading parameters 1 %d \n", nread);
            exit(1);
          }

          for (int ith = 0; ith < nParTheta; ith++) {
            int nread = fscanf(fp, "%lf %lf %lf ", &TimeModelPar_Stage2[icomp][ir][ith][iMS][0], &TimeModelPar_Stage2[icomp][ir][ith][iMS][1],
                               &TimeModelPar_Stage2[icomp][ir][ith][iMS][2]);
            if (nread != nPar_Stage2)  {
              printf("Error reading parameters 2 %d\n", nread);
              exit(1);
            }

            nread = fscanf(fp, "%lf %lf %lf ", &eTimeModelPar_Stage2[icomp][ir][ith][iMS][0], &eTimeModelPar_Stage2[icomp][ir][ith][iMS][1],
                           &eTimeModelPar_Stage2[icomp][ir][ith][iMS][2]);
            if (nread != nPar_Stage2)  {
              printf("Error reading parameters 3 %d \n", nread);
              exit(1);
            }
          }
        }
    fclose(fp);
  }

  // Read Binomial normalization
  {
    std::string DAT = DATA_DIR; // from UnivParamTimeConfig.h
    DAT += "/";
    DAT += "Binomial.norm";

    FILE* f = fopen(DAT.c_str(), "r");
    int iN, ifrac;
    double N, frac, val;
    double fq[3] = {0., 0., 0.};
    while (fscanf(f, "%d %d %lf %lf %lf %lf %lf %lf", &iN, &ifrac, &N, &frac, &val, &fq[0], &fq[1], &fq[2]) != EOF) {
      if (iN < 0 || iN > 199 || ifrac < 0 || ifrac > 98) {
        printf("Error reading Binomial coefficients \n");
        exit(1);
      }
      BinomialNorm[iN][ifrac] = val;
      Binomial_q[iN][ifrac][0] = fq[0];
      Binomial_q[iN][ifrac][1] = fq[1];
      Binomial_q[iN][ifrac][2] = fq[2];
    }
    fclose(f);
  }
}



double UnivParamTime::GetMS_ir_ith(double DL, double logE, int      ith, int   ir, double psi , int icomp, int iMS)
{
  double r    = TimeModelParDist[ir];
  double theta = TimeModelParTheta[ith];

  double* par_stage1 = TimeModelPar_Stage1[icomp][ir][iMS];
  double* par_stage2 = TimeModelPar_Stage2[icomp][ir][ith][iMS];

  return GetMS(DL, logE,  theta,  r,  psi ,  icomp,  iMS,  par_stage1, par_stage2);
}


double UnivParamTime::GetMS_ir(double DL, double logE, double theta, int   ir, double psi , int icomp, int iMS)
{
  double r = TimeModelParDist[ir];
  double theta_deg = theta * 180. / M_PI;
  if (theta_deg < theta_Low) theta_deg = theta_Low ;
  if (theta_deg > theta_High) theta_deg = theta_High;

  double* par_stage1 = TimeModelPar_Stage1[icomp][ir][iMS];
  double par_stage2[nPar_Stage2];

  for (int ipar = 0; ipar < nPar_Stage2; ipar++) {
    double A = TimeModelPar_Stage2_Theta[icomp][ir][iMS][ipar][0];
    double B = TimeModelPar_Stage2_Theta[icomp][ir][iMS][ipar][1];

    if (ipar > 0) par_stage2[ipar] = sin(theta_deg * M_PI / 180. / B) * A;
    else          par_stage2[ipar] = A + (B - A) * theta_deg / 60.;
  }

  return GetMS(DL, logE,  theta,  r,  psi ,  icomp,  iMS,  par_stage1, par_stage2);

}


double
UnivParamTime::GetMS(double DL, double logE, double theta, double r, double psi , int icomp, int /*iMS*/, double* par_stage1, double* par_stage2)
{
  double ms0 = par_stage1[0];
  double ms1 = par_stage1[1];
  double ms2 = par_stage1[2];
  double ms_logE = par_stage1[3];

  double ms_slope = par_stage2[0];
  double ms_psi_down = par_stage2[1];
  double ms_psi_up  = par_stage2[2];

  //  DL dependence at Xmax=Xmax_ref using psi=90 [deg] and different zenith angles
  double val;

  double val0 = ms0;
  double val1 = ms0 - ms1;
  double B = ms2;

  double C = (val1 - val0) / (exp(-B) - 1.0);
  double A = val0 - C;

  val = A + C * exp(-B * DL / 1000.);

  // Local modification of the DL dependence for psi=90 [deg]
  AtmosphereNS::Atmosphere fatm;
  double DL_ref = fatm.Get_DX_DL(r,  psi,  Xmax_ref     , theta,  AtmosphereNS::hGroundRef,
                                 UnivParamTimeNS::UseDL_time[icomp], UnivParamTimeNS::UseDiffusive_time[icomp]);  //Avg over the 12 atmospheres
  val +=  ms_slope * (DL - DL_ref) / 1000.;

  // Psi modulation
  double ms_psi = (cos(psi) > 0. ? ms_psi_up * cos(psi) : ms_psi_down * cos(psi));
  val += ms_psi;

  //Energy depedence
  val += ((logE - logE_ref) * ms_logE);

  return val;
}


double UnivParamTime::interpol(double v1, double v2, double x1, double x2, double x)
{
  if (x1 == x2) return v1;
  return v1 + (v2 - v1) * (x - x1) / (x2 - x1);
}


double UnivParamTime::GetMS(double DL, double logE, double theta, double r, double psi , int icomp, int iMS)
{
  if (r < r_Low)         r = r_Low;
  if (r > r_High)         r = r_High;
  if (logE < logE_Low) logE = logE_Low;
  if (logE > logE_High)  logE = logE_High;

  double val = 0.;

  //--------
  int ir1 = 0, ir2 = 0;
  int ir_max_c = ir_max;
  if (fUseThetaInterpolation && (icomp == 1 || icomp == 3) && theta > 55.*M_PI / 180. && ir_max >= 5) ir_max_c = 5; // WARNING : exception

  if (r >= TimeModelParDist[ir_max_c])     {
    ir1 = ir_max_c - 1;
    ir2 = ir_max_c;
  } else if (r <= TimeModelParDist[ir_min + 1]) {
    ir1 = ir_min;
    ir2 = ir_min + 1;
  } else {
    ir1 = ir_min;
    while (r > TimeModelParDist[ir2]) ir2++;
    ir1 = ir2 - 1;
  }

  //-------------------------------
  //--- Interpolation in (r,theta)
  //-------------------------------
  if (fUseThetaInterpolation) {
    int ith1 = 0, ith2 = 0;
    if (theta >= TimeModelParTheta[ith_max])   {
      ith1 = ith_max - 1;
      ith2 = ith_max;
    } else if (r <= TimeModelParTheta[ith_min + 1])   {
      ith1 = ith_min;
      ith2 = ith_min + 1;
    } else {
      ith1 = ith_min;
      while (theta > TimeModelParTheta[ith2]) ith2++;
      ith1 = ith2 - 1;
    }

    double val_ith1[2] = { GetMS_ir_ith(DL, logE, ith1, ir1, psi , icomp, iMS), GetMS_ir_ith(DL, logE, ith1, ir2, psi , icomp, iMS) };
    double val0 = interpol(val_ith1[0], val_ith1[1], TimeModelParDist[ir1], TimeModelParDist[ir2], r);

    double val_ith2[2] = { GetMS_ir_ith(DL, logE, ith2, ir1, psi , icomp, iMS), GetMS_ir_ith(DL, logE, ith2, ir2, psi , icomp, iMS) };
    double val1 = interpol(val_ith2[0], val_ith2[1], TimeModelParDist[ir1], TimeModelParDist[ir2], r);

    val = interpol(val0, val1, sin(TimeModelParTheta[ith1]), sin(TimeModelParTheta[ith2]), sin(theta));
  }
  //-------------------------------
  //--- Interpolation in r
  //-------------------------------
  else {
    double val_i[2] = { GetMS_ir(DL, logE, theta, ir1, psi , icomp, iMS), GetMS_ir(DL, logE, theta, ir2, psi , icomp, iMS) };
    val = interpol(val_i[0], val_i[1], TimeModelParDist[ir1], TimeModelParDist[ir2], r);
  }


  //--- Offset in <m> for Muon/Em Muon component
  if (iMS == 0) val += GetOffsetM_r(r, icomp);

  if (val < ms_minVal) val = ms_minVal;

  return val;

}


bool UnivParamTime::GetShapeParameters(double* DL, double r, double logE, double psi, double theta, double* m, double* s)
{
  for (int iMS = 0; iMS < 2; iMS++)
    for (int icomp = 0; icomp < 4; icomp++) {
      double val = GetMS(DL[icomp], logE, theta, r ,  psi,  icomp, iMS);
      if (iMS == 0) m[icomp] = val;
      else          s[icomp] = val;
    }
  return true;
}



//------------------------------------------------------------------
//----  CDF/PDF
//------------------------------------------------------------------


double UnivParamTime::GetFraction(double* t0, double* m, double* s, double* fcomp, double ti, bool isCDF)
{
  double frac = 0., sum = 0.;
  bool ok = true;
  for (int icomp = 0; icomp < 4; icomp++) {
    if (fcomp[icomp] < 0.) {
      ok = false;
      continue;
    }
    if (fcomp[icomp] > 0. && (m[icomp] <= 0. ||  s[icomp] <= 0.)) {
      ok = false;
      continue;
    }
    sum += fcomp[icomp];
  }
  if (sum == 0. || !ok) return UNDEF;

  for (int icomp = 0; icomp < 4; icomp++) {
    if (fcomp[icomp] == 0.) continue;
    if (ti < t0[icomp])     continue;

    double frac_i;
    if (isCDF) frac_i = ROOT::Math::lognormal_cdf(ti - t0[icomp], m[icomp], s[icomp]) * (fcomp[icomp] / sum);
    else        frac_i = ROOT::Math::lognormal_pdf(ti - t0[icomp], m[icomp], s[icomp]) * (fcomp[icomp] / sum);
    frac += frac_i;
  }
  return frac;
}

double UnivParamTime::GetFraction(double* DL, double r ,  double logE, double psi, double theta, double* fcomp, double* t0, double ti, bool isCDF)
{

  double m[4], s[4];
  bool ok = GetShapeParameters(DL, r, logE, psi, theta, m, s);
  if (!ok) return UNDEF;

  return GetFraction(t0 , m,  s, fcomp, ti, isCDF);
}

//-------

double UnivParamTime::GetCDF(double* t0, double* m, double* s, double* fcomp, double ti)
{
  return GetFraction(t0, m, s, fcomp,  ti, true);
}

double UnivParamTime::GetCDF(double* DL, double r ,  double logE, double psi, double theta, double* fcomp, double* t0, double ti)
{
  return GetFraction(DL, r, logE, psi, theta, fcomp, t0, ti, true);
}

//-------

double UnivParamTime::GetPDF(double* t0, double* m, double* s, double* fcomp,  double ti)
{
  return GetFraction(t0, m, s, fcomp,  ti, false);
}
double UnivParamTime::GetPDF(double* DL, double r ,  double logE, double psi, double theta, double* fcomp, double* t0, double ti)
{
  return GetFraction(DL, r, logE, psi, theta, fcomp, t0, ti, false);
}

//------------------------------------------------------------------
//----  Time quantiles
//------------------------------------------------------------------

double UnivParamTime::GetTime(double* DL, double r ,  double logE, double psi, double theta, double* fcomp, double* t0,  double fi)
{
  double m[4], s[4];

  bool ok = true;
  for (int icomp = 0; icomp < 4; icomp++)
    if (isnan(DL[icomp]) || isinf(DL[icomp]) || isnan(fcomp[icomp]) || isinf(fcomp[icomp]) || isnan(t0[icomp]) || isinf(t0[icomp]) || isnan(fi)) ok = false;
  if (!ok) return 0.;
  ok = GetShapeParameters(DL, r, logE, psi, theta, m, s);
  if (!ok) return 0.;

  return GetTime(t0 , m,  s, fcomp,  fi);
}


double UnivParamTime::GetTime(double* t0, double* m, double* s, double* fcomp, double fi)
{
  // NOTE: the standard algorithm would be to increase in one direction with a fixed step and once the
  //       we are over the desire fraction, go back and decrease the step and continue upwards.
  //       However, as it is it seems to be faster

  bool ok = true;
  for (int icomp = 0; icomp < 4; icomp++)
    if (isnan(t0[icomp]) || isinf(t0[icomp]) || isnan(m[icomp]) || isinf(m[icomp]) || isnan(s[icomp]) || isinf(s[icomp]) || isnan(fi)
        || isnan(fcomp[icomp]) || isinf(fcomp[icomp]) || fi < 0 || fi > 1.) ok = false;
  if (!ok) return 0.;


  double Tolerance = 1.e-8;

  //----

  double tmin = 1.e10, tmax = 0.;
  for (int icomp = 0; icomp < 4; icomp++) {
    if (fcomp[icomp] == 0.) continue;

    double tmin_i = t0[icomp];
    double tmax_i = t0[icomp] + exp(m[icomp] + 5.*s[icomp]);

    if (tmin > tmin_i) tmin = tmin_i;
    if (tmax < tmax_i) tmax = tmax_i;
  }
  if (tmin < 1.e-3) tmin = 1.e-3;

  //-----

  double xmin = log(tmin);
  double xmax = log(tmax);

  double fmin = GetCDF(t0 , m,  s, fcomp,  exp(xmin));
  double fmax = GetCDF(t0 , m,  s, fcomp,  exp(xmax));

  if (fmin >= fi) return exp(xmin);
  if (fmax <= fi) return exp(xmax);

  double x = (xmin + xmax) / 2.;
  double step = (xmax - xmin) / 2.;
  double sign = 1.;
  int iter = 0;
  const int nIterMax = 1000;
  while (true) {
    double frac = GetCDF(t0 , m,  s, fcomp, exp(x));
    if (fabs(frac - fi) < Tolerance) {
      double a[2] = {x, x + sign * step};
      double b[2] = {frac, GetCDF(t0 , m,  s, fcomp, exp(x + sign * step))};
      x = a[0] + (a[1] - a[0]) / (b[1] - b[0]) * (fi - b[0]);
      break;
    }
    step /= 2.;
    if (frac > fi)   {
      sign = -1.;
      x -= step;
    } else           {
      sign = 1.;
      x += step;
    }
    iter++;
    if (iter > nIterMax) break;
  }

  if (iter > nIterMax) {
    printf("ERROR: time quantile could not be calculated tq=%lf fi=%lf (%lf)  | %lf %lf \n", exp(x), fi, GetCDF(t0 , m,  s, fcomp, exp(x)), fmin, fmax);
    printf(" t0=(%lf,%lf,%lf,%lf) m=(%lf,%lf,%lf,%lf) s=(%lf,%lf,%lf,%lf) fcomp=(%lf,%lf,%lf,%lf) \n",
           t0[0], t0[1], t0[2], t0[3],
           m[0], m[1], m[2], m[3],
           s[0], s[1], s[2], s[3],
           fcomp[0], fcomp[1], fcomp[2], fcomp[3]);
    return 0.;
  }
  if (isnan(x) || isinf(x)) return 0.;

  return exp(x);
}


//------------------------------------------------------------------
//----  Time quantile PDF
//          Use f=0.01,0.02,...0.99
//          The model is only valid for integer (npart,N)
//------------------------------------------------------------------

double UnivParamTime::GetPoissonFactor(double fem, double f)
{
  if (fDetectorType == 2) return 1.;

  if (fem < 0.) fem = 0.;
  if (fem > 1.) fem = 1.;

  // double A = 0.15;
  double B = 0.93;
  B = 1. + (f - 0.5);
  double pf = 1. + B * pow(fem / 0.5, 2.);
  if (pf < 1.) pf = 1.;

  return pf;
}

double UnivParamTime::tQuantilePDF(double* t0, double* m, double* s, double* fcomp, double vemTot, double ti, double f,
                                   bool UseApprox, double pf, double& mean, double& rms)
{
  double pdf = 1.e-100;
  mean = UNDEF;
  rms = UNDEF;
  double N = vemTot * pf;
  if (isnan(N)) return pdf;

  //----------------------------------------------------------------
  //    Approximation (so mean/rms can be calculated)
  //----------------------------------------------------------------

  // Quantiles of P(fi)
  double fq[3];
  QuantilePDF(N , -100. ,  f,  fq);

  // Quantiles of P(ti)
  double tq[3];
  for (int iq = 0; iq < 3; iq++)
    tq[iq] = GetTime(t0, m, s, fcomp, fq[iq]);

  const double ndev = 1.28156;
  // Normal approximation
  if ((tq[2] - tq[1]) < (tq[1] - tq[0])) {
    mean = tq[1];
    rms = (tq[1] - tq[0]) / ndev;
    if (UseApprox) pdf = ROOT::Math::normal_pdf(ti, rms, mean);
  }
  // Log normal
  else {
    double t0_i = tq[1]  - (tq[2] - tq[1]) * (tq[1] - tq[0]) / (tq[2] + tq[0] - 2.*tq[1]);
    double m_i = log(tq[1] - t0_i);
    double s_i = (log(tq[1] - t0_i) - log(tq[0] - t0_i)) / ndev;

    mean = exp(m_i + s_i * s_i / 2.) + t0_i;
    rms = sqrt(exp(s_i * s_i + 2.*m_i) * (exp(s_i * s_i) - 1.));

    if (UseApprox && ti > t0_i) pdf = ROOT::Math::lognormal_pdf(ti - t0_i, m_i, s_i);
  }

  //----------------------------------------------------------------
  //
  //----------------------------------------------------------------
  if (!UseApprox) {
    double fi = GetCDF(t0, m, s, fcomp, ti);
    if (fi > 0.) {
      double derivative = (fi - GetCDF(t0, m, s, fcomp, ti * 0.99)) / ti / 0.01;
      pdf = QuantilePDF(N, fi, f);
      pdf *= derivative;
    }
  }

  return pdf;
}

double UnivParamTime::tQuantilePDF(double* t0, double* m, double* s, double* fcomp, double vemTot, double ti, double f,
                                   bool UseApprox, double& mean, double& rms)
{
  double fem0 = fcomp[1] / (fcomp[0] + fcomp[1] + fcomp[2] + fcomp[3]);
  double pf = GetPoissonFactor(fem0, f);
  return tQuantilePDF(t0, m, s, fcomp, vemTot, ti, f, UseApprox, pf, mean, rms);
}

double UnivParamTime::tQuantilePDF(double* DL, double r, double logE, double psi, double theta, double* fcomp, double* t0, double vemTot, double ti, double f,
                                   bool UseApprox, double& mean, double& rms)
{
  mean = 0., rms = 0.;
  double m[4], s[4];
  bool ok = GetShapeParameters(DL, r, logE, psi, theta, m, s);
  if (!ok) return UNDEF;

  double pdf = tQuantilePDF(t0, m, s, fcomp, vemTot, ti, f, UseApprox, mean, rms);
  return pdf;
}


//----------------------------------------

double UnivParamTime::GetBinomialNorm(double N, double f , double* fq)
{
  double stepN = log10(40000. / 100.) / 100.;


  // iN=0,......,98 |  99,....., 199
  //  N=1,......,99 | 100,     , 40000
  double N1, N2;
  int iN;
  if (N < 99.5) {
    iN = int(N) - 1;
    if (iN < 0) iN = 0;
    N1 = double(iN + 1);
    N2 = double(iN + 2);
  } else {
    iN = int(log10(N / 100.) / stepN + 0.5);
    if (iN > 99) iN = 99;
    N1 = pow(10., double(iN) * stepN + 2.);
    N2 = pow(10., double(iN + 1) * stepN + 2.);
    iN += 99;
  }


  // (ifrac)
  f = double(int(f * 100)) / 100.;
  double stepFrac = 0.01;
  int ifrac = int((f - stepFrac) / stepFrac + 0.5);
  if (ifrac < 0) ifrac = 0;
  if (ifrac > 98) ifrac = 98;

  double norm = interpol(BinomialNorm[iN][ifrac], BinomialNorm[iN + 1][ifrac], N1, N2, N);
  for (int iq = 0; iq < 3; iq++)
    fq[iq] = interpol(Binomial_q[iN][ifrac][iq], Binomial_q[iN + 1][ifrac][iq], N1, N2, N);

  return norm;
}


double UnivParamTime::QuantilePDF(double N , double fi, double f)
{
  double fq[3] = {0., 0., 0.};
  return QuantilePDF(N, fi, f, fq);
}


double UnivParamTime::QuantilePDF(double N , double fi, double f, double* fq)
{
  // Do not extrapolate
  if (N > 40000.) N = 40000.;

  // Change N if ipos<0 so <f> is reproduced.
  if (f * (N + 1) - 1 < 0.) N = 1. / f - 1.;
  double ipos = f * (N + 1) - 1;

  double norm = GetBinomialNorm(N, f, fq);

  if (fi < 0. || fi > 1.) return 0.;

  double logval = (N - ipos - 1) * log((1. - fi) / (1. - f)) + (ipos) * log(fi / f);
  double val = logval - log(norm);
  if (val < -100.) val = -100.;
  return exp(val);
}

//------------------------------------------------------------------
//----  PDF of first particle time
//------------------------------------------------------------------

double UnivParamTime::tFirstPDF(double DL_mu, double r, double logE, double psi, double theta, double t0_mu, double npart, double t)
{
  double mean = 0., rms = 0.;
  return tFirstPDF(DL_mu, r, logE, psi, theta, t0_mu, npart,  t, false, mean, rms);
}

double UnivParamTime::tFirstPDF(double DL_mu, double r, double logE, double psi, double theta, double t0_mu, double npart, double t,
                                bool UseApprox, double& mean, double& rms)
{
  mean = 0., rms = 0.;

  int icomp = 0;
  double m_mu = GetMS(DL_mu, logE, theta, r , psi,   icomp, 0);
  double s_mu = GetMS(DL_mu, logE, theta, r,  psi,   icomp, 1);

  if (npart <= 1.) npart = 1.;
  double prob = tFirstPDF(t0_mu, m_mu, s_mu, npart, t, UseApprox, mean, rms);

  return prob;
}

double UnivParamTime::tFirstPDF(double t0_mu, double m_mu, double s_mu, double npart, double ti, bool UseApprox, double& mean, double& rms)
{
  double t0[4] = {t0_mu, 0., 0., 0.};
  double m[4] = {m_mu, 0., 0., 0.};
  double s[4] = {s_mu, 0., 0., 0.};
  double fcomp[4] = {1., 0., 0., 0.};

  double tFirstPDF = 1.e-100;
  mean = UNDEF;
  rms = UNDEF;

  //----------------------------------------------------------------
  //    Approximation (so mean/rms can be calculated)
  //----------------------------------------------------------------

  //-----
  double q1[3] = {0.1, 0.5, 0.9};
  double t1[3];
  for (int i = 0; i < 3; i++) {
    double q = 1. - pow(1. - q1[i], 1. / npart);
    t1[i] = GetTime(t0, m, s, fcomp, q);
  }

  //-----

  const double ndev = 1.28156;
  // Normal
  if ((t1[2] - t1[1]) < (t1[1] - t1[0])) {
    mean = t1[1];
    rms = (t1[1] - t1[0]) / ndev;
    if (UseApprox) tFirstPDF = ROOT::Math::normal_pdf(ti, rms, mean);
  }
  // Log normal
  else {
    double t0_i = t1[1]  - (t1[2] - t1[1]) * (t1[1] - t1[0]) / (t1[2] + t1[0] - 2.*t1[1]);
    double m_i = log(t1[1] - t0_i);
    double s_i = (log(t1[1] - t0_i) - log(t1[0] - t0_i)) / ndev;

    mean = exp(m_i + s_i * s_i / 2.) + t0_i;
    rms = sqrt(exp(s_i * s_i + 2.*m_i) * (exp(s_i * s_i) - 1.));
    if (UseApprox && ti > t0_i) tFirstPDF = ROOT::Math::lognormal_pdf(ti - t0_i, m_i, s_i);
  }


  //----------------------------------------------------------------
  //
  //----------------------------------------------------------------
  if (!UseApprox && ti > t0[0]) {
    double cdf = GetCDF(t0, m, s, fcomp,  ti);
    double pdf = GetPDF(t0, m, s, fcomp,  ti);
    tFirstPDF = npart * pow((1 - cdf), npart - 1) * pdf;
  }

  return tFirstPDF;

}


//---------------------------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------------------------------
//
// [iq] = [0] T1vem  [1] T10  [2] T50
//

//-----  Corrections due failure of : (a) Log normal ansatz (b) Time variance model shifts (c) Discrete pulse effects
//----  This quantity must be substracted to the reconstructed quantile
//   From MultiShower/time-corrections.cc
//       [NanoSecPerVEM] =   (T50-T10)/vem/0.4  [ns]/[vem]
double UnivParamTime::tQuantileCorrection(double NanoSecPerVEM, int iq)
{
  if (fDetectorType == 2) return 0.;
  if (iq < 0 || iq > 2) return 0.;

  double A[3] = {150.0, 70.0, 0.0};
  double corr = A[iq] * NanoSecPerVEM / 200.;

  return corr;
}

//---------------------------------------------------------------------------------------------------------------------------------
//----  This quantity should be substracted to the reconstructed start time
//               Valid only r<1000 [m].
//  From MultiShower/time-corrections.cc
//---------------------------------------------------------------------------------------------------------------------------------

double UnivParamTime::tStartCorrection(double r, double logE, double theta, bool Is8nsFADC)
{
  if (fDetectorType == 2) return 0.;
  double logE_i[3] = {19.0, 19.5, 20.0};
  int iE;
  if (logE <= logE_i[0])      iE = 0;
  else if (logE <= logE_i[1]) iE = 1;
  else                      iE = 1;

  double A1_25ns[3] = { 0.00, 0.00, 0.00};
  double B1_25ns[3] = {34.65, 29.17, 21.78};
  double A2_25ns[3] = { 0.00, 0.00, 0.00};
  double B2_25ns[3] = { -41.13, -33.07, -26.58};


  double A1_8ns[3] = { 0.00, 0.00, 0.00};
  double B1_8ns[3] = {41.03, 30.93, 23.21};
  double A2_8ns[3] = { 0.00, 0.00, 0.00};
  double B2_8ns[3] = { -48.71, -38.09, -30.47};


  double* A1_i, *A2_i, *B1_i, *B2_i;
  if (Is8nsFADC) {
    A1_i = A1_8ns;
    A2_i = A2_8ns;
    B1_i = B1_8ns;
    B2_i = B2_8ns;
  } else             {
    A1_i = A1_25ns;
    A2_i = A2_25ns;
    B1_i = B1_25ns;
    B2_i = B2_25ns;
  }

  double A1 = interpol(A1_i[iE], A1_i[iE + 1], logE_i[iE], logE_i[iE + 1], logE);
  double B1 = interpol(B1_i[iE], B1_i[iE + 1], logE_i[iE], logE_i[iE + 1], logE);
  double A2 = interpol(A2_i[iE], A2_i[iE + 1], logE_i[iE], logE_i[iE + 1], logE);
  double B2 = interpol(B2_i[iE], B2_i[iE + 1], logE_i[iE], logE_i[iE + 1], logE);

  double sectheta = 1. / cos(theta);
  double A = A1 * (sectheta - 1.) + A2;
  double B = B1 * (sectheta - 1.) + B2;

  if (r < 200.e2) return 0.;
  double corr = A + B * (r - 200.e2) / 800.e2;
  return corr;
}


//---------------------------------------------------------------------------------------------------------------------------------
//----  This quantity should be substracted to the reconstructed quantile to get the expected value for AreaOverPeak=3.60
//---------------------------------------------------------------------------------------------------------------------------------
//
// [iq] = [0] T1vem  [1] T10  [2] T50
// [tq]
//
//   From MultiShower/time-corrections.cc

double UnivParamTime::tQuantileCorrection_AoP(double RiseTime, int iq, double AreaOverPeak)
{
  if (iq < 0 || iq > 2) return 0.;
  if (AreaOverPeak < 0.) return 0.;
  if (fDetectorType > 0) return 0.; // Not implemented

  const double RiseTimePulse = 40.;
  if (RiseTime < RiseTimePulse) RiseTime = RiseTimePulse;

  const double AsymptoticCorr[3] = { 7.5, 13., 17. };

  const double A = AsymptoticCorr[iq] * (AreaOverPeak - 3.60) / (4.30 - 3.60);
  const double B[3] = { 1.0, 1.30 , 4.6 };

  double corr = A * (1. - exp(-B[iq] * RiseTime / 200.));
  return corr;
}


//---------------------------------------------------------------------------------------------------------------------------------
//----  D_TO  (proton QGSJet II)      r=800 [m] psi=90 [deg]   Convolved with the detector response
//---------------------------------------------------------------------------------------------------------------------------------

double UnivParamTime::GetD_TO(double X, double theta, double logE, int icomp)
{
  if (fDetectorType == 2 && icomp > 0) return UNDEF;

  double secTheta = (theta > 60.*M_PI / 180. ?  2. :  1. / cos(theta));

  const double D0_TO_Par_WCD[4][3] = { {  27.16 ,  36.33 ,   0.00   } , {  27.01 ,  -3.59 ,   0.00   } , {  29.54 ,  34.38 ,   0.00   } , {  18.42 ,  -5.34 ,   0.00   } };
  const double D0_TO_Par_Scin[4][3] = { {  17.79 ,  33.85 ,   0.00   } , {  17.95 ,  -11.00 ,   0.00   } , {  23.48 ,  19.71 ,   0.00   } , {  19.07 ,  -12.06 ,   0.00   } };
  const double D0_TO_Par_MD[4][3]  = { {  85.40 ,  46.60 ,   0.00   } , {   0.00 ,   0.00 ,   0.00   } , {   0.00 ,   0.00 ,   0.00   } , {   0.00 ,   0.00 ,   0.00   } };


  const double D1_TO_Par_WCD[4] =  { 309.14, -4.56, 280.30 , -1.68 };
  const double D1_TO_Par_Scin[4] = { 338.33 , 1.98, 707.64 ,  1.39 };
  const double D1_TO_Par_MD[4] =   { 181.82 , 0., 0., 0. };

  const double D2_TO_Par_WCD[4] = {  3.58,  -0.07 ,  5.73,   1.75 };
  const double D2_TO_Par_Scin[4] = {  2.13 , -1.34 ,  2.17 ,  2.56 };
  const double D2_TO_Par_MD[4] =  { 25.65 ,  0.00 ,  0.00 ,  0.00 };

  double D0_TO, D1_TO, D2_TO;
  if (fDetectorType == 0) {
    D0_TO = D0_TO_Par_WCD[icomp][0] + D0_TO_Par_WCD[icomp][1] * (secTheta - 1.) +  D0_TO_Par_WCD[icomp][2] * pow(secTheta - 1, 2.);
    D1_TO = D1_TO_Par_WCD[icomp];
    D2_TO = D2_TO_Par_WCD[icomp];
  } else if (fDetectorType == 1) {
    D0_TO = D0_TO_Par_Scin[icomp][0] + D0_TO_Par_Scin[icomp][1] * (secTheta - 1.) +  D0_TO_Par_Scin[icomp][2] * pow(secTheta - 1, 2.);
    D1_TO = D1_TO_Par_Scin[icomp];
    D2_TO = D2_TO_Par_Scin[icomp];
  } else if (fDetectorType == 2) {
    D0_TO = D0_TO_Par_MD[icomp][0] + D0_TO_Par_MD[icomp][1] * (secTheta - 1.) +  D0_TO_Par_MD[icomp][2] * pow(secTheta - 1, 2.);
    D1_TO = D1_TO_Par_MD[icomp];
    D2_TO = D2_TO_Par_MD[icomp];
  } else return UNDEF;

  double X0 = 750.;
  double D_TO;
  if (icomp == 1 || icomp == 3) D_TO = D0_TO  + (X - X0) / 200.*D1_TO  +  D2_TO * (logE - 19.);
  else                        D_TO = (D0_TO + D2_TO * (logE - 19.)) * exp(-(X - X0) / D1_TO);



  return D_TO;
}

//---------------------------------------------------------------------------------------------------------------------------------
//----
//---------------------------------------------------------------------------------------------------------------------------------

double UnivParamTime::GetOffsetM_r(double r, int icomp)
{
  if (fDetectorType == 2)      return 0.;
  if (icomp == 1 || icomp == 3)  return 0.;

  double SysPar1 = 2.65;
  double C = fOffsetM_Mu / (1. - exp(-SysPar1));
  double shift = C * (1. - exp(-SysPar1 * r / 2000.e2));

  return shift;
}