#include "fQChrgPDF.h"

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <bitset>
#include <stdio.h>

using namespace std;

//ClassImp(fQChrgPDF);

fQChrgPDF* fQChrgPDF::staticthis = NULL;

bool fQChrgPDF::flgUseGlParArr=false;

int fQChrgPDF::iPMTType_gl;
Double_t fQChrgPDF::qval_gl;// Current charge q
Double_t fQChrgPDF::param_gl[fQChrgPDF::nParamMax];

fQChrgPDF::fQChrgPDF(){
  
  flgLoaded=false;
  
  flgStrict = false;
  
  flgUseGlParArr=false;
  
  cPDFparf = NULL;
  
}

fQChrgPDF::~fQChrgPDF(){
  
  if (cPDFparf) delete cPDFparf;
  
  cout << "Destructing fQChrgPDF!!" << endl;
}

void fQChrgPDF::LoadParams(char *ifName) {
  
  cout << "fQChrgPDF: Loading parameters from " << ifName << endl;
  
  if (cPDFparf) delete cPDFparf;
  cPDFparf = new TFile(ifName);//read charge PDF parameter graphs
  if (cPDFparf->IsZombie()) {
    cout << "Charge PDF parameter file does not exist!" << endl;
//    DeleteInstances();
    exit(-1);
  }
  
  for (int iPMTType=0; iPMTType<nPMTType; iPMTType++) {
    TH1D *h_PunhitPar = (TH1D*)(cPDFparf->Get(Form("hPunhitPar_type%d",iPMTType)));
    if (h_PunhitPar) {
      cout << Form("P_unhit coeff. for Type %d PMT: c_n=[",iPMTType);
      for (int i=0; i<nPunhitOrd; i++) {
        Punhit_cn[iPMTType][i] = h_PunhitPar->GetBinContent(i+1);
        cout << Form("%f, ",Punhit_cn[iPMTType][i]);
      }
      cout << "]" << endl;
    }
    else continue;
    
    hCPDFrange[iPMTType] = new TH1D(*(TH1D*)(cPDFparf->Get(Form("hCPDFrange_type%d",iPMTType))));
    
    int nRang = hCPDFrange[iPMTType]->GetXaxis()->GetNbins();
    
    for (int iRang=0; iRang<=nRang; iRang++) {
      int nparam = (int)(hCPDFrange[iPMTType]->GetBinContent(iRang+1)+1e-8);
      
      for (int i=0; i<nparam; i++) {
        gParam[iPMTType][iRang][i] = new TGraph(*(TGraph*)(cPDFparf->Get(Form("gParam_type%d_Rang%d_%d",iPMTType,iRang,i))));
      }
      
      for (int k=0; k<2; k++) {
        gmuthr[iPMTType][iRang][k] = new TGraph(*(TGraph*)(cPDFparf->Get(Form("gmuthr_type%d_Rang%d_%d",iPMTType,iRang,k))));
      }
    }
    
  }
  
  flgLoaded=true;
  
}


// Wrapper function for TF1 ###########################
Double_t fQChrgPDF::flogcPDFptr(Double_t *x, Double_t *par){//log of charge PDF
  
//  cout << "flogcPDFptr is called!" << endl;
  
  if (!staticthis) {
    cout << "fQChrgPDF: Not ready!" << endl;
    return 0.;
  }
  
  for (int i=0; i<nParamMax; i++) {
    param_gl[i] = par[i];
  }
  
  Double_t muval=*x;// 1d function of mu
  
  flgUseGlParArr=true;// use param_gl instead of gParam
  Double_t tmp = staticthis->flogcPDF(muval,qval_gl,iPMTType_gl);
  flgUseGlParArr=false;
  
  return tmp;
  
}


double fQChrgPDF::flogcPDF(double mu, double q, int iPMTType){//log of charge PDF
  
  int nparam;
  double par[nParamMax];
  double muthr[2];
  double coeff[6];
  
  SetCPDFParams(q,iPMTType,nparam,par,muthr,coeff);
  
  return flogchrgPDF(mu,nparam,par,muthr,coeff);
}


double fQChrgPDF::flogchrgPDF(double mu, int nparam, double *par, double *thr, double *coeff){//log of charge PDF
  
  double tmp=0.;
  
  double lth=thr[0];
  double hth=thr[1];
  
  if (mu>hth) {//high mu
    tmp=coeff[3]+(mu-hth)*(coeff[4]+(mu-hth)*coeff[5]);
  }
  else if (mu<lth) {//low mu
    tmp=coeff[0]+(mu-lth)*(coeff[1]+(mu-lth)*coeff[2]);
  }
  else {
    double mupow=1.;
    for (int i=0; i<nparam; i++) {
      tmp+=par[i]*mupow;
      mupow*=mu;
    }
  }
  
  return tmp;
}

int fQChrgPDF::GetCPDFRange(double tmpchrg, int iPMTType, int &iRang, int &nparam) {// c++ indexing!
  
  int nRang = hCPDFrange[iPMTType]->GetXaxis()->GetNbins();// # of q range(excl. ovfl. bin)
  
  iRang = hCPDFrange[iPMTType]->GetXaxis()->FindBin(tmpchrg);
  iRang--;
  
  if (iRang<0) iRang=0;
  if (iRang>=nRang) iRang=nRang;// overflow!
  
  nparam = (int)(hCPDFrange[iPMTType]->GetBinContent(iRang+1)+1e-8);
  
  return nRang;
  
}

void fQChrgPDF::SetCPDFParams(double tmpchrg, int iPMTType, int &nparam, double *par, double *muthr, double *coeff){
  
  //  int icab=ihitPMT[ih];
  //  int iPMTType = fqshared->PMTtype[icab];
  //  double tmpchrg=chrg[icab];
  //  int &nparam = cpdfnpar[icab];
  //  double *par = cpdfpar[icab];
  //  double *muthr = cpdfthr[icab];
  //  double *coeff = cpdfcoeff[icab];
  
  if (!flgLoaded && !flgUseGlParArr) {
    cout << "fQChrgPDF: Parameters not loaded!" << endl;
    exit(-1);
  }
  
  if (!(tmpchrg>0.05)) tmpchrg=0.05;
  
  int iRang;
  int nRang = GetCPDFRange(tmpchrg,iPMTType,iRang,nparam);
  
  for (int i=0; i<nparam; i++) {
    if (flgUseGlParArr) {// for fitting the parameters via flogcPDFptr
      par[i] = param_gl[i];
    }
    else {
      par[i] = gParam[iPMTType][iRang][i]->Eval(tmpchrg,0,"");
    }
  }
  
  for (int i=0; i<6; i++) coeff[i]=0.;
  
  for (int k=0; k<2; k++) {//0:lower tail, 1:upper tail
    muthr[k] = gmuthr[iPMTType][iRang][k]->Eval(tmpchrg);
    
    for (int i=0; i<nparam; i++) {
      coeff[0+3*k]+=par[i]*pow(muthr[k],i);//0th order
    }
    if (k==0 || iRang<nRang) {// for overflow bin, flat upper tail(saturation)
      for (int i=1; i<nparam; i++) {
        coeff[1+3*k]+=i*par[i]*pow(muthr[k],i-1);//1st order
      }
      for (int i=2; i<nparam; i++) {
        coeff[2+3*k]+=i*(i-1)*par[i]*pow(muthr[k],i-2)/2;//2nd order
      }
    }
    if (coeff[2+3*k]>0.) {// make sure the tails don't concave up!
      if (flgStrict) {
        cout << Form("fQChrgPDF::SetCPDFParams: Tail concaves up! q=%f, a_2=%f",tmpchrg,coeff[2+3*k]) << endl;
//        exit(-1);
      }
      coeff[2+3*k]=0.;
    }
  }
}