#include "TInterpreter.h"
#include "TRandom3.h"
#include "TDatime.h"
#include "TFile.h"
#include "TH1F.h"
#include "TMath.h"
#include "TF1.h"
#include "TCanvas.h"
void particleProduction(){
	
	// This macro simulates nucleus nucleus collisions and tracks the number of participants, 
	// collisions, and impact parameter for each collision	
	float p0 = 1; 		//nucleon density at center of nucleus 	[0]
	float w = 0; 		//deviations from spherical shape		[1]	
	float R = 6.62; 	//nuclear radius (fm)					[2]
	float a = .546; 	//skin depth (fm)						[3]
	const float A=208;  //mass of Au
	float r = 0; 		//radius
	float b; 			//impact parameter (fm) ((6,5,4) is a temporary value. Value of b is assigned later)
	int n_inc = 59448;	//number of incident collisions
	int n_col = 0;		//number of nucleon-nucleon collisions for each incident nucleus
	int n_part = 0;		//number of nucleon participants in each Au-Au collision	
	float mu = 1.6;		//binomial coefficient mu				[0]
	float k_param=1.00;	//binomial coefficient k 				[1]
	float n = 100; 		// binomial coefficient n 				
	float E_t = 0;
	
	//declaration of initial histos
	TH1F *bi_col_histo = new TH1F("E_{T}", "E_{T}",50, 0, 2500);


	//define negative binomial distribution

	TF1 *bi_distru = new TF1("Negative Binomial Distribution","TMath::Binomial(x+[1]-1,x)*TMath::Power([0]/[1],x)/(TMath::Power([0]/[1]+1,x+[1]))",0,150); //for Au
	bi_distru->SetParameters(mu,k_param);
	TCanvas* c1 = new TCanvas("c1","Practice1",800,600);
	//bi_distru->Draw();


	//define nuclear density function
	TF1 *nuc_density = new TF1("Nuclear Density Function","[0]*(1+[1]*TMath::Power(x/[2],2))/(1+TMath::Exp((x-[2])/[3]))",0,10); //for Au
	nuc_density->SetParameters(p0,w,R,a);
	nuc_density->GetHistogram()->SetTitle("Nuclear Density");
	//c1->SaveAs("/home/christos/public_html/protected/c1.png");

	//b parameter for hard nucleus
	TH1F *bsphere_histo = new TH1F("Impact Parameter Assuming Hard Nucleus", "# of Collisions Vs. Impact Parameter",50, 0, 15.5);
	float b_param = 0.0; //b in units of fm
	float db = 1.0; //db in units of fm
	float lumi = 15.0; //lumi in units of /fb

	//Random # generator
	TDatime *my_clock = new TDatime();
	TRandom3 rgen(my_clock->GetTime());
	
	for(int i = 0; i<10000; i++){
		b_param = rgen.Uniform(0,15);
		if (1) {
			for(int j = 0; j<b_param*lumi*2*TMath::Pi()*db; j++){
				bsphere_histo->Fill(b_param);
			}
		}
	}
	//c0->SaveAs("/home/christos/public_html/protected/c0.png");

	//arrays to store points of all nucleons (used in part 3)
	float nuc1x[A];
	float nuc1y[A];
	float nuc2x[A];
	float nuc2y[A];

	float cosTheta, phi, theta;
	TF1 *nuc_prob = new TF1("Nuclear Probability Function","[0]*(1+[1]*TMath::Power(x/[2],2))/(1+TMath::Exp((x-[2])/[3]))*TMath::Power(x,2)",0,10); //for Au
	nuc_prob->SetParameters(p0,w,R,a);
	std::cout << "--- Processing: " << n_inc << " collisions" << std::endl;
	TStopwatch sw;
   	sw.Start();

	for(int k = 0; k<n_inc; k++){
		if(k%100 == 0){
			cout<<"--- Processing collision: "<<k<<endl;
		}

		//Assign random value to B
		/*b(0) = bsphere_histo->GetRandom();
		b(1) = bsphere_histo->GetRandom();
		b(2) = bsphere_histo->GetRandom();*/

		b = bsphere_histo->GetRandom();
		

		//Assign Spatial coordinates to nucleus 1
		for (int i = 0; i <A; i++){
			r = nuc_prob->GetRandom();
			cosTheta = rgen.Uniform(-1,1);
			phi = rgen.Uniform(0,2*TMath::Pi());
			theta = TMath::ACos(cosTheta);
			nuc1x[i] = (r*TMath::Sin(theta)*TMath::Cos(phi));
			nuc1y[i] = (r*TMath::Sin(theta)*TMath::Sin(phi));
		}
		
		//Assign spacial coordinates to nucleus 2
		for (int i = 0; i < A; i++){
			r = nuc_prob->GetRandom();
			cosTheta = rgen.Uniform(-1,1);
			phi = rgen.Uniform(0,2*TMath::Pi());
			theta = TMath::ACos(cosTheta);
			nuc2x[i] = (b+r*TMath::Sin(theta)*TMath::Cos(phi));
			nuc2y[i] = (r*TMath::Sin(theta)*TMath::Sin(phi));
		}
		
		//Fill graphs, count collision, count partipants
		for(int i = 0; i<A; i++){
			for(int j = 0; j<A; j++){
				if(TMath::Power((nuc1x[i]-nuc2x[j]),2)+TMath::Power((nuc1y[i]-nuc2y[j]),2) <2.23){ //assuming 70 mb xsec
					//Fill impact parameter histogram
					E_t += bi_distru->GetRandom();
				}
			}
		}


    	bi_col_histo->Fill(E_t);
    	E_t = 0;
		n_col = 0;
	}
	sw.Stop();
    std::cout << "--- End of collision loop: "; sw.Print();
	TCanvas* c2 = new TCanvas("c2","Practice2",800,600);
	//c2->SetLogy();
	bi_col_histo->Draw();
	TFile *outfile = new TFile("/root_v5.34.32/test/10-004_Sum_ET_HF_Energy.root");
}