#include <iostream>
#include <iomanip>
using namespace std;
// for standalone//
/*
  //lut
*/
int fw_binned(double *aob, double *aobe, double &R0, double &sigma, double &phi0, double &birks, double initR0){
  // int my_precision=16; // NODIF IN 6 DIG for integ1 //def
  // int fw_nintegbin=4; // def But dont use calc by hand
  int my_precision=8; // NODIF IN 6 DIG for integ1 //def
  int fw_nintegbin=1; // def But dont use calc by hand
  double fw_precision = 1e-4; // a parameter for numerical integration, reduce to speed up calculation

  int fw_nbin=32; //just number of bin of data 
  //HHHH LUT gen HHHHH//
  // hist infor e.g. layer widht of detector
  double bins_lut[fw_nbin+1]=
    { 0,  2.90158,  5.7519,  8.78677,  11.6986,  14.5387,  17.5428,  20.4444,
      23.3152,  26.3193,  29.2209,  32.1122,  35.1061,  38.0179,  40.899,  43.8723,
      46.7739,  49.6755,  52.8437, 55.6837,  58.5033,  61.4561,  64.3885,  67.2798,
      70.2224,  73.1445,  75.9538,  78.6606,  81.3879,  84.4125,  87.109,  89.785,
      92.4508};
  // Calc from bin width divided by # of division to average in a bin for expected 
  double binw_lut[fw_nbin]=
    { 0.580316, 0.570064, 0.606974, 0.582366, 0.56802, 0.60082, 0.58032, 0.57416,
      0.60082, 0.58032, 0.57826, 0.59878, 0.58236, 0.57622, 0.59466, 0.58032,  0.58032,
      0.63364, 0.568, 0.56392, 0.59056, 0.58648, 0.57826, 0.58852, 0.58442, 0.56186,
      0.54136, 0.54546, 0.60492, 0.5393, 0.5352, 0.53316};


  //HHHH Minimize parameters HHHH//
  int debug_lv=0;
  double rx=0;
  int max_test=128;
  //  int max_test=64;  
  //HHHHH parameter min max setting
  // double cR0=50; //current val
  // double rR0=50; //range
  // double minR0=50; //current R0 of minimum -2logL
  double cR0=initR0; //current val
  double rR0=20; //range
  double minR0=initR0; //current R0 of minimum -2logL
  // double cR0=70; //current val
  // double rR0=20; //range
  // double minR0=70; //current R0 of minimum -2logL
  // double cR0=50; //current val
  // double rR0=50; //range
  // double minR0=50; //current R0 of minimum -2logL

  // double csigma=4; //current val
  // double rsigma=4; //range
  // double minsigma=4; //current sigma of minimum -2logL
  double csigma=2; //current val
  double rsigma=1; //range
  double minsigma=2; //current sigma of minimum -2logL
  // double csigma=5; //current val double rsigma=5; //range double
  // minsigma=5; //current sigma of minimum -2logL

  double cphi0=20; //current val *1e-7
  double rphi0=20; //range
  double minphi0=20; //current phi0 of minimum -2logL
  // double cphi0=10e-7; //current val
  // double rphi0=10e-7; //range
  // double minphi0=10e-7; //current phi0 of minimum -2logL

  double cbirks=0.07; //current val
  double rbirks=0.03; //range
  double minbirks=0.07; //current birks of minimum -2logL
  // double cbirks=0.07; //current val
  // double rbirks=0.03; //range
  // double minbirks=0.07; //current birks of minimum -2logL
  // double cbirks=0.25; //current val
  // double rbirks=0.25; //range
  // double minbirks=0.25; //current birks of minimum -2logL

  int tlv=0, stagelv=0, paran=0;
  double minLogL=9999999;
  int npar=4;
  for(int t=0; t<max_test; t++){
    //HHHH Stepping HHHH//
    tlv++;
    paran=(tlv/2)%npar;
    stagelv=(tlv-2)/(npar*2);
    if(tlv==1){      //init
    }else if(paran==0){    
      if(tlv%2==0)cphi0+=rphi0/(pow(2,stagelv)*2);
      else if(tlv%2==1)cphi0-=rphi0/(pow(2,stagelv));
    }else if(paran==1){    
      if(tlv%2==0)cbirks+=rbirks/(pow(2,stagelv)*2);
      else if(tlv%2==1)cbirks-=rbirks/(pow(2,stagelv));
    }else if(paran==2){    
      if(tlv%2==0)cR0+=rR0/(pow(2,stagelv)*2);
      else if(tlv%2==1)cR0-=rR0/(pow(2,stagelv));
    }else if(paran==3){    
      if(tlv%2==0)csigma+=rsigma/(pow(2,stagelv)*2);
      else if(tlv%2==1)csigma-=rsigma/(pow(2,stagelv));
    }
    //HHHH LogL calc HHHH//
    double logL=0;
    for(int ibin=1; ibin<fw_nbin+1; ibin++){
      double logP = 0;
      //HH observed HH//
      double ob=aob[ibin-1];
      double obe=aobe[ibin-1];
      //HH expected HH//
      double sum_integ1=0, sum_integ2=0;
      double bin_step = binw_lut[ibin];    //double bin_width=bins[ibin]-bins[ibin-1]; //double bin_step = bin_width/(nintegbin+1);
      for(int integbin=0; integbin<nintegbin; integbin++){
	rx=bins[ibin-1]+bin_step*(integbin+1);  // tierd @ here.... bellow is not optimized to FW.
	double zDef = 5.*csigma; //This constant determines in which area the integrator of fInt1 and fInt2 gonna be defined and integrated, see definition of a and b
	double a = min(rx-zDef,(1.-fw_precision-1e-3)*cR0); // integration start. Subtract constant to make sure a is smaller than b
	double b = min(rx+zDef,(1.-fw_precision)*cR0); // integration end
	// HH Bo integral  part HH// ---------------------------------------------------------------------------------------------------------
	double integral1=0;//integ 1 par
	double width1=(b-a)/my_precision; 
	double xx1=a+width1/2.;
	double integral2=0;//integ 2 par
	double pF2 = gamma_H2O*beta/(1+beta*cR0);
	double LofZ = 0;  
	double myLofZ=0;
	double width2=((1-fw_precision)*cR0-rx)/my_precision;
	double xx2=rx+width2/2.;
	for(int i=0; i<my_precision; i++){
	  // integ 1
	  double pF1 = (1.+beta*(cR0-xx1))/(1+beta*cR0);
	  double foldingTerm = 1./(2.507*csigma)*exp(-(rx-xx1)*(rx-xx1)/(2*csigma*csigma));
	  double tdLdz=1./(0.2117*pow(cR0-xx1,0.425)+cbirks);
	  double val1 = pF1 * tdLdz * foldingTerm;
	  integral1+=val1*width1;
	  xx1+=width1;
	  // integ 2
	  if(rx<(1-fw_precision)*cR0){
	    double val2  = 1./(0.2117*pow(cR0-xx2,0.425)+cbirks);
	    myLofZ+=val2*width2;
	    xx2+=width2;
	  }
	}
	integral2 = pF2 * myLofZ;
	sum_integ1+=integral1;
	sum_integ2+=integral2;
      }
      double integral = (sum_integ1 + sum_integ2)/nintegbin; // Put all parts together and add scaling factor    
      // Combine prefactor and integral of all parts
      double preFactor = cphi0; //S*cphi0/rho //Put all parts together and add scaling factor
      double ex = preFactor*integral;
      //HH -2 Log gaussian likelihood calculatio to estimate mini L position for minimizer
      if(ex>0){
	logP=(ob-ex)*(ob-ex)/(obe*obe)+2*log(obe*2.507);
      }
      logL+=logP;
    } // end loop for bins for minL
    if(debug_lv>3){
      cout<<setprecision(4)
	  <<"stage="<<setw(3)<<stagelv<<", paran="<<paran<<", tlv="<<setw(5)<<tlv
	  <<", R0="<<setw(8)<<cR0<<", minR0="<<setw(8)<<minR0
	  <<", sigma="<<setw(8)<<csigma<<", minsigma="<<setw(8)<<minsigma
	  <<", phi0="<<setw(8)<<cphi0<<", minphi0="<<setw(8)<<minphi0
	  <<", birks="<<setw(8)<<cbirks<<", minbirks="<<setw(8)<<minbirks
	  <<",  logL="<<setw(8)<<logL<<", MinLogL="<<setw(8)<<minLogL<<endl;
    }
    //HHHH minimization HHHH//
    if(logL<minLogL){ //mini
      minLogL=logL;
      minR0=cR0;
      minsigma=csigma;
      minphi0=cphi0;
      minbirks=cbirks;
    }
    if(paran==0 and tlv%2==1){
      cphi0=minphi0;
    }else if(paran==1 and tlv%2==1){
      cbirks=minbirks;
    }else if(paran==2 and tlv%2==1){
      cR0=minR0;
    }else if(paran==3 and tlv%2==1){
      csigma=minsigma;
    }
  }
  R0=minR0; sigma=minsigma; phi0=minphi0; birks=minbirks;

  return 0;
}






double GetLL(double *aob, double *aobe, double cR0, double csigma, double cphi0, double cbirks){

  int my_precision=16; // NODIF IN 6 DIG for integ1 //def
  int fw_nintegbin=4; // def But dont use calc by hand
  double fw_precision = 1e-4; // a parameter for numerical integration, reduce to speed up calculation

  int fw_nbin=32; //just number of bin of data 
  //HHHH LUT gen HHHHH//
  // hist infor e.g. layer widht of detector
  double bins_lut[fw_nbin+1]=
    { 0,  2.90158,  5.7519,  8.78677,  11.6986,  14.5387,  17.5428,  20.4444,
      23.3152,  26.3193,  29.2209,  32.1122,  35.1061,  38.0179,  40.899,  43.8723,
      46.7739,  49.6755,  52.8437, 55.6837,  58.5033,  61.4561,  64.3885,  67.2798,
      70.2224,  73.1445,  75.9538,  78.6606,  81.3879,  84.4125,  87.109,  89.785,
      92.4508};
  // Calc from bin width divided by # of division to average in a bin for expected 
  double binw_lut[fw_nbin]=
    { 0.580316, 0.570064, 0.606974, 0.582366, 0.56802, 0.60082, 0.58032, 0.57416,
      0.60082, 0.58032, 0.57826, 0.59878, 0.58236, 0.57622, 0.59466, 0.58032,  0.58032,
      0.63364, 0.568, 0.56392, 0.59056, 0.58648, 0.57826, 0.58852, 0.58442, 0.56186,
      0.54136, 0.54546, 0.60492, 0.5393, 0.5352, 0.53316};

  //HHHH Minimize parameters HHHH//
  int debug_lv=0;
  double rx=0;

  //  int max_test=64;  
    //HHHH LogL calc HHHH//
    double logL=0;
    for(int ibin=1; ibin<fw_nbin+1; ibin++){
      double logP = 0;
      //HH observed HH//
      double ob=aob[ibin-1];
      double obe=aobe[ibin-1];
      //HH expected HH//
      double sum_integ1=0, sum_integ2=0;
      double bin_step = binw_lut[ibin];    //double bin_width=bins[ibin]-bins[ibin-1]; //double bin_step = bin_width/(nintegbin+1);
      for(int integbin=0; integbin<nintegbin; integbin++){
	rx=bins[ibin-1]+bin_step*(integbin+1);  // tierd @ here.... bellow is not optimized to FW.
	double zDef = 5.*csigma; //This constant determines in which area the integrator of fInt1 and fInt2 gonna be defined and integrated, see definition of a and b
	double a = min(rx-zDef,(1.-fw_precision-1e-3)*cR0); // integration start. Subtract constant to make sure a is smaller than b
	double b = min(rx+zDef,(1.-fw_precision)*cR0); // integration end
	// HH Bo integral  part HH// ---------------------------------------------------------------------------------------------------------
	double integral1=0;//integ 1 par
	double width1=(b-a)/my_precision; 
	double xx1=a+width1/2.;
	double integral2=0;//integ 2 par
	double pF2 = gamma_H2O*beta/(1+beta*cR0);
	double LofZ = 0;  
	double myLofZ=0;
	double width2=((1-fw_precision)*cR0-rx)/my_precision;
	double xx2=rx+width2/2.;
	for(int i=0; i<my_precision; i++){
	  // integ 1
	  double pF1 = (1.+beta*(cR0-xx1))/(1+beta*cR0);
	  double foldingTerm = 1./(2.507*csigma)*exp(-(rx-xx1)*(rx-xx1)/(2*csigma*csigma));
	  double tdLdz=1./(0.2117*pow(cR0-xx1,0.425)+cbirks);
	  double val1 = pF1 * tdLdz * foldingTerm;
	  integral1+=val1*width1;
	  xx1+=width1;
	  // integ 2
	  if(rx<(1-fw_precision)*cR0){
	    double val2  = 1./(0.2117*pow(cR0-xx2,0.425)+cbirks);
	    myLofZ+=val2*width2;
	    xx2+=width2;
	  }
	}
	integral2 = pF2 * myLofZ;
	sum_integ1+=integral1;
	sum_integ2+=integral2;
      }
      double integral = (sum_integ1 + sum_integ2)/nintegbin; // Put all parts together and add scaling factor    
      // Combine prefactor and integral of all parts
      double preFactor = cphi0; //S*cphi0/rho //Put all parts together and add scaling factor
      double ex = preFactor*integral;
      //HH -2 Log gaussian likelihood calculatio to estimate mini L position for minimizer
      if(ex>0){
	logP=(ob-ex)*(ob-ex)/(obe*obe)+2*log(obe*2.507);
      }
      logL+=logP;
    } // end loop for bins for minL

    return logL;
}