#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"
#include "TFile.h"
#include "TLegend.h"
#include "TH2F.h"
#include "TTree.h"
#include "TStyle.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;
}

int main() {
	TTree * PrScn = new TTree("PrScn", "Parameter Scan");



	TCanvas *c1 = new TCanvas("c1", "Histogram Test", 1280, 720);
	
	c1->Divide(3,1);

	TCanvas *c2 = new TCanvas("c2", "Energy Production", 600, 800);


	gRandom->SetSeed(0);

	TFile * data = TFile::Open("10-001_Sum_ET_HF_Energy.root");

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

	NBD->SetParameter(0, 1);
	NBD->SetParameter(1, 2);
	NBD->SetParameter(2, 0.6);


	// "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 = 100000;

	// 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;
	double databin = 0.0;
	double MCbin = 0.0;
	double Chi2 = 0.0;
	double Ch_g = 0.0;
	double error = 0.0;

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

	// Simple Structs for TTree
	//typedef struct {float k,p,chd;} PARSCAN;

	/*static PARSCAN ps1;
	static PARSCAN ps2;
	static PARSCAN ps3;
	static PARSCAN ps4;
	static PARSCAN ps5;
	static PARSCAN ps6;
	static PARSCAN ps7;
	static PARSCAN ps8;
	static PARSCAN ps9;*/

	const int NUM_TESTS = 4;

	float psK[NUM_TESTS*NUM_TESTS];
	float psP[NUM_TESTS*NUM_TESTS];

	PrScn->Branch("K", psK,"psK[NUM_TESTS*NUM_TESTS]/F");
	PrScn->Branch("P", psP, "psP[NUM_TESTS*NUM_TESTS]/F");
//	PrScn->Branch("CHI", psCHI, "psCHI[NUM_TESTS*NUM_TESTS]/F");
//	PrScn->Branch("CH_G", psCH_G, "psCH_G[NUM_TESTS*NUM_TESTS]/F");

TH2 * hps = new TH2D("hps", "Parameter Scan",4,0.51,0.59,4,0.72,0.88);	

int ia = 0;//for (int ia = 0;ia<NUM_TESTS;++ia) {
int ib = 0;//for (int ib = 0;ib<NUM_TESTS;++ib) {

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);

TH1 *h4 = new TH1F("h4", "Monte Carlo Glauber", 100, 0, 5);
TH1 * hdata = 0;
data->GetObject("hSumEnergy_nominal", hdata);


memset(Bx, 0, sizeof(Bx));
memset(Ncollx, 0, sizeof(Ncollx));
memset(Npartx, 0, sizeof(Npartx));
memset(NBDx, 0, sizeof(NBDx));


NBD->SetParameter(1, 0.75+(ia/25.0));
NBD->SetParameter(2, 0.55+(ib/50.0));
int l = 0;
for (int n = 0; n<numEvents; n++) {
	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 << "ia: " << ia << "\tib: " << ib << "\tn: " << n << std::endl;
//if (Ch_g != 0) {
	std::cout << "Chi2: " << Chi2 << "\tGetBinContents: " << h4->GetBinContent(50) << "\t" << h4->GetBinContent(40) << "\t" << h4->GetBinContent(30) << 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()) * 1.0/438.0;
	}





		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;
	if (Num_N1C>0) ++l;
}

//for (int n = 0; n < 1600; ++n) {
//	h4->Fill(NBDy[n]);
//}
	for (int m = 0; m < (numEvents/3)*2; ++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->SetLineColor(kRed);
	h4->Draw("");
	hdata->Scale(h4->GetEntries());
	hdata->Draw("same");

	TLegend * legend = new TLegend(0.2,0.6,0.70,0.90);
	legend->AddEntry(h4, "Monte Carlo Glauber", "l");
	legend->AddEntry(hdata, "Sum_ET_HF_Energy.root", "l");
	legend->Draw();

	c2->SaveAs("/home/jmerges/public_html/par_scan/ChiFit_2.png");
//Chi Square:
Chi2 = 0.0;
	for (int i=1;i<=100;++i) {
		databin = hdata->GetBinContent(i);
		MCbin = h4->GetBinContent(i);
		error = hdata->GetBinError(i);
		if (error>0) Chi2 += (pow(databin-MCbin,2)/pow(error,2));
	}
//Chi2 = hdata->Chi2Test(h4, "CHI2");
//Ch_g = hdata->Chi2Test(h4, "WW CHI2/NDF");
//psK[(NUM_TESTS*ia)+ib] = NBD->GetParameter(1);
//psP[(NUM_TESTS*ia)+ib] = NBD->GetParameter(2);
//psCHI[(NUM_TESTS*ia)+ib] = Chi2;
//psCH_G[(NUM_TESTS*ia)+ib] = Ch_g;
hps->SetBinContent(hps->GetBin(ib+1,ia+1),Chi2);

printf("Chi2: %lf\n", Chi2);

//destroy histograms
h1->~TH1();
h2->~TH1();
h3->~TH1();
h4->~TH1();
hdata->~TH1();
//}
//}//end parameter scan for loops
PrScn->Fill();
//PrScn->SaveAs("/home/jmerges/public_html/tests/PrScn2.root");
TCanvas * c3 = new TCanvas("c3", "Parameter Scan", 1200,900);
c3->cd();
//gPad->SetLogz(1);
hps->SetStats(0);
hps->Draw("COLZ");
c3->SaveAs("/home/jmerges/public_html/par_scan/K_074_086_P_052_58.png");
std::cout << "Chi2: " << Chi2 << std::endl;
std::cout << "Done!" << std::endl;

return 0;
}