#include "NewTypes.h"

//Functions for the fit
double normal(double *x, double *par); //Normal
double parabolic(double *x, double *par); //Parabolic
double double_gaussian(double *x, double *par); //Double Gaussian

int main(int argc, char **argv)
{	
	//Create the application
	TApplication myApp("myApp", &argc, argv);
	
	//Display options
        gStyle->SetOptFit(1111);
	
	//Create the histograms
	TH1F histNormal("histNormal", "Gaussian", 100, -5., 5.);
	TH1F histParabolic("histParabolic", "Parabolic", 100, -2., 2.);
	TH1F histDoubleGaussian("histDoubleGaussian", "Double Gaussian", 100, -6., 6.);
	
	//Fill the histograms with 10000 entries
	for(int i(0); i < 10000; ++i)
	{
		//Gaussian distribution
		double numberGaussian(0.);
		rand_normal(numberGaussian, numberGaussian);
		histNormal.Fill(numberGaussian);
		
		//Parabolic distribution
		double numberParabolic(0.);
		rand_parabolic(numberParabolic, numberParabolic);
		histParabolic.Fill(numberParabolic);
		
		//Double Gaussian distribution
		double numberDoubleGaussian(0.);
		rand_double_gauss(numberDoubleGaussian, numberDoubleGaussian);
		histDoubleGaussian.Fill(numberDoubleGaussian);
	}
	
	//Fit the histograms
	TF1 fitFunctionNormal("fitFunctionNormal", normal, -5., 5., 3);
	fitFunctionNormal.SetParNames("N","m","#sigma");
	fitFunctionNormal.SetParameters(200, histNormal.GetMean(), 1);
	TF1 fitFunctionParabolic("fitFunctionParabolic", parabolic, -1., 1., 3);
	fitFunctionParabolic.SetParNames("N","m","#sigma");
	fitFunctionParabolic.SetParameters(200, histParabolic.GetMean(), 1);
	TF1 fitFunctionDoubleGaussian("fitFunctionDoubleGaussian", double_gaussian, -6., 6., 4);
	fitFunctionDoubleGaussian.SetParNames("N","m","#sigma","d");
	fitFunctionDoubleGaussian.SetParameters(200, histParabolic.GetMean(), 1,1);
	
	//Draw the histograms
	TCanvas *canvasDistributions = new TCanvas("canvasDistributions","Distributions",0,0,600,600);
	canvasDistributions->Divide(2,2);
	
	//Gaussian
	canvasDistributions->cd(1);
	histNormal.Draw();
	histNormal.Fit(&fitFunctionNormal,"RQ","");
	
	//Parabolic
	canvasDistributions->cd(2);
	histParabolic.Draw();
	histParabolic.Fit(&fitFunctionParabolic,"RQ","");

	//Double Gaussian
	canvasDistributions->cd(3);
	histDoubleGaussian.Draw();
	histDoubleGaussian.Fit(&fitFunctionDoubleGaussian,"RQ","");
		
	//Run the application
	myApp.Run();
	
	//Delete canvas
	delete canvasDistributions;
	
	return 0;
}

//Functions for the fit
double normal(double *x, double *par) //Normal
{
	double N = par[0]; //Normalization constant
	double m = par[1]; //Mean
	double s = par[2]; //Sigma
	double x_ = x[0];
   
   	if(s != 0.)
   	{
		return (N/(sqrt(2.*M_PI)*s)) * exp(-((x_-m)*(x_-m))/(2.*s*s));
	}
	else
	{
		return 1e30;
	}
}

double parabolic(double *x, double *par) //Parabolic
{
	double N = par[0]; //Normalization constant
	double m = par[1]; //Mean
	double s = par[2]; //Sigma
	double x_ = x[0];
    
	return N*(-(x_-m)*(x_-m)+s*s);
}

double double_gaussian(double *x, double *par) //Double Gaussian
{
	double N = par[0]; //Normalization constant
	double m = par[1]; //Mean
	double s = par[2]; //Sigma
	double d = par[3]; //Disatnce betwwen the peaks
	double x_ = x[0];
   
   	if(s != 0.)
   	{
   		double gaussianLeft((N/(sqrt(2.*M_PI)*s)) * exp(-((x_-(m-d/2.))*(x_-(m-d/2.)))/(2.*s*s)));
   		double gaussianRight((N/(sqrt(2.*M_PI)*s)) * exp(-((x_-(m+d/2.))*(x_-(m+d/2.)))/(2.*s*s)));
   		
		return gaussianLeft + gaussianRight;
	}
	else
	{
		return 1e30;
	}
}