#include "NewTypes.h"

//Define PI for Windows
double M_PI = 3.14159265;

//Global variables for the simulation
double QS = 0.; 
double QMOD = 0.005;
double BX = 1.0; //The x-beta function
double BY = 1.0; //The y-beta function

//Create a generator of random numbers
TRandom3 random_number = 0;

//Matrix multiplication
Vector MV(const Matrix & matrix, const Vector & vector)
{
	unsigned int size(vector.size());
	Vector result(size,0.0);
	
	for(unsigned int i(0); i < size; ++i)
	{
		for(unsigned int j(0); j < size; ++j)
		{
			result[i] += matrix[i][j] * vector[j];
		}
	}
	
	return result;
}

//Generate two random numbers with a Normal distribution
void rand_normal(double & r1, double & r2)
{	
	//Generate two numbers with a uniform distribution ]0,1]
	double x1(random_number.Rndm());
	double x2(random_number.Rndm());
	
	//Use the Box-Muller transformation to get two variables with a Normal distribution
	r1 = sqrt(-2.*log(x1))*cos(2.*M_PI*x2);
	r2 = sqrt(-2.*log(x1))*sin(2.*M_PI*x2);
}

//Generate a random numbers with a parabolic distribution 3/4*(-x^2 + 1)
void rand_parabolic(double & r1, double & r2)
{
	//Von Neumann - to do !
	//Parabolic function
	TF1 parabola("parabola","(3/4.)*(-pow(x,2)+1.)",-1.,1.);
	
	//Generate a number with the parabolic distribution
	double r(parabola.GetRandom());
	
	//Generate a uniform angle between 0 and 2pi
	double theta(2.*M_PI*random_number.Rndm());
	
	//Generate the random numbers in polar coordinates
	r1 = sqrt(abs(r))*cos(theta);
	r2 = sqrt(abs(r))*sin(theta);
}

//Generate two random numbers with a double Gaussian distribution
void rand_double_gauss(double & r1, double & r2, double average, double sigma, double distance) //distance is the distance between the peaks, average is the position of the hole and sigma the sigma of the two Gaussians
{
	//Generate two numbers with a uniform distribution ]0,1]
	double x1(random_number.Rndm());
	double x2(random_number.Rndm());
	double x3(random_number.Rndm());

	//Use the Box-Muller transformation to get two variables with a Normal distribution
	r1 = sqrt(-2.*log(x1) + x3*distance)*cos(2.*M_PI*x2)*sigma + average;
	r2 = sqrt(-2.*log(x1) + x3*distance)*sin(2.*M_PI*x2)*sigma + average;
}

//Compute and plot a FFT from an array of real numbers
TGraph* ComputeFFT(TCanvas *c1, int color_plot, int size_results, double *input_array, double tune)
{
	//Do the FFT
	int size_fft(size_results); //Size of the array
	TVirtualFFT *fft = TVirtualFFT::FFT(1, &size_fft, "R2C"); //Real input
	fft->SetPoints(input_array); //Input
	fft->Transform(); //Do the FFT
	double re(0.); //Store the tempory real part
	double im(0.); //Store the tempory imaginary part
	double *fft_output = new double[size_results]; //Norm of the output sqrt(re*re + im*im)
	double *frequency = new double[size_results]; //Array of the FFT frequencies
	
	//Recover the results
	for(int i(0); i < size_results/2; ++i)
	{
		fft->GetPointComplex(i, re, im); //Get the real and complex part
		fft_output[i] = sqrt(re*re + im*im); //Compute the norm
		frequency[i]=i/((double)size_results); //Correspondant frequency
	}
	
	//Plot the results
	TGraph *graph_fft = new TGraph(size_results/2, frequency, fft_output);
	graph_fft->GetXaxis()->SetTitle("Q");
	graph_fft->GetYaxis()->SetTitle("[]");
	graph_fft->SetLineColor(color_plot);
	graph_fft->GetXaxis()->SetLimits(tune - 0.02,tune + 0.01);
	//graph_fft->GetHistogram()->SetMinimum(-8.0);
	//graph_fft->GetHistogram()->SetMaximum(8.0);
	c1->cd();
	graph_fft->Draw("AL");
	c1->Update();
	
	//Clean the memory
	delete fft_output;
	delete frequency;
	
	//Return the graph
	return graph_fft;
	
}

//Find the two peaks of a FFT with their error
void FFT_find_peaks(TCanvas *c1, TGraph *graph_fft, double first_hypothesis, double second_hypothesis, double & x_1, double & FWHM_1, double & x_2, double & FWHM_2)
{
	//Clean data
	x_1 = 0.;
	x_2 = 0.;
	
	//Replot
	c1->cd();
	graph_fft->Draw("AL");
	
	//Find the real peaks and the FWHM by fitting the results with a Gaussian
	//The first peak
	TF1 *gauss_1 = (TF1*)gROOT->GetFunction("gaus"); //Create the function
	graph_fft->Fit(gauss_1, "QN", "", first_hypothesis-0.0005, first_hypothesis+0.0005); //Fit
	x_1 = gauss_1->GetParameter(1); //Get the results
	FWHM_1 = gauss_1->GetParameter(2);
	gauss_1->SetLineColor(kYellow); //Plot
	gauss_1->DrawCopy("same");
	
	//The second peak
	TF1 *gauss_2 = (TF1*)gROOT->GetFunction("gaus"); //Create the function
	graph_fft->Fit(gauss_2, "QN", "", second_hypothesis-0.0005, second_hypothesis+0.0005); //Fit
	x_2 = gauss_2->GetParameter(1); //Get the results
	FWHM_2 = gauss_2->GetParameter(2);
	gauss_2->SetLineColor(kGreen);//Plot
	gauss_2->DrawCopy("same");
	
	c1->Update();
		
}