//	filename: "bwfit02_2.C"	03/12/2008
//	modified from testfit04.C by moneta
//	try pi+ particles only first.


#include <iostream>
#include <fstream>
#include "TF1.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"

#include "TGraphErrors.h"

using namespace std;

#ifdef __CINT__
//gROOT -> Reset();
#endif

Double_t radius;
TCanvas * c1 = 0; 

//	TF2 *integralfunc;


Double_t x[50], y[50], errx[50], erry[50], par[10], err[10];
Int_t npts = 0;		// # of data

void readfile()
{
//	fills up the x and y reading file
	ifstream infile("piplus.txt");

	Int_t index = 0;
//	while(index <= 50)
	while(!infile.eof())
	{
		infile >> x[index] >> y[index] >> errx[index] >> erry[index];
		index++;
	}

	infile.close();
	index--;
	cout << "\n\tNumber of data" << "\n\tindex = " << index << endl << endl;

	npts = index;
	for(Int_t i=0; i<index; i++)
	{	
		cout << "\t" << x[i] << "\t" << y[i] << "\t\t" << errx[i] << "\t" << erry[i] << endl; 
	}
	cout << endl;
}


void drawgraph()
{
//	creates a TGraph object to display the computed curves

	TGraph *gr  = new TGraphErrors(npts, x, y, errx, erry);

	gr->SetFillColor(19);
	gr->SetLineColor(4);
	gr->SetLineWidth(1);
	gr->SetMarkerColor(4);
	gr->SetMarkerStyle(21);
	gr->SetTitle("piplus.txt");
	gr->Draw("AC*");

	c1->cd();
	gPad->SetLogy();
	gPad->Modified();

//	TF1 *fit = new TF1("fit", "gaus", 0.0, 8.0);

//	cout << "\t" << "par[0]" << "  " << "par[1]" << "  " << "par[2]" << endl;
//	cout << "\t" << par[0] << "  " << par[1] << "  " << par[2] << endl;

//	fit->SetParameters( par[0], par[1], par[2]);
//	fit->Draw("same");

}


//	blastwave function from E.J.Kim "flow_c12.C"
Double_t bwfitfunc(Double_t *x, Double_t *par)
{
	//	*x is Pt in this case.
	//	par[0] = alpha;
	//	par[1] = Tf;
	//	par[2] = beta;

	Double_t arg = 0.0;
	Double_t out;

	Double_t mass = 0.1395679;
	Double_t mt = x[0] + mass;  
	Double_t con = 1023.0; 

	Double_t rho = TMath::ATanH(par[2]*pow(x[0]/13.0, par[0]));
	Double_t pt = sqrt(mt*mt - mass*mass);
//	Double_t temp = param[5];

	out = con * x[0]
			* TMath::BesselI0(pt * TMath::SinH(rho)/par[1])
			* TMath::BesselK1(mt * TMath::CosH(rho)/par[1]);
	return out;
}



/*
//	structure representing the integral of a function between 0 and radius
struct MyIntegFunc()
{
	//	constructor using the TF1 pointer
	MyIntegFunc(TF1 *f):
		fFunc(f){}

	//	evaluate the function (is it independent of x ???)
	double operator() (double *x, double *p) const
	{
		double radius = 13.0;	//	radius = 13.0 fm (Rmax)
		return fFunc->Integral(a, radius, p);
	}

	TF1 *fFunc; // pointer to the integral function
};


//	integral function
TF1 *funcx = new TF1("funcx", bwfitfunc, 0.0, radius+2.0, 3);


// integral function of funcx
TF1 *ifuncx = new TF1("ifuncx", new MyIntegFunc(funcx), 0.0, radius +2.0, 3, "MyIntegFunc");
*/




/*
//	need the integral for bwfitfunc
Double_t integralfunc(Double_t *x, Double_t *par)
{
	Double_t integral;
	Double_t radius = 13.0;	//	radius = 13.0 fm (Rmax)	


	TF1 *funcx = new TF1("funcx", bwfitfunc, 0.0, radius+2.0, 3);

	funcx->SetParameter(2, 0.75);	//	beta
	funcx->SetParameter(1, 0.118);	//	temp.
	funcx->SetParameter(0, 0.398);	//	alpha

	funcx->SetParLimits(2, 0.2, 0.95);	//	beta
	funcx->SetParLimits(1, 0.01, 0.2);	//	temp.
	funcx->SetParLimits(0, 0.390, 0.400);	//	alpha

//	cout << "\nin the integralfunc" << endl;
//	cout << "par[0] = " << par[0] << endl;
//	cout << "par[1] = " << par[1] << endl;
//	cout << "par[2] = " << par[2] << endl;
	

	integral = funcx->Integral(0.0, radius);
	if(integral <= 0.0) integral = 0.00000001;

	return integral;
}
*/


//	structure representing the integral of a function between 0 and radius
struct MyIntegFunc
{
	//	constructor using the TF1 pointer
	MyIntegFunc(TF1 *f):
		fFunc(f) {}

	//	evaluate the function (is it independent of x ???)
	double operator() (double * /*x */, double *p) const
	{
		double radius = 13.0;	//	radius = 13.0 fm (Rmax)
		return fFunc->Integral(0., radius, p);
	}

   TF1 *fFunc; // pointer to the integral function
};




void fitexample()
{
   readfile();
   

//	TF1 *func = new TF1("fit", fitf, -3, 3, 3);
//	TF1 *func = new TF1("bwfitfunc", bwfitfunc, 0.0, 3.0, 3);
//	TF1 *func = new TF1("integralfunc", integralfunc, 0.0, 3.0, 3);

//	func->SetParameters(5.0, 3.5, 1.5);
//	func->SetParNames("Constant", "Mean_value", "Sigma");


//	creates a TGraph object to display the computed curves

	TGraph *gr  = new TGraphErrors(npts, x, y, errx, erry);

//	gr->SetXaxis()->SetLimits(0.0, 3.0);
	
	gr->SetFillColor(19);
	gr->SetLineColor(4);
	gr->SetLineWidth(1);
	gr->SetMarkerColor(4);
	gr->SetMarkerStyle(21);
	gr->SetTitle("piplus.txt");
	gr->Draw("ACP");

        c1 = new TCanvas("c1");

	c1->cd();
	gPad->SetLogy();
	gPad->Modified();


//	============================================================


	gStyle->SetOptFit(1111);
//	gr->Fit("gaus");



	//	have to set the parameters limits.
//	++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/*
	func->SetParameter(2, 0.75);	//	beta
	func->SetParameter(1, 0.118);	//	temp.
	func->SetParameter(0, 0.398);	//	alpha
 
	//Need parameter limits. Otherwise I end up with negative Bessel
	//functions.
	func->SetParLimits(2, 0.2, 0.95);	//	beta
	func->SetParLimits(1, 0.01, 0.2);	//	temp.
//	func->SetParLimits(0, 0.398, 0.398);	//	alpha
	func->SetParLimits(0, 0.390, 0.400);	//	alpha

//	gr->Fit("bwfitfunc");
//	gr->Fit("integralfunc");
*/

//	++++++++++++++++++++++++++++++++++++++++++++++++++++++++



//	integral function
        TF1 *funcx = new TF1("funcx", bwfitfunc, 0.0, radius+2.0, 3);



        MyIntegFunc *integ = new MyIntegFunc(funcx);

// integral function of funcx
//	TF1 *ifuncx = new TF1("ifuncx", new MyIntegFunc(funcx), 0.0, radius +2.0, 3, "MyIntegFunc");

        TF1 *ifuncx = new TF1("ifuncx", integ,  0.0, radius +2.0, 3, "MyIntegFunc");


//	++++++++++++++++++++++++++++++++++++++++++++++++++++++++





	//	have to set the parameters limits.
//	++++++++++++++++++++++++++++++++++++++++++++++++++++++++

	ifuncx->SetParameter(2, 0.75);	//	beta
	ifuncx->SetParameter(1, 0.118);	//	temp.
	ifuncx->SetParameter(0, 0.398);	//	alpha
 
	//Need parameter limits. Otherwise I end up with negative Bessel
	//functions.
	ifuncx->SetParLimits(2, 0.2, 0.95);	//	beta
	ifuncx->SetParLimits(1, 0.01, 0.2);	//	temp.
//	func->SetParLimits(0, 0.398, 0.398);	//	alpha
	ifuncx->SetParLimits(0, 0.390, 0.400);	//	alpha

//	++++++++++++++++++++++++++++++++++++++++++++++++++++++++





	gr->Fit(ifuncx);

/*
	TMinuit *gMinuit = new TMinuit(3);

	gMinuit->SetFCN("integralfunc");
	Double_t arglist[4];
	Int_t ierflg = 0;

	arglist[0] = 1;
	gMinuit->mnexcm("SET ERR", arglist,1 ,ierflg);

	Double_t vstart[3] = {0.398, 0.118, 0.75};
	Double_t step[3] = {0.00001, 0.00001, 0.0001};
	gMinuit->mnparm(0, "alpha", vstart[0], step[0], 0, 0, ierflg);
	gMinuit->mnparm(1, "temp", vstart[1], step[1], 0, 0, ierflg);
	gMinuit->mnparm(2, "beta", vstart[2], step[2], 0, 0, ierflg);

// Now ready for minimization step
	arglist[0] = 1.0;
	arglist[1] = 1.0;
	gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);

// Print results
	Double_t amin,edm,errdef;
	Int_t nvpar,nparx,icstat;
	gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
	gMinuit->mnprin(3, amin);
*/

	gr->Draw("same");



	c1->Update();
/*
	for(Int_t i=0; i<3; i++)
	{
		cout << "\tpar[" << i << "] = " << par[i] << endl;
	}
*/
//	Gets the fitted parameters from minuit
	Double_t parameter, erro;

	for (Int_t n=0 ; n<3 ; n++)
    {
//		gMinuit->GetParameter(n, parameter, erro) ;
       par[n] = ifuncx->GetParameter(n) ;
       err[n] = ifuncx->GetParError(n) ;
    }

	cout << "In the fitexample()" << endl;

//	cout << "\t" << "par[0]" << "  " << "par[1]" << "  " << "par[2]" << endl;

	for(int n=0; n<3; n++)
	{
		cout << "\t" << "par[" << n << "] = " << par[n] << " +- " << err[n] << endl;
	}

//	testgminuit();

}


#ifdef USE_GMINUIT

//void  fcn(Double_t &f, Double_t *par)
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
//	calculate chisquare
	Double_t chisq = 0.0;
	Double_t delta = 0.0;

//	cout << "\t" << "par[0]" << "  " << "par[1]" << "  " << "par[2]" << endl;

	for (Int_t i=0; i<npts; i++)
	{
           //  what is fitf ? need to be defined _ comment next line
           //	delta  = (y[i] - fitf(x[i], par));
		chisq += delta*delta;
	}

//	cout << "\tchisq = " << chisq << endl;
	f = chisq;
}


void testgminuit()
{
	cout << "\n\nIn the testgminuit()" << endl;

	
	for(Int_t i=0; i<npts; i++)
	{	
		cout << "\t" << x[i] << "\t" << y[i] << "\t" << errx[i] << "\t" << erry[i] << endl;
	}

//	TMinuit *gMinuit = new TMinuit(4);  //initialize TMinuit with a maximum of 5 params
	TMinuit *gMinuit = new TMinuit(3);

	gMinuit->SetFCN(fcn);
	Double_t arglist[4];
	Int_t ierflg = 0;

	arglist[0] = 1;
	gMinuit->mnexcm("SET ERR", arglist,1 ,ierflg);

	static Double_t vstart[4] = {5.0, 3.0, 1.0};
	static Double_t step[4] = {0.001, 0.001, 0.001};
	gMinuit->mnparm(0, "a1", vstart[0], step[0], 0, 0, ierflg);
	gMinuit->mnparm(1, "a2", vstart[1], step[1], 0, 0, ierflg);
	gMinuit->mnparm(2, "a3", vstart[2], step[2], 0, 0, ierflg);

//   gMinuit->mnparm(3, "a4", vstart[3], step[3], 0,0,ierflg);

// Now ready for minimization step
	arglist[0] = 5000;
	arglist[1] = 1.0;
	gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);

// Print results
	Double_t amin,edm,errdef;
	Int_t nvpar,nparx,icstat;
	gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
	gMinuit->mnprin(3, amin);

 // Gets the fitted parameters from minuit
	Double_t parameter, erro;
	for (Int_t n=0 ; n<3 ; n++)
    {
		gMinuit->GetParameter(n, parameter, erro) ;
		par[n] = parameter ;
		err[n] = erro;
    }

	cout << "\t" << "par[0]" << "  " << "par[1]" << "  " << "par[2]" << endl;
	cout << "\t" << par[0] << "\t  " << par[1] << "\t  " << par[2] << endl;
/*
	cout << "\t cnt = " << cnt << endl;

//	creates a TGraph object to display the computed curves
	TGraph *gr  = new TGraphErrors(npts, x, y, errx, erry);

	gr->SetMinimum(-2.0);
	gr->SetMaximum(8.0);
//	gr->GetXaxis()->SetRangeUser(-2.0, 10.0);
	gr->GetXaxis()->SetLimits(-2.0, 10.0);

	gr->SetFillColor(19);
	gr->SetLineColor(4);
	gr->SetLineWidth(1);
	gr->SetMarkerColor(4);
	gr->SetMarkerStyle(4);
	gr->SetTitle("testdata01.txt");
	gr->Draw("ACP");

	TF1 *fit = new TF1("fit", "gaus", 0.0, 8.0);

	for (Int_t i = 0; i < 3; i++)
	{
		gMinuit->GetParameter(i, par[i], err[i]);
		cout << "Parameter " << i << " = " << par[i] << " +- " << err[i] << endl;
	}

	cout << "\t" << "par[0]" << "  " << "par[1]" << "  " << "par[2]" << endl;
	cout << "\t" << par[0] << "  " << par[1] << "  " << par[2] << endl;

//	fit->SetParameters( 4.93835e+00, 3.17093e+00, 1.29764e+00);
	fit->SetParameters( par[0], par[1], par[2]);
	fit->Draw("same");
*/

//	gStyle->SetOptFit();
//	gr->Fit("gaus");

//	c1->Update();

//	cout << "\tchisq = " << f << endl;

}
#endif