#ifndef _utl_MuonArrivalTime_h_
#define _utl_MuonArrivalTime_h_

#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/Integrator.h>
#include <utl/Math.h>
#include <utl/NormalDistribution.h>
#include <limits>

// based on original code and analysis of Lorenzo Cazon
// rewritten and extended by Hans Dembinski

// not available in some compiler collections
namespace std {
  inline double pow(const double x, const unsigned int i) { return std::pow(x, int(i)); }
}

namespace utl {

  /// all time delays are relative to arrival time of shower front plane
  class MuonArrivalTime {
  public:

    MuonArrivalTime()
      : fTheta(0.0),
        fR(0.0),
        fDelta(0.0),
        fLogLambda(0.0),
        fLogMean(0.0),
        fLogSigma(0.0),
        fNorm(1.0),
        fMean(0.0),
        fSigma(0.0),
        fTimePDFArg(*this),
        fTimePDF(fTimePDFArg),
        fTimeCDFArg(*this),
        fTimeCDF(fTimeCDFArg),
        fSuperArg(*this),
        fSuper(fSuperArg),
        fApproxMomentArg(*this),
        fApproxMoment(fApproxMomentArg)
    {
      // do not change these
      fTimePDF.SetAccuracy(1e-2);
      fTimeCDF.SetAccuracy(1e-5);
      fSuper.SetAccuracy(1e-5);
      fApproxMoment.SetAccuracy(1e-2);
    }

    /// setup of model for given zenith angle
    void
    SetTheta(const double theta)
    {
      fTheta = theta;
      const double c = std::cos(theta);

      const double pars1[] = { 5.682, -5.524, 9.180, -8.277, 2.797 };
      fLogMean = EvalPoly(c, 5, pars1);

      const double pars2[] = { 0.03266, 0.456, -1.084, 1.199, -0.4366 };
      fLogSigma = EvalPoly(c, 5, pars2);

      const double pars3[] = { -0.01111, 0.2581, -0.7385, 3.362, -2.261 };
      fLogLambda = EvalPoly(c, 5, pars3);
      if (fLogLambda < 0.03) fLogLambda = 0.03;
    }

    /// setup model for given zenith angle and mean production distance
    void
    SetThetaAndDistance(const double theta, const double distance)
    {
      SetTheta(theta);
      fLogMean = std::log10(distance);
    }

    /// set coordinates of detector (shower front plane coordinate system)
    /// r = radial distance to shower axis, delta = distance to shower front plane
    void
    SetCoordinates(const double r, const double delta)
    {
      fR = r;
      fDelta = delta;

      // compute normalisation and 1rst and 2nd moment
      double result[4] = { 0, 0, 0, 0 };
      fNorm = 1.0;
      fStart = -20;
      for (int i = fStart; i < 20; ++i)
      {
        double delta[4];
        fSuper.GetIntegral(delta, i, i+1);
        result[0] += delta[0];
        result[1] += delta[1];
        result[2] += delta[2];
        result[3] += delta[3];
        if (result[0] == 0) fStart = i+1;
        if (delta[0] < 1e-5*result[0]) break;
      }

      fNorm = result[0];
      fMean = result[2]/result[1];
      fSigma = std::sqrt(result[3]/result[1]-fMean*fMean);
    }

    /// PDF of arrival time of muons
    double
    TimePDF(const double t)
      const
    {
      if (t <= 0) return 0;
      fTimePDFArg.fT = t;
      return fTimePDF.GetIntegral(0, t);
    }

    /// CDF of arrival time of muons
    double
    TimeCDF(const double t)
      const
    {
      if (t <= 0) return 0;
      const double lnt = std::log(t);
      double result = 0;
      for (int i = fStart; i < lnt; ++i)
      {
        const double delta = fTimeCDF.GetIntegral(i, i+1 < lnt ? i+1 : lnt);
        result += delta;
        if (delta < 1e-5*result) break;
      }
      return result;
    }

    /// extreme value distribution of arrival time of first particle
    double
    FirstTimePDF(const double t, const unsigned int nmuons)
      const
    {
      if (nmuons <= 1) return TimePDF(t);
      return nmuons*std::pow(1.0 - TimeCDF(t), nmuons-1)*TimePDF(t);
    }

    /// approximate PDF of arrival time of muons using normal approximation in log(t)
    double
    ApproximateTimePDF(const double t)
      const
    {
      if (t <= 0) return 0;
      // beware: normal distribution in log(t), not log-normal distribution!
      return utl::NormalPDF(std::log(t),fMean,fSigma)/std::exp(fMean+0.5*fSigma*fSigma);
    }

    /// CDF of arrival time of muons using normal approximation in log(t)
    double
    ApproximateTimeCDF(const double t)
      const
    {
      if (t <= 0) return 0;
      const double z = (std::log(t)-fMean)/fSigma - fSigma;
      return 0.5*(1.0+erf(z/kSqrt2));
    }

    /// extreme value distribution of arrival time of first particle using normal approximation in log(t)
    double
    ApproximateFirstTimePDF(const double t, const unsigned int nmuons)
      const
    {
      if (nmuons <= 1) return ApproximateTimePDF(t);
      return nmuons*std::pow(1.0 - ApproximateTimeCDF(t), nmuons-1)*ApproximateTimePDF(t);
    }

    void
    GetApproximateMeanAndStDev(double& mean, double& stdev, const unsigned int nmuon)
      const
    {
      fApproxMomentArg.fNMuon = nmuon;
      double result[2] = { 0, 0 };
      for (int i = fStart; i < 20; ++i)
      {
        double delta[2];
        fApproxMoment.GetIntegral(delta, i, i+1);
        result[0] += delta[0];
        result[1] += delta[1];
        if (delta[1] < 1e-5*result[1]) break;
      }

      mean = result[0];
      stdev = std::sqrt(result[1] - mean*mean);
    }

  private:

    /// PDF of geometric time delay due to distribution of muon production distance
    double
    PDF_dNdt_G(const double t)
      const
    {
      if (t <= 0) return 0.0;

      const double c = utl::kSpeedOfLight;
      const double r = fR;

      // t is time delay after shower front (for ct << z?)
      const double z0 = 0.5*(r*r/(c*t)-c*t);
      const double z = z0 + fDelta;
      if (z <= 0) return 0.0;

      const double dzdt = 0.5*(r*r/(c*t*t)+c);

      const double cosaDa = std::sqrt(1.0-r*r/(z0*z0 + r*r));

      const double lgZ = std::log10(z);
      const double dNdz = 1.0/(z*utl::kLn10)*PDF_dNdlgZ(lgZ);

      return cosaDa*dNdz*dzdt;
    }

    /// PDF of kinetic time delay due to distribution of muon energy
    double
    PDF_dNdt_E(const double t)
      const
    {
      if (t <= 0) return 0.0;

      const double zMean = std::pow(10, fLogMean);
      const double z0 = zMean - fDelta;
      const double x = std::sqrt(z0*z0 + fR*fR);

      const double m2 = 0.011;
      const double pk = 0.0002;
      const double c = utl::kSpeedOfLight;

      const double energy = 0.5*pk*x*(1.0+std::sqrt(1.0+2.0*m2/(pk*pk*x*c*t)));
      const double dEdt = 0.5*pk*x*1.0/std::sqrt(1.0+2.0*m2/(pk*pk*x*c*t))*m2/(pk*pk*x*c*t*t);

      return PDF_dNdE(energy, x, fR)*dEdt;
    }

    /// PDF model of muon production distance
    double
    PDF_dNdlgZ(const double logz)
      const
    {
      const double f1 = (logz-fLogMean)/fLogLambda;
      const double f2 = 0.5*(fLogSigma*fLogSigma)/(fLogLambda*fLogLambda);
      const double f3 = (logz-fLogMean)/(utl::kSqrt2*fLogSigma);
      const double f4 = fLogSigma/(utl::kSqrt2*fLogLambda);

      return std::exp(f1+f2)*erfc(f3+f4)/(2.0*fLogLambda);
    }

    /// PDF model of muon energy distribution
    double
    PDF_dNdE(const double energy, const double x, const double r)
      const
    {
      const double pk=0.0002;
      const double kappa=0.8;
      const double lamb=1;
      const double ptc=0.17;
      const double gam=2.6;
      return (1.0-r*r/(x*x))
            *std::pow(r/x,lamb)
            *std::pow(energy,-gam-kappa+lamb+1)
            *std::pow(energy-pk*x,kappa)
            *std::exp(-energy*r/(x*ptc));
    }

    // Utiltities
    double
    EvalPoly(const double x, const int n, const double* pars)
      const
    {
      double result = pars[n-1];
      for (int i = 1; i < n; ++i)
      {
        result *= x;
        result += pars[n-1-i];
      }
      return result;
    }

    struct TimePDFArg {
      const MuonArrivalTime& fParent;
      double fT;

      TimePDFArg(const MuonArrivalTime& parent)
        : fParent(parent)
      {}

      double
      operator()(const double t)
        const
      {
        return fParent.PDF_dNdt_E(t)*fParent.PDF_dNdt_G(fT-t)/fParent.fNorm;
      }
    };

    struct TimeCDFArg {
      const MuonArrivalTime& fParent;

      TimeCDFArg(const MuonArrivalTime& parent)
        : fParent(parent)
      {}

      double
      operator()(const double lnt)
        const
      {
        const double t = std::exp(lnt);
        fParent.fTimePDFArg.fT = t;
        return t*fParent.fTimePDF.GetIntegral(0, t);
      }
    };

    struct SuperArg {
      const MuonArrivalTime& fParent;

      SuperArg(const MuonArrivalTime& parent)
        : fParent(parent)
      {}

      void
      operator()(double* result, const double lnt)
        const
      {
        const double t = std::exp(lnt);
        fParent.fTimePDFArg.fT = t;
        const double f = fParent.fTimePDF.GetIntegral(0, t);

        result[0] = f*t;       // norm of f(t) d t
        result[1] = f;         // norm of f(t) d ln(t)
        result[2] = f*lnt;     // mean of f(t) d ln(t)
        result[3] = f*lnt*lnt; // var  of f(t) d ln(t)
      }
    };

    struct ApproxMomentArg {
      const MuonArrivalTime& fParent;
      unsigned int fNMuon;

      ApproxMomentArg(const MuonArrivalTime& parent)
        : fParent(parent)
      {}

      void
      operator()(double* result, const double lnt)
        const
      {
        const double t = std::exp(lnt);
        const double f = fParent.ApproximateFirstTimePDF(t, fNMuon)*t;

        result[0] = f*t;
        result[1] = f*t*t;
      }
    };

  private:
    double fTheta;
    double fR;
    double fDelta;
    double fLogLambda;
    double fLogMean;
    double fLogSigma;
    double fNorm;
    double fStart;
    double fMean;
    double fSigma;
    mutable TimePDFArg fTimePDFArg;
    utl::Integrator<TimePDFArg> fTimePDF;
    mutable TimeCDFArg fTimeCDFArg;
    utl::Integrator<TimeCDFArg> fTimeCDF;
    mutable SuperArg fSuperArg;
    utl::VectorIntegrator<SuperArg, 4> fSuper;
    mutable ApproxMomentArg fApproxMomentArg;
    utl::VectorIntegrator<ApproxMomentArg, 2> fApproxMoment;
  };

}

#endif