#include <iostream>
#include <cstdlib>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <math.h>

#include "TGraph.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH1S.h"
#include "TH1I.h"
#include "TRandom3.h"
#include "Math/PdfFuncMathCore.h"
#include "TF1.h"


double Binomial(double k, double mu, float* pointer, int i) {	
	double n = gRandom->Uniform(5);
	pointer[i] = n;
	double a = (tgamma(n+k)/(tgamma(n+1)*tgamma(k)))*(pow(mu/k, n)/pow(mu/(k+1), n+k));
return a;
}

double Prob_IP() {
	double b = (gRandom->Rndm()) * 20.0;
	if (gRandom->Rndm() <= b/20.0) {
	return b;
	}
	else {
	Prob_IP();
	}
}

int Probability(double R) {
	double e = 2.71828182846;
	if (gRandom->Rndm() <= (0.1693)/(pow(e, (R-6.38)/0.535) + 1)) {
		return 1;
	}
	else {
	return 0;
	}
}

double Random() {
	double a = ((gRandom->Rndm())*14) - 7;
	return a;
}

double Round(double IP) {
	double c = IP * 5;
	int d = (int)(c + 0.5);
	c = d/(5.0);
return c;
}

void GMx106() {
	TCanvas *c1 = new TCanvas("c1", "Histogram Test", 1280, 720);
	TH1* h1 = new TH1D("h1", "Impact Parameter",100,0,20);
	TH1* h2 = new TH1I("h2", "Number of Collisions",100,0,1600);
	TH1* h3 = new TH1I("h3", "Number of Participants",100,0,400);
	
	c1->Divide(3,1);

	TCanvas *c2 = new TCanvas("c2", "Energy Production", 600, 800);
	TH1 *h4 = new TH1D("h4", "Energy Production", 100, 0, 100);

	gRandom->SetSeed(0);

	TF1 * NBD = new TF1("NBD", "[0]*ROOT::Math::negative_binomial_pdf([1],[2], x)", 0,5);

	NBD->SetParameter(0, 1.0/20.0);
	NBD->SetParameter(1, 2);
	NBD->SetParameter(2, 0.57);

	// "At separation <B>, # of Participants, # of Collisions"
	int Num_P = 0;
	int Num_temp_P = 0;
	int Num_C = 0;
	double B = 6.0;
	int Num_N1C = 0;
	int Num_N2C = 0;

	
	//Constants
	double pi = 3.14159265358979;
	double mu;
	double K;
	const	int numEvents = 1000;

	// Arrays for Nucleus 1
	int N1 [197] = {0};
	double X1 [197] = {0};
	double Y1 [197] = {0};
	double Z1 [197] = {0};
	
	// Arrays for Nucleus 2
	int N2 [197] = {0};
	double X2 [197] = {0};
	double Y2 [197] = {0};
	double Z2 [197] = {0};
	
	// Temp Variables for calculation
	double x = 0.0;
	double y = 0.0;
	double z = 0.0;
	double Rmag = 0.0;
	double d2 = 0.0;

	// More arrays for 10^6
	float Bx [numEvents] = {0};
	//Bx = malloc(sizeof(float) * numEvents);
	//int By [101] = {0};
	short Ncollx [numEvents] = {0};
	//Ncoll = malloc(sizeof(short) * numEvents);
	//int Ncolly [2000] = {0};
	short Npartx [numEvents] = {0};
	//Npartx = malloc(sizeof(short) * numEvents);
	//int Nparty [400] = {0};
	float NBDx [numEvents] = {0};
	//NBDx = malloc(sizeof(short) * numEvents);

	//srand( time( NULL ));

for (int l = 0; l<numEvents; l++) {
	B = Prob_IP();
	Num_temp_P = 0;
	Num_N1C = 0;
	Num_N2C = 0;
	memset(N1, 0, sizeof(N1));
	memset(N2, 0, sizeof(N2));

std::cout << l << std::endl;
	//Initializes Nucleus #1 and stores in X1 Y1 and Z1 Arrays
	for (int counter = 0; counter < 197; ) {
		x = Random();
		y = Random();
		z = Random();
		Rmag = sqrt(x*x + y*y + z*z);
		if (Rmag <= 7.0 && Probability(Rmag) == 1) {
			X1[counter] = x;
			Y1[counter] = y;
			Z1[counter] = z;
			counter += 1;
		}
	}

	//Initializes Nucleus #2 and stores in X2 Y2 and Z2 Arrays
	for (int counter2 = 0; counter2 < 197; ) {
		x = Random();
		y = Random();
		z = Random() - B;
		Rmag = sqrt(x*x + y*y + (z + B)*(z + B));
		if (Rmag <= 7.0 && Probability(Rmag) == 1) {
			X2[counter2] = x;
			Y2[counter2] = y;
			Z2[counter2] = z;
			counter2 += 1;
		}
	}
	
	//Compares each pair of Nucleons to see if a collision occurred, and to count the participants
	for (int i = 0; i < 197; i++) {
		for (int j = 0; j < 197; j++) {
			d2 = (pow(Y1[i] - Y2[j], 2.0) + pow(Z1[i] - Z2[j], 2.0));
			if ( d2 < 4.0/pi ) {
				N1[i] += 1;
				N2[j] += 1;
			}
		}
	}

	//For loops to count number of participants and number of collisions
	for (int k = 0; k < 197; k++) {
		Num_N1C += N1[k];
		Num_N2C += N2[k];
		if (N1[k] > 0) {
			Num_P += 1;
			Num_temp_P += 1;
		}
		if (N2[k] > 0) {
			Num_P += 1;
			Num_temp_P += 1;
		}
	}
	if (Num_N1C != Num_N2C) {
		std::cout << "Num_N1C = " << Num_N1C <<" and Num_N2C = " << Num_N2C << std::endl;
continue;
	}
	for (int h = 0; h<Num_N1C; ++h) {
		NBDx[l] += NBD->GetRandom();
	}





		Num_C += Num_N1C;
	
		Ncollx[l] = Num_N1C;
		//Ncolly[Num_N1C] += 1;
		Npartx[l] = Num_temp_P;
		//Nparty[Num_temp_P] += 1;
		Bx[l] = (float)B;
		//By[(int)(Bx[l]*5+0.5)] += 1;
}
//for (int n = 0; n < 1600; ++n) {
//	h4->Fill(NBDy[n]);
//}
	for (int m = 0; m < numEvents; ++m) {
		if (Ncollx[m] > 0 && Npartx[m] > 0) {	
			h1->Fill(Bx[m]);
			h2->Fill(Ncollx[m]);
			h3->Fill(Npartx[m]);
			h4->Fill(NBDx[m]);
		}
}
	c1->cd(1);
	h1->Draw();
	c1->cd(2);
	gPad->SetLogy(1);
	h2->Draw();
	c1->cd(3);
	gPad->SetLogy(1);
	h3->Draw();
//	c1->SaveAs("/home/jmerges/public_html/tests/Final3Histos.png");

	c2->cd();
	gPad->SetLogy(1);
	h4->Draw();
//	c2->SaveAs("/home/jmerges/public_html/tests/FinalEnergy.png");
std::cout << "Done!" << std::endl;

return;
}