#ifndef __UnivParamTime_h_
#define __UnivParamTime_h_

#include <math.h>

#ifndef UNDEF
#define UNDEF -100
#endif

//----------------------------------------------------------
//  WARNINGS:
//          GetAtmosphericCorrection/ tQuantileCorrection / tStartCorrection : calculated for WCD and used for ASCII
//          tQuantileCorrection_AoP                                          : not implemented for ASCII
//
//----------------------------------------------------------

namespace UnivParamTimeNS {


  //  Variables needed to calculate attenuation parameter
  //-----------------------------------------------------

  const bool     XmaxShift_time[4]   = {    0.,   0. ,    0.,   0.  };
  const bool     UseDiffusive_time[4] = {  true, true ,  true, true  };
  const bool     UseDL_time[4]       = {  true, true ,  true, true  };

  //  Parameterization distances/angles
  //-----------------------------------------------------

  static const int nParDist = 8;
  static const double TimeModelParDist[nParDist] = {100.e2, 200.e2, 400.e2, 800.e2, 1000.e2, 1500.e2, 2000.e2, 2500.e2};
  static const int nParTheta = 7;
  static const double TimeModelParTheta[nParTheta] = { 0., 12.*M_PI / 180., 25.*M_PI / 180., 36.*M_PI / 180., 45.*M_PI / 180., 53.*M_PI / 180., 60.*M_PI / 180.};

  const int nPar_Stage1 = 4, nPar_Stage2 = 3, nPar_Stage2_Theta = 2;

  //
  const double ndev = 1.28156; // The number of sigma that the reference quantiles correspond to

  // Ranges outside of which the variable is fixed to the lower/upper range
  //-----------------------------------------------------
  static const double theta_Low = 0., theta_High = 65.*180. / M_PI;
  static const double logE_Low = 17.5, logE_High = 21.;
  static const double r_Low = 10.e2;
  static const double r_High = 2200.e2;
  static const double ms_minVal = 0.01;

  //-------------------------
  //-------------------------
  static const double OffsetM_Mu_0deg [4] = { 0.00,  0.03, 0.04, 0.07 };
  static const double OffsetM_Mu_60deg[4] = { 0.00,  0.02, 0.03, 0.04 };
  //-------------------------
  //-------------------------

  class UnivParamTime {

  private:

    int fDetectorType; // [0] WCD
    double fOffsetM_Mu;
    bool fUseThetaInterpolation;

    int    ith_min, ith_max;                      //  Zenith angle ranges  -> Interpolation limits
    int    ir_min, ir_max;                        //  Core distance ranges -> Interpolation limits
    double logE_ref, Xmax_ref;

    // ms0,ms1,ms2,ms_logE
    double TimeModelPar_Stage1[4][nParDist][2][nPar_Stage1];
    // ms_slope  ms_psi_up , ms_psi_down
    double TimeModelPar_Stage2[4][nParDist][nParTheta][2][nPar_Stage2], eTimeModelPar_Stage2[4][nParDist][nParTheta][2][nPar_Stage2];
    double TimeModelPar_Stage2_Theta[4][nParDist][2][nPar_Stage2][nPar_Stage2_Theta];

    double BinomialNorm[200][99];
    double Binomial_q[200][99][3];


  public:

    UnivParamTime(int DetectorType);
    ~UnivParamTime() {};


    double GetD_TO(double X, double theta, double logE, int icomp);

    double GetMS(double DL, double logE, double theta, double r, double psi , int icomp, int iMS, double* par_stage1, double* par_stage2);
    double GetMS_ir(double DL, double logE, double theta, int   ir, double psi , int icomp, int iMS);
    double GetMS_ir_ith(double DL, double logE, int      ith, int   ir, double psi , int icomp, int iMS);

    double interpol(double v1, double v2, double x1, double x2, double x);
    double GetMS(double DL, double logE, double theta, double r, double psi , int icomp, int iMS);

    bool GetShapeParameters(double* DL, double r, double logE, double psi, double theta,  double* m,  double* s);

    double GetFraction(double* t0, double* m, double* s, double* fcomp, double ti, bool isCDF);
    double GetFraction(double* DL, double r ,  double logE, double psi, double theta, double* fcomp, double* t0, double ti, bool isCDF);

    double GetCDF(double* t0, double* m, double* s, double* fcomp, double ti);
    double GetCDF(double* DL, double r ,  double logE, double psi, double theta, double* fcomp,  double* t0, double ti);

    double GetPDF(double* t0, double* m, double* s, double* fcomp, double ti);
    double GetPDF(double* DL, double r ,  double logE, double psi, double theta, double* fcomp,  double* t0, double ti);

    double GetTime(double* DL, double r ,  double logE, double psi, double theta, double* fcomp,  double* t0, double fi);
    double GetTime(double* t0, double* m, double* s, double* fcomp, double fi);

    double GetPoissonFactor(double fem, double f);
    double tQuantilePDF(double* t0, double* m, double* s, double* fcomp, double vemTot,  double ti, double f, bool UseApprox, double pf, double& mean, double& rms);
    double tQuantilePDF(double* t0, double* m, double* s, double* fcomp, double vemTot,  double ti, double f, bool UseApprox, double& mean, double& rms);
    double 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);

    double GetBinomialNorm(double N, double f , double* fq);
    double QuantilePDF(double N,  double fi, double f);
    double QuantilePDF(double N,  double fi, double f, double* fq);

    double tFirstPDF(double t0_mu, double m_mu, double s_mu, double npart, double t, bool UseApprox, double& mean, double& rms);
    double tFirstPDF(double DL_mu, double r, double logE, double psi, double theta, double t0_mu, double npart, double t);
    double 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);

    double tQuantileCorrection(double NanoSecPerVEM, int iq);
    double tStartCorrection(double r, double logE, double theta, bool Is8nsFADC);

    double tQuantileCorrection_AoP(double RiseTime, int iq, double AoP);

    void   SetOffsetM_Mu(double Offset)
    {
      fOffsetM_Mu = Offset;
    }
    double GetOffsetM_r(double r, int icomp);
  };
}

#endif // __UnivParamTime_h_