//	filename: "bwfit03.C"	03/07/2008
//	modified from testfit04.C and bwfit02.C (moneta)
//	try pi+ particles only first.


#include <iostream.h>
#include <fstream.h>
#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
};


void fitexample()
{
	readfile();

//	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");


//	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", 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.
	ifuncx->SetParLimits(0, 0.390, 0.400);	//	alpha

	gr->Fit("ifuncx");
	gr->Draw("same");

	c1->Update();

//	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;
    }

	for(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++)
	{
		delta  = (y[i] - bwfitfunc(x[i], par));
		chisq += delta*delta;
	}

	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;
}
#endif