//	filename: "bwfitpiplus02.C" 03/27/2008
//	add extra parameter for the integral constant (con)
//	from bwfitpiplus01.C
//
//	filename: "bwfit02_5.C"	03/21/2008 was wrong x[],y[] convert Pt, mt
//	modified from testfit04.C by moneta
//	copy from bwfit02_4.C
//	try pi+ particles only first.
//	(convert from pt to mt)

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

//#include "TGraph.h"
#include "TGraphErrors.h"

using namespace std;

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

Double_t radius;
TCanvas * c1 = 0; 
TGraph * gr = 0; 


void drawgraph(int npts, const double * x, const double * y, const double * errx, const double * erry)
{
//	creates a TGraph object to display the computed curves
	gr  = new TGraphErrors(npts, x, y, errx, erry);
	c1 = new TCanvas();

	gr->SetFillColor(19);
	gr->SetLineColor(4);
	gr->SetLineWidth(1);
	gr->SetMarkerColor(4);
	gr->SetMarkerStyle(21);
	gr->SetTitle("piplus.txt");
//	gr->SetXTitle("mt - m0 [GeV/c^2]");
//	gr->GetXaxis()->CenterTitle();

	gr->Draw("AC*");
	c1->cd();
	gPad->SetLogy();
	gPad->Modified();
	c1->Update();

}


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

//	define the vector in the routine (do not use globals, can have conflict with other definition of x !!!)
	const Int_t npts = 17;		// # of data
	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
*/

	Double_t mass = 0.1395679;			//	pi+ mass = 0.1395679 GeV
	for(Int_t i=0; i<npts; i++)			//	convert from Pt to (mt - m0)
	{
		x[i] = sqrt( x[i]*x[i] + mass*mass ) - mass;
	}

	cout << "\n\tNumber of data" << "\n\tindex = " << npts << endl << endl;

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

	drawgraph(npts, x, y, errx, erry); 
}


//	blastwave function from E.J.Kim "flow_c12.C" modified
Double_t bwfitfunc(Double_t *x, Double_t *par)
{
	//	*x is r in this case:	x[0]
	//	par[0] = alpha;
	//	par[1] = Tf;
	//	par[2] = beta;
	//	par[3] = con.;	:	integral constant is par 3
	//	par[4] = mt;	:	mt is par 4

	Double_t out;
	Double_t mass = 0.1395679;
	Double_t mt = par[4];
	Double_t pt = sqrt( mt*mt - mass*mass );

//	Double_t con = 1023.0; 
//	Double_t con = 55.0; 

	Double_t rho = TMath::ATanH(par[2]*pow(x[0]/13.0, par[0]));

	out = par[3] * x[0] * mt 
			* 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 integral of fFunc (mt, r) in r
	double operator() (double *x , double *p) 
	{
		// *x is mt
		// p[] is alpha, Tf, beta, and con
		double radius = 13.0;	//	radius = 13.0 fm (Rmax)
		std::copy(p, p + 4, param);
		param[4] = *x;    // set value of mt for integrand function (fFunc)
		return fFunc->Integral(0.0, radius, param);
	}

	TF1 *fFunc; // pointer to the integral function
//	double param[4]; 
	double param[5];
};


void fitexample()
{
	//	read the file, fill and  draw the graph
	readfile();

	//	write the fitting parameter values on the graph
	gStyle->SetOptFit(1111);

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

	MyIntegFunc *integ = new MyIntegFunc(funcx);

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

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

	ifuncx->SetParameter(0, 0.398);	//	alpha
	ifuncx->SetParameter(1, 0.118);	//	temp.
	ifuncx->SetParameter(2, 0.75);	//	beta
	ifuncx->SetParameter(3, 100.0);	//	con.
 
	//Need parameter limits. Otherwise I end up with negative Bessel functions.

	//	ifuncx->SetParLimits(0, 0.398, 0.398);	//	alpha
	ifuncx->SetParLimits(0, 0.390, 0.400);		//	alpha
	ifuncx->SetParLimits(1, 0.01, 0.2);			//	temp.
	ifuncx->SetParLimits(2, 0.2, 0.95);			//	beta
	ifuncx->SetParLimits(3, 10.0, 10000.0);		//	con.
//	ifuncx->SetParLimits(3, 100.0, 100.0);		//	con.

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

	ifuncx->SetLineColor(kBlue); 
	ifuncx->Draw("same");

	c1->Update();

	gr->Fit(ifuncx);

	gr->Draw("same");

	c1->Update();

	//	Gets the fitted parameters from minuit
	//	Double_t parameter, erro;

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

	for(Int_t n=0; n<5; 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;
	}
#endif