#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 {

  //  Parameterization
  static const int nPar=8;
  static const int nParAng=3;
  static const double TimeModelParAng[nParAng]={0.,M_PI/2.,M_PI};
  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};
  
  // 
  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=18.,logE_High=21.;
  static const double r_Low=10.e2;
  static const double r_High=2200.e2;

  // Systematics   (for the WCD)
  static const double Systematic_ShiftM_Models[4]={  0.00  , -0.015 ,  0.03  , 0.015  }; 
  static const double Systematic_OffsetT1Mu_Models[4]={ 0.0,  3.00,  1.50,  6.50 };    // Intrinsic
  static const double Systematic_C_XmaxMu_Models[4]={ 0.85 , 1.00, 0.85 , 1.00 };  
  
  //-------------------------
  //-------------------------
  //-------------------------

  class UnivParamTime {
  
  private:
    int fDetectorType; // [0] WCD 
    int fMonth; // Month in the year

    double TimeModelPar_mod1[4][2][2],TimeModelPar_mod2[4][2][2],TimeModelPar_logE[4][2][2]; // [icomp][iMS][ipar_r]
    double TimeModelParDX[4][nParDist][2][2];                                                // [icomp][ir][iMS][ms0/ms1]

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

  
  public:

    UnivParamTime(int DetectorType);
    ~UnivParamTime(){};
    
    void SetMonth(int month_i){ fMonth=month_i;}

    double GetT1(     double Xfirst, double DX, int icomp, int im, int ip);  
    double GetT1_Ref( double Xmax,   double logE,  double DX, int icomp, double OffsetMu );

    int GetPos(const double *pars,double dist);
    double interpol(double v1,double v2,double x1,double x2, double x);

    bool   InterpolateDXParameters(double DX, double r, int iMS, int icomp, double &ms0, double &ms1);
    double GetMS(double DX, double r, double logE, double psi, double theta, int iMS, int icomp);
    bool GetShapeParameters(double *DX, 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 *DX, 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 *DX, 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 *DX, double r ,  double logE, double psi, double theta, double *fcomp,  double *t0, double ti);

    double GetTime(double *DX, 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 *DX, 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 DX_mu, double r, double logE, double psi, double theta, double t0_mu, double npart, double t);
    double tFirstPDF( double DX_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);

    double GetSystematicShift_r(double r, int icomp, double SystematicShift );
    double GetAtmosphericCorrection(double r, double imonth);
    
    bool   GetPar_FIP_Xmax(double logE,int im0,int ip0, double &alpha, double &beta );
    double GetMoyalProb(double X, double alpha, double beta);
  };
}