//	filename: "bwfit02_3.C"	03/12/2008
//	modified from testfit04.C by moneta
//	try pi+ particles only first.
//	changing the way of reading data

#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--;
*/

	Int_t index = 17;

	Double_t x[17] = {0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85};
	Double_t y[17] = {277.15, 190.05, 112.20, 70.24, 45.53, 29.16, 19.09, 12.84, 8.54, 5.73, 3.99, 2.76, 1.84, 1.23, 0.92, 0.66, 0.48};
	Double_t errx[17] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
	Double_t erry[17] = {6.27, 2.72, 1.36, 0.75, 0.46, 0.32, 0.23, 0.17, 0.13, 0.10, 0.08, 0.06, 0.05, 0.04, 0.03, 0.03, 0.02};

//	data from "piplus.txt" file	//
/*	  x			y		errx	erry
	0.25	277.15		0.0		6.27
	0.35	190.05		0.0		2.72
	0.45	112.20		0.0		1.36
	0.55	 70.24		0.0		0.75
	0.65	 45.53		0.0		0.46
	0.75	 29.16		0.0		0.32
	0.85	 19.09		0.0		0.23
	0.95	 12.84		0.0		0.17
	1.05	  8.54		0.0		0.13
	1.15	  5.73		0.0		0.10
	1.25	  3.99		0.0		0.08
	1.35	  2.76		0.0		0.06
	1.45	  1.84		0.0		0.05
	1.55	  1.23		0.0		0.04
	1.65	  0.92		0.0		0.03
	1.75	  0.66		0.0		0.03
	1.85	  0.48		0.0		0.02
*/

	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