#include <iostream>
#include <iomanip>
using namespace std;
#include "TH1D.h"
#include "TF1.h"



const int nbin=32;
int mytest(){
  
  double bins[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};

  TH1D *hd = new TH1D("hd","hd", nbin, bins);
		   
  hd->SetBinContent(0  , 0          ); 
  hd->SetBinContent(1  , 4.53902    ); 
  hd->SetBinContent(2  , 4.70824    ); 
  hd->SetBinContent(3  , 4.73104    ); 
  hd->SetBinContent(4  , 4.8716  	  ); 
  hd->SetBinContent(5  , 4.94311    ); 
  hd->SetBinContent(6  , 4.92585    ); 
  hd->SetBinContent(7  , 5.0351  	  ); 
  hd->SetBinContent(8  , 5.25363    ); 
  hd->SetBinContent(9  , 5.30603    ); 
  hd->SetBinContent(10 ,  5.44443   ); 
  hd->SetBinContent(11 ,  5.53738   ); 
  hd->SetBinContent(12 ,  5.54719   ); 
  hd->SetBinContent(13 ,  5.76489   ); 
  hd->SetBinContent(14 ,  5.94505   ); 
  hd->SetBinContent(15 ,  6.15808   ); 
  hd->SetBinContent(16 ,  6.43717   ); 
  hd->SetBinContent(17 ,  6.73969   ); 
  hd->SetBinContent(18 ,  6.9798    ); 
  hd->SetBinContent(19 ,  7.25404   ); 
  hd->SetBinContent(20 ,  7.69753   ); 
  hd->SetBinContent(21 ,  8.0544    ); 
  hd->SetBinContent(22 ,  8.73745   ); 
  hd->SetBinContent(23 ,  9.82025   ); 
  hd->SetBinContent(24 ,  11.7717   ); 
  hd->SetBinContent(25 ,  13.485    ); 
  hd->SetBinContent(26 ,  7.77343   ); 
  hd->SetBinContent(27 ,  1.36615   ); 
  hd->SetBinContent(28 ,  1.25825   ); 
  hd->SetBinContent(29 ,  0.556537  );  
  hd->SetBinContent(30 ,  0.316932  );  
  hd->SetBinContent(31 ,  0.348092  );  
  hd->SetBinContent(32 ,  0.472119  );  
  hd->SetBinContent(33 ,  0  	  ); 








  
  hd->Draw("");

  return 0;
}




//get dL/dz without normalisation
double dLdz(double zp, double R0, double birk) {
  double alpha=0.02543;
  double p=1.742;
  double dEdx_inv = p*pow(alpha,1./p)*pow((R0-zp),(1.-1./p));
  double val = 1./(dEdx_inv + birk);
    return val;
}


// Return term which includes dL/dz
double dLTerm(double *zp, double *par) {
  double pi = 3.14159265358979323846;
  double beta=0.0012;
  double R0     = par[0];
    double sigma  = par[1];
    double birk   = par[2];
    double z      = par[3];
    double preFactor = (1.+beta*(R0-zp[0]))/(1+beta*R0);
    double foldingTerm = 1./(sqrt(2*pi)*sigma)*exp(-(z-zp[0])*(z-zp[0])/(2*sigma*sigma));
    double val = preFactor * dLdz(zp[0],R0,birk) * foldingTerm;
    return val;
}


// Return term which includes L(z)
struct LTerm {
  double gamma_H2O = 0.6;
  double precision = 1e-7; // a parameter for numerical integration, reduce to speed up calculation
  double beta=0.0012;
  TF1* fdLdz = nullptr;

    LTerm(TF1 &f) : fdLdz(&f) {}

    double operator() (double *zp, double *par) {
        double R0   = par[0];
        double birk = par[1];

        fdLdz->SetParameters(R0,birk);
        double LofZ = 0;
        if(zp[0]<(1-precision)*R0) LofZ = fdLdz->Integral(zp[0],(1-precision)*R0,precision);

        double preFactor = gamma_H2O*beta/(1+beta*R0);
        double foldingTerm = 1.;//approximation without folding. It is absolutely identical to with folding since LTerm is pretty flat

        double val = preFactor * LofZ * foldingTerm;
        return val;
    }
};







//This returns the actual quenched Bragg fit curve
struct QuenchedBragg {
  double rho     = 1.0e-3;
  double S       = 10000;
  double precision = 1e-7; // a parameter for numerical integration, reduce to speed up calculation    
    std::unique_ptr<TF1> fdLTerm;
    std::unique_ptr<TF1> fLTerm;
    std::unique_ptr<TF1> fdLdz;
    
    QuenchedBragg() {
        fdLdz   = std::unique_ptr<TF1>(new TF1("fdLdz",dLdz,0.,1.,2));
        fdLTerm = std::unique_ptr<TF1>(new TF1("fdLTerm",dLTerm,0.,1.,4));
        fLTerm  = std::unique_ptr<TF1>(new TF1("fLTerm",LTerm(*fdLdz),0.,1.,2));
    }
      
    double operator() (double* z, double* par) {
        double R0      = par[0];
        double sigma   = par[1];
        double phi0    = par[2];
        double xOffset = par[3];
        double birks   = par[4];

        if (R0<=0.) return 0.;
        if (sigma<=0.) return 0.;
        if (phi0<=0.) return 0.;
        if (xOffset<0.) return 0.;
        if (birks<=0.) return 0.;
        if (R0 != R0) return 0.; //Check if R0 is nan
        if (z[0] != z[0]) return 0.; //Check if z is nan
        z[0] = z[0] + xOffset; //add xOffset
        if (z[0] < 0) return 0.;

        // Some numerical parameters. They are optimized in such a way that the fit returns accurate results without throwing segfaults or unreached precision errors
        double zDef = 5.*sigma; //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(z[0]-zDef,(1.-precision-1e-3)*R0); // integration start. Subtract constant to make sure a is smaller than b
        double b = min(z[0]+zDef,(1.-precision)*R0); // integration end

        fdLTerm->SetParameters(R0,sigma,birks,z[0]);
        double integral1 = fdLTerm->Integral(a,b,precision);
	;
        fLTerm->SetParameters(R0,birks);
        double integral2 = fLTerm->Eval(z[0]); // approximation without folding

        // Put all parts together and add scaling factor
        double preFactor = S*phi0/rho;
        double integral = integral1 + integral2;

        // Combine prefactor and integral of all parts
        double val = preFactor*integral;
        return val;
    }
};