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

	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"
Double_t bwfitfunc(Double_t *x, Double_t *par)
{
   
        // *x is r in this case
        // pt is par 3
	//	par[0] = alpha;
	//	par[1] = Tf;
	//	par[2] = beta;
        //      par[3] = pt

	//Double_t arg = 0.0;
	Double_t out;
        Double_t pt = par[3]; 
	Double_t mass = 0.1395679;
        Double_t mt = sqrt( mass*mass + pt*pt);

	Double_t con = 1023.0; 

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

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

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




void fitexample()
{

   // read the file, fill and  draw the graph
   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




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


	gStyle->SetOptFit(1111);



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

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



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



        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

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


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


	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;

        double par[3]; 
        double err[3]; 
	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