#include <fstream.h>
#include <iostream.h>
#include <stdlib.h>
#include <TROOT.h>
#include <stdio.h>
#include <TGraph.h>
#include <TF2.h>
#include <TF3.h>
#include <TFile.h>
#include <TMath.h>
#include <TGraph2DErrors.h>
#include <TCanvas.h>

#define data_length 12  // number of points
#define Mq 0.01
#define Mc 1.40

// ****************************************** Global variables **************************************

double x_glob = 0.0, Q_glob = 0.0, x0 =0.02, omega_s = 0.0;

// **************************************************************************************************

#define R       0.89
#define c       ( (2.828*b)/R )
#define f		 4.0
#define b_0		 ((11-(2.0/3)*f))
#define gamma	(12.0/b_0)
#define nu		0.5
#define file_DIS	"Fit_DIS_1.txt"

#define aa    	4.75

#define a		TMath::Sqrt( (qu2)*z*(1-z)+mq*mq )
#define K12 	TMath::Power(TMath::BesselK1(a*r*aa),2)
#define K02 	TMath::Power(TMath::BesselK0(a*r*aa),2)
#define Jr		2*TMath::Pi()*r
#define Jb		2*TMath::Pi()*b 

#define r_min 	0
#define r_max  	50				// limits of triple integration (r,b,z)
#define b_min 	0
#define b_max 	20
#define z_min 	0
#define z_max 	1	

// *********************************** Internal functions of integration that work properly *********************
double S(double b){return c*TMath::BesselK1(c);}

double psiphq(double qu2, double r, double z, double mq, double ch)
{
	return (  ((3*ch*ch)/(137*2*TMath::Power(TMath::Pi(),2)))* ( (z*z + (1-z)*(1-z))*a*a*K12 + mq*mq*K02 + 
4*qu2*z*z*(1-z)*(1-z)*K02 ) );
}

double photonwf(double qu2, double r, double z)
{
	return ( psiphq(qu2,r,z,Mc,2.0/3) +
			 psiphq(qu2,r,z,Mq,2.0/3) +
			 2*psiphq(qu2,r,z,Mq,1.0/3) 
			);
}

double BesselI2(float z){return (TMath::BesselI0(z) - 2*TMath::BesselI1(z)/z);}
double BesselI3(float z){return (TMath::BesselI1(z) - 4*BesselI2(z)/z);}

double o_less_s(float dif_Y, float t, float omega_s, double omega_o, double gam, double A1)
{
	float arg = 2*TMath::Sqrt(t*dif_Y);

	return (gam*(TMath::BesselI0(arg) 
		 + TMath::BesselI1(arg)*(omega_o/omega_s)  
		 + BesselI2(arg)*TMath::Power(omega_o/omega_s,2)         
		 + BesselI3(arg)*TMath::Power(omega_o/omega_s,3))
		 + A1*(TMath::BesselI0(arg)/(omega_s + nu) + TMath::Exp(-(t/nu + nu*dif_Y)) )
		   );
}
double o_larger_s(float dif_Y, float t, float omega_s, double omega_o, double gam, double A1)
{
	float arg = 2*TMath::Sqrt(t*dif_Y);

	return (gam*(TMath::Exp(omega_o*dif_Y + t/omega_o)
		 - TMath::BesselI1(arg)*(omega_s/omega_o)  
		 - BesselI2(arg)*TMath::Power(omega_s/omega_o,2)
		 - BesselI3(arg)*TMath::Power(omega_s/omega_o,3))
	 	 + A1*(TMath::BesselI0(arg)/(omega_s + nu) + TMath::Exp(-(t/nu + nu*dif_Y)))
		   );
}     

double xG(double r,double x, double MU0, double omega_o, double gam, double A1, double sss)
{
	double Y   = TMath::Log(1/x);
	double Y_0 = TMath::Log(1/x0);
	double dif_Y = (Y-Y_0);

	double w     = TMath::Log( (Q_glob + MU0)/0.04 );	
	double w_0   = TMath::Log( MU0/0.04 );
	double t	 = gamma*TMath::Log(w/w_0)/2;
	
	omega_s		 = TMath::Sqrt(t/dif_Y);

	double pref  = TMath::Exp(-t);

	if (dif_Y == 0) return 1;

	if (omega_s >= omega_o)
		return pref*o_less_s(dif_Y, t, omega_s, omega_o, gam, A1);
	
	if (omega_s < omega_o)
		return pref*o_larger_s(dif_Y, t, omega_s, omega_o, gam, A1);
}

double alpha_s(double r, double MU0, double omega_o, double A1, double gam, double sss)
{
	return (4*TMath::Pi())/(2*b_0*TMath::Log( (Q_glob + MU0)/0.04 )); 
}


double xG_Y_0(double r, double MU0, double omega_o, double gam, double A1)
{

	double Y_0			= TMath::Log(1/x0);
	double w		    = TMath::Log( (Q_glob + MU0)/0.04 );
	double w_0			= TMath::Log( MU0/0.04 );
	double t			= gamma*TMath::Log(w/w_0)/2;
	double pref			= TMath::Exp(-t);
	
	double omega_S		= TMath::Sqrt(t/Y_0);


	if (omega_s >= omega_o)
		return pref*o_less_s(Y_0, t, omega_S, omega_o, gam, A1);
				
	if (omega_s < omega_o)
		return pref*o_larger_s(Y_0, t, omega_S, omega_o, gam,A1);
}

double amplitude(double r, double b, double MU0, double omega_o,double gam, double A1, double sss)
{
	double xm   = x_glob*(1 + (4*Mc*Mc)/(Q_glob) );
	
	double amp  = r*r*alpha_s(r, MU0, omega_o, A1, gam, sss)*xG_Y_0(r, MU0, omega_o, gam, A1)
		*xG(r, xm, MU0, omega_o, gam, A1, sss)*S(b);
		
		
	return ( amp/(amp + 1) );

}


// *************************************** Integrtion (z,r,b,) **********************************************

double F2DIS(double* var, double* par)
{
		double r  = var[0];
		double z  = var[1];
		double b  = var[2];

		double q2  = Q_glob;
		
		double MU0		= par[0];
		double omega_o	= par[1];
		double gam		= par[2];
		double A1		= par[3];
		double xuy      = par[4];
		
		
		return ( (q2*TMath::Power(aa,4)*137)/(2*TMath::Power(TMath::Pi(),2)) )*Jr*Jb*photonwf(q2, r, z)*amplitude(r,b, MU0, omega_o, gam, A1,xuy);
	
}


double integrated_F2(double MU0, double omega_o, double gam, double A1, double sss)
{
	double F2dis = 0.0;

	TF3 *dis = new TF3("struct", F2DIS, r_min,r_max, z_min, z_max, b_min, b_max,5);
	dis->SetParameters(MU0, omega_o, gam, A1, sss);
	return (F2dis = dis->Integral(r_min, r_max, z_min, z_max,b_min, b_max, 1e-20));
}

double fitt(double* x, double* par)
{
	x_glob = x[0];    // global variables defined above used for inserting the values x,y in the fit
	Q_glob = x[1];
	return integrated_F2(par[0], par[1], par[2], par[3], par[4]); // parameters to be fitted
}



// ********************************** Main ********************************************

int main()
{

	double param[5] = {1.0, 1.0, 1.0, 1.0, 1.0}; // initial values of fitted parameters: MU0  omega_o  gam  A1 sss
	
	float Q2[data_length+1], x_b[data_length+1], F2[data_length+1], 
		er_x[data_length+1], er_y[data_length+1], er_z[data_length+1];


	float q2,xb,f2,er1,er2,Er;

	int file_length = 0;

	FILE * fdata;
    fdata = fopen ("data.txt","r");   // reading the data points from the file
	rewind (fdata);

// *************************************** Loading the data points *************

    while (!feof(fdata)) {
	fscanf (fdata, "%f %e %f %f %f", &q2, &xb, &f2, &er1, &er2);
	Er = TMath::Sqrt(er1*er1 + er2*er2);

	Q2[file_length]	 = q2;
	x_b[file_length] = xb;    // x, y, z, values
	F2[file_length]	 = f2;
	
	er_x[file_length] = 0.0;
	er_y[file_length] = 0.0;  // their errors, errors are only for value z (f2)
    er_z[file_length] = Er;
	file_length++;
	}

// ****************************************************************************

	TGraph2DErrors *dte = new TGraph2DErrors(data_length);
	for(int i=0;i<data_length;i++){
	dte->SetPoint(i,Q2[i],x_b[i],F2[i]);					// loading the points to TGraph
	dte->SetPointError(i,er_x[i],er_x[i],er_x[i]);
	};

	// output, to check whether it works
	for(int i=0;i<data_length;i++){
	cout<<" Q2 "<<Q2[i]<<" x "<<x_b[i]<<" f2 "<<F2[i]<<" er x "<<er_x[i]<<" er y "<<er_y[i]<<" er f2 "<<er_z[i]<<endl;}
	cout<<endl<<endl;	


	TF2 *func = new TF2("fit", fitt, 0,10, 0,10,5); // <- If you change the number of parameters to 5-> no convergence 
													//(you should wait a couple of minutes)
													// for numb_par = 1 it is converged but still doesn't work.		

	func->SetParName(0,"MU0");
	func->SetParName(1,"Omega_0");
	func->SetParName(2,"gam");							// setting  the parameters
	func->SetParName(3,"A1");
	func->SetParName(4,"sss");


	func->SetParameters(param);
	dte->Fit("fit","r");								// the fit
	

return 0;

}