//////////////////////////////////////////////////////////////////////////////////////////////////////
///
/// Ashley Rose, June 2008
///
//////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>

using namespace std;

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TStyle.h"
#include <math.h>


// int oldmain()
// {
//   //gaussian fit
//   gROOT->SetBatch();

//   //TFile *file_op = new TFile("/twist/data28/olin/tzero_2007/split1/set84/set84.root");
//   //TFile *file_op = new TFile("/home/arose/e614soft/triumf/mofia/devel/run/root/MC_IonClusterTest/001001.root");
//   TFile *file_op = new TFile("/home/arose/e614soft/triumf/mofia/ICT/run/root/MC_IonCluster_Tuning/008216.root");


//   //TH1D *H = (TH1D*)file_op->Get("H511096");
//   //TH1D *H = (TH1D*)file_op->Get("H99996");
//     TH1D *H = (TH1D*)file_op->Get("H09052");
//   H->Fit("gaus");
//   double xbar = H->GetMean();
//   cout<<H<<endl;
//   TCanvas *c1 = new TCanvas();
//   cout<<"The mean of the first gaussian is "<<xbar<<endl;

//   H->Draw();
//   c1->Print("GausFitOutput.ps","landscape");

// }

// special fit


//define function
double fitf( double *x, double *par)
{

		//_______Art________//
		//   double result = 0.0;
		//   double xx = x[0];
		//   result = xx-par[0];
		//   result /= 1.4142*par[1];
		//   result = erf(result)+1.0;
		//   result *= 0.5*par[2];
		//   if( xx >= par[3] )
		//     {
		//       result *= exp((par[3]-xx)/par[4]);
		//     }
		//   return result;
		// };

		//   double fitglen( double *x, double *par )
		//   //______Glen_______//

		//   double result = 0.0;
		//   double xx = x[0];
		//   result = xx-par[0];
		//   result /= 1.4142*par[1];
		//   result = erf(result)+1.0;
		//   result *= 0.5*par[2];
		//   par[3] = par[0];      //testing

		//   if( xx >= par[3] )
		//     {
		//       result *= exp((par[3]-xx)/par[4]);
		//     }
		//   return result;
		// };

		// double fitroman( double *x, double *par )
		// {
		//________Roman _______//
	double result, t, a, b, c, d, f, h;   //b = mean, c = std, a = normalisation cst
	b = par[0];  
	c = par[1];
	a = par[2];
	h = par[3];
	
	d = a*exp(-0.5*((h-b)/c)*((h-b)/c));
	f = c*c/(h-b);
	
	t = x[0];
	
	if ( t <= h )
			result = a*exp(-0.5*((t-b)/c)*((t-b)/c));
	if ( t > h )
	result = d*exp(-1*(t-h)/f);
	return result;
	};

//use the function to fit the histogram
int main(int argc, char **argv)
{
	gROOT->SetBatch();

	//open file and get histogram:
	//TFile *f = new TFile("/twist/data28/olin/tzero_2007/split1/set84/set84.root");
	TFile *file_op;
	if (argc < 2 )  
			file_op = new TFile("/home/arose/e614soft/triumf/mofia/devel/run/root/Data_IonClusterTest/051000.root");
	else
			file_op = new TFile(argv[1]);

	cout<<file_op<<endl;

	TFile *Out = new TFile("Output.root","RECREATE");

	//First fit the gaussian
	TH1D *H2;
	if (argc < 3)
			H2 = (TH1D*)file_op->Get("H99996");
	else
			H2 = (TH1D*)file_op->Get(argv[2]);
	cout<<H2<<endl;

	H2->Fit("gaus");

	//Extraction of the parameters
	double peakcont = H2->GetMaximum();
	int peakxbin = H2->GetMaximumBin();
	double xval, sigma, chi2, xFWHM, FWHM; 
	double xmin, xmax;


	xval = H2->GetBinCenter(peakxbin); //this will be the "mean"
	xmin = xval - 7.0;  
	xmax = xval + 10.0; 
	H2->SetAxisRange(xmin, xmax);

	cout<<"The bin with the most entries is "<<peakxbin<<endl;
	cout<<"xval (pseudomean) is "<<xval<<endl;

	//Pick a fit to use with above parameters
	//  TF1 *func = new TF1("mysuperfit", fitf, xmin, xmax, 5);  // ART's 5 par
	TF1 *func = new TF1("mysuperfit", fitf, xmin, xmax, 4);  //ROMAN's 4 par

	sigma = H2->GetRMS();
	//double sigma = H2->GetRMS()/sqrt(H2->Integral());
	// FWHM = func->GetX(func, xFWHM, xFWHM);
	//cout<<"The full width half maximum is "<<FWHM<<endl;
	//func->SetParLimits(0, xval-5, xval+5);

	// func->SetParameters( xval, sigma, peakcont, xval+10, 5.171); //for ART's fit
	func->SetParameters( xval, sigma, peakcont, xval+10); //for ROMAN's fit

	//func->FixParameter(4, 5.171); 

	//fit the gaussian
	//H2->Fit( "gaus","","", xmin, xmax);  
	//TCanvas *c3 = new TCanvas();  //plot of fit
	//H2->Draw();
	//c3->Print("GausFitOutput2.ps","landscape");

	H2->Fit( func,"","", xmin, xmax);
	xFWHM = func->GetParameter(0) - (func->GetParameter(1) * 1.1775);
	cout<<"The x value at the full width half maximum is "<<xFWHM<<endl;

	cout<<"sigma = "<<func->GetParameter(0)<<endl;
	cout<<"height = "<<func->GetParameter(1)<<endl;

	chi2 = (func->GetChisquare())/(func->GetNDF());

	TStyle* gStyle = new TStyle("gStyle","gStyle");
	gStyle->SetOptFit(111);
	gStyle->cd();
	gStyle->SetLegendBorderSize(30);

	H2->UseCurrentStyle();
	func->SetParNames("Mean", "Sigma", "Peak", "10");

	Out->cd();

	H2->Write();

	// //plot of fit function
	// TCanvas *c2 = new TCanvas();
	// H2->SetXTitle("Time (ns)");
	// H2->SetYTitle("Entries");
	// H2->Draw();
	// c2->Print("FitOutput.ps","landscape");


}


//make && ./GausFit "/twist/data28/olin/tzero_2007/split1/set84/set84.root" "H511096"