#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Random/RandGauss.h>
#include <CLHEP/Random/RandGeneral.h>
#include <CLHEP/Random/RandFlat.h>
using namespace CLHEP;

#include <G4UImessenger.hh>

#include <RAT/PMTCharge.hh>
#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/DS/DataNotFound.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/ChanHWStatus.hh>

#include <TH1.h>
#include <TMath.h> //Call Function: TMath::Gamma(a,b);

#include <cmath>
#include <algorithm>
#include <vector>

namespace RAT {

void PMTCharge::BeginOfRun()
{
  const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();

  fQlo = 0.001;
  fQhi = 10;

  fPmtCalib = DB::Get()->GetLink("PMTCALIB");
  fMeanHHP = fPmtCalib->GetD("mean_hhp");
  float min_hhp = fPmtCalib->GetD("min_hhp");
  float max_hhp = fPmtCalib->GetD("max_hhp");

  // Default is to use measured PMT-specific HHP values
  fUseQHSHHP = true;
  DBLinkPtr ldaqsettings = DB::Get()->GetLink("DAQ");
  try {
    fUseQHSHHP = (bool) ldaqsettings->GetI("use_qhs_hhp");
  }
  catch (DBNotFoundError& e) {}

  if (fUseQHSHHP) {
    fQhshhp = fPmtCalib->GetFArrayFromD("QHS_hhp");
    int length = fQhshhp.size();

    // Ensure that all channels have HHPs
    DBLinkPtr lpmtinfo = DB::Get()->GetLink("PMTINFO");
    std::vector<int> pmttype = lpmtinfo->GetIArray("pmt_type");
    Log::Assert(length == (int)pmttype.size(),
                "PMTCharge: Number of HHP values in PMTCALIB doesn't match the number of PMTs in PMTINFO");

    for(int ip=0;ip<length;++ip){
      fQhshhp[ip] = (float)(int)(fQhshhp[ip]+0.5);
      if(fQhshhp[ip]==0.0)fQhshhp[ip]=DS::INVALID;
      if(fQhshhp[ip]<min_hhp&&fQhshhp[ip]>0.0)fQhshhp[ip]=min_hhp;
      if(fQhshhp[ip]>max_hhp)fQhshhp[ip]=max_hhp;
      bool hitpmt = true;
      if( channelHardwareStatus.IsEnabled() )
        hitpmt = channelHardwareStatus.IsDAQEnabled( ip );
      if(!hitpmt)fQhshhp[ip]=0.0;
      if(fQhshhp[ip]<0.0){
        if(hitpmt)fQhshhp[ip]=fMeanHHP;
      }
    }
  }

  // Variables below only needed for old charge distribution
  const int Qintervals = 200;
  double Qstep = (fQhi - fQlo) / Qintervals;

  double dist_array[Qintervals];
  for (int i = 0; i < Qintervals; i++)
      dist_array[i] = DoublePolyaQSpectrum(fQlo + Qstep * i);

  fQdist = new CLHEP::RandGeneral(dist_array, Qintervals);
  fM = 8.1497;
  fA1 = 0.27235E-01;
  fB1 = 0.7944E-03;
  fC1 = 0.81228E-01;
  fQP0 = .59;
  fQPa = 4.41;
  fGamM = TMath::Gamma(fM);
  fBeta = pow(fM,fM)/fGamM;
}

PMTCharge::~PMTCharge()
{
  delete fQdist;
}

double PMTCharge::Gamma(const double A, const double X)
{
  double d = TMath::Gamma(A);
  double g = TMath::Gamma(A,X);
  return d*(1.0-g);
}

double PMTCharge::PIntegrate(const double b, const int ipmt)
{
  const double eps = 1.0e-6;

  // In practice, only the first 8 or so elements are used.  This lets
  // it go all the way up to the point where the below code would break,
  // and only wastes 32k of memory.
  static double v[64][64];

  v[0][0] = 0.5*b*(DoublePolyaQSpectrum(0, GetHHP(ipmt))+
                   DoublePolyaQSpectrum(b, GetHHP(ipmt)));

  int n = 1;
  while(true){
    int nraise = (int) 1 << n;
    const double h = b/nraise;

    v[n][0] = 0.5*v[n-1][0];
    nraise = (int) 1 << (n - 1);
    for( int i = 1; i < nraise + 1; i++ )
      v[n][0] += h * DoublePolyaQSpectrum( (2*i-1)*h, GetHHP(ipmt));

    for (int m = 1; m < n+1; m++) {
      nraise = (int) 1 << ( 2 * ( m - 1 ) );
      v[n][m] = v[n][m-1] + (v[n][m-1]-v[n-1][m-1])/nraise;
    }

    if (fabs(v[n][n-1]-v[n][n]) < eps)
      return v[n][n];

    n++;
  }
}


float PMTCharge::Bisection(const double flat, const int ipmt, const double norm)
{
  float xl = 0.000000001;
  float xh =500;
  if (flat-norm*PMTCharge::PIntegrate(20.0, ipmt)<0){
     xh=20.;
  }
  else if (flat-norm*PMTCharge::PIntegrate(30.0, ipmt)<0){
     xl=20.;
     xh=30.;
  }
  else if (flat-norm*PMTCharge::PIntegrate(50.0, ipmt)<0){
     xl=30.;
     xh=50.;
  }
  else if (flat-norm*PMTCharge::PIntegrate(90.0, ipmt)<0){
     xl=50.;
     xh=90.;
  }
  else{
     xl=90.;
  }
  double currentq = .5*(xh+xl);
  double dxold = fabs(xh-xl);
  double dx = dxold;
  double f = -flat+norm*PMTCharge::PIntegrate(currentq, ipmt);
  double df = norm*DoublePolyaQSpectrum(currentq, GetHHP(ipmt));
  for (int i=0;i<100;i++) {
      if (((currentq - xh)*df-f)*((currentq-xl)*df-f) >= 0) {
         dxold = dx;
         dx=0.5*(xh-xl);
         currentq = xl+dx;
         if (xl == currentq) return currentq;
      }
      else if ( fabs(2.*f) > fabs(dxold*df)) {
         dxold = dx;
         dx=0.5*(xh-xl);
         currentq = xl+dx;
         if (xl == currentq) return currentq;
      }
      else{
         dxold = dx;
         dx = f/df;
         currentq = currentq-dx;
         if (dx==0) return currentq;
      }
      if (fabs(dx)<0.01) { // 1% accuracy is good enough (same assumption as for SNO)
         return currentq;
      }
      f = -flat+norm*PMTCharge::PIntegrate(currentq, ipmt);
      df = norm*DoublePolyaQSpectrum(currentq, GetHHP(ipmt));
      if (f < 0) {
         xl = currentq;
      }
      else{
         xh = currentq;
      }
  }
  //cout << "Exceeded maximum iterations in PMTCharge::Bisection" << endl;
  return 0.0;
}

// Calculate the charge induced by a photon striking a PMT, using the dipolya charge spectrum

float PMTCharge::CalcCharge(int /*npe*/, bool /*Polya*/, const int ipmt) // New charge in ADC counts including variable gain
{
  const double norm = 1.0/PMTCharge::PIntegrate(500, ipmt);
  const double rFlat = CLHEP::RandFlat::shoot(1);
  return PMTCharge::Bisection(rFlat, ipmt, norm);
}

// Old charge calculation

float PMTCharge::CalcCharge(const int npe) // Old charge in pe
{
  float pedwidth = 0.1;
  float qtot = 0.;

  if (npe>0){
     double* qran = new double[npe];
     double* qsmear = new double[npe];
     CLHEP::RandGauss::shootArray(npe, qsmear);
     fQdist->shootArray(npe, qran); // Random numbers in [0,1) range
     for (int ipe=0; ipe<npe; ipe++){
         qtot += /* charge distribution */ fQlo + qran[ipe] * (fQhi - fQlo)
                 /* pedestal smearing */   + pedwidth * qsmear[ipe];
     }
     delete[] qran;
     delete[] qsmear;
     if (qtot<0) qtot=0;
  }
  else {
     qtot += pedwidth * CLHEP::RandGauss::shoot();
  }
  return qtot;
}

float PMTCharge::DoublePolyaQSpectrum(const double Q)
{
  return PMTCharge::DoublePolyaQSpectrum(Q, 45.0);
}

float PMTCharge::DoublePolyaQSpectrum(const double Q, double hhp)
{
  if (hhp < 0 || hhp > 100)
    hhp=45;

  static const float m1 = 8.1497;
  static const float gmratio = m1/exp(lgamma(m1));
  static const float ng1 = 0.2735E-01;
  static const float ng2 = 0.7944E-03;
  static const float nexp = 0.81228E-01;
  static const float iQ0 = 1/20.116;

  const float G1 = 1.2525 + 0.71726*hhp;
  const float G2 = 118.21 + 0.71428*hhp;
  const float Qpe1 = Q/G1;
  const float Qpe2 = Q/G2;

  const float funpol1 = pow(m1*Qpe1, m1-1)*exp(-m1*Qpe1);
  const float funpol2 = pow(m1*Qpe2, m1-1)*exp(-m1*Qpe2);
  return gmratio*(ng1*funpol1 + ng2*funpol2) + nexp*exp(-Q*iQ0);
}

} // namespace RAT