void Draft()
{
	
	TFile *input = new TFile("AuAu200GeV_SM_NT3_DecaysOff_0_20fm_file0.root", "read");
	
	TTree *tree = (TTree*)input->Get("tr");
	
	Int_t Event, Mult;
	Float_t Imp; //Impact parameter
	
	tree->SetBranchAddress("Event", &Event);
	tree->SetBranchAddress("Mult", &Mult);
	tree->SetBranchAddress("Imp", &Imp);
		 
	Float_t Array[10000]; //Impact parameter array
	Int_t MArray[10000];	// Multiplicity array
	
	
	/* For multiplicities*/
	
	Int_t PID[12000]; //Particle ID	
	Float_t Mass[12000]; //Mass
	Float_t E[12000]; //Energy
	Float_t y[12000]; //Rapidity
	Float_t Pt[12000]; //Transverse Momentum 
	Float_t P[120000]; // |P| mod P
	Float_t eta[12000]; //Pseudo-rapidity
	Float_t phi[12000]; //Azimuthal angle
	
	tree->SetBranchAddress("PID",PID);
	tree->SetBranchAddress("Mass",Mass);
	
	Float_t Px[12000];
	Float_t Py[12000];
	Float_t Pz[12000];
	
	tree->SetBranchAddress("Px",Px);
	tree->SetBranchAddress("Py",Py);
	tree->SetBranchAddress("Pz",Pz);
	
	
	Int_t Stat; //Stat = 0 means Spectators Stat>0 means participants
	
	Float_t Nx;
	Float_t Ny;
	Float_t Nz;
	
	Float_t Nx_arr[394] = {};
	Float_t Ny_arr[394] = {};
	Float_t Nz_arr[394] = {};
	Int_t Stat_arr[394] = {};
	
	tree->SetBranchAddress("Nx", &Nx_arr);
	tree->SetBranchAddress("Ny", &Ny_arr);
	tree->SetBranchAddress("Nz", &Nz_arr);
	tree->SetBranchAddress("Stat", &Stat_arr);
	
	
	
	Int_t entries = tree->GetEntries();
		
	
	/* Rapidity (y), Pseudorapidity (eta), Transverse momentum (Pt), Azimuthal angle (phi) */ 
	
	TCanvas *c8 = new TCanvas();
	TH1F *hist1 = new TH1F("hist1", "Pt for S = 0", 100, 0, 4);
		
		
	for(Int_t i = 0; i<entries; i++)
	{	
		tree->GetEntry(i);
		for(Int_t j = 0; j<Mult; j++)
		{	
			E[j] = TMath::Sqrt((TMath::Power(Px[j],2)) + (TMath::Power(Py[j],2)) + (TMath::Power(Pz[j],2)) + (TMath::Power(Mass[j],2)));	// Energy
			y[j] = 0.5 * TMath::Log((E[j] + Pz[j]) / (E[j] - Pz[j]));	// Rapidity
			Pt[j] = TMath::Sqrt(TMath::Power(Px[j],2) + TMath::Power(Py[j],2));	// Transverse Momentum
			P[j] = TMath::Sqrt(TMath::Power(Px[j],2) + TMath::Power(Py[j],2) + TMath::Power(Pz[j],2));	// |P| Magnitude of Momentum
			eta[j] = 0.5 * TMath::Log((P[j] + Pz[j]) / (P[j] - Pz[j]));	// Pseudo-rapidity
			phi[j] = TMath::ATan2(Py[j],Px[j]);	// Azimuthal angle
			hist1->Fill(Pt[j]);
		}
			 
	}
	c8->cd();
	hist1->Draw();
	
	
	/* Getting Collision Centrality from Impact Parameter Histogram */
	
	TCanvas *c9 = new TCanvas();
	
	TH1F *myhist = new TH1F("myhist", "Impact Parameter b", 100, 0., 0.); // auto binning
	for (Int_t i=0; i<entries; i++){
		tree->GetEntry(i);
		myhist->Fill(Imp);
	}
	
	c9->cd();
	myhist->Draw();
	
	Int_t TBins = myhist->GetNbinsX(); // Total bins = 100 	
	
	Int_t sum;
	Int_t binnumber;
	Float_t bc = 0; // bin center value (Impact parameter) of a bin
	
	/*Int_t sum_1;
	Int_t binnumber_1;
	Float_t bc_1 = 0; // bin center value (Impact parameter) of a bin
	*/
	Double_t binentries[12000];	
	
	Int_t A[5] = {100, 200, 300, 400, 500};
	
	for(Int_t j = 0;j<5;j++){
		sum = 0;
		for(Int_t i = 0; i<TBins; i++){
		binentries[i] = 0;
		binnumber = 0;
		binentries[i] = myhist->GetBinContent(i);
		sum = sum + binentries[i];
		if(sum>A[j]){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			//cout << binnumber << "  " << bc << "  " << endl;
			break;
		}
	}}
	
		
	Float_t r_2[394] = {}; // r^2 = Nx^2 +Ny^2 
	Float_t r_2_sum[10000]; // r_2[0] + r_2[1] + .... + r_2[Mult].	It is for each collision. Index runs over i. r_2_sum[i]
	Float_t r_2_avg[10000]; // <r^2> .	It is for each collision. Index runs over i. r_2_avg[i]
		
	
	Float_t e2_cos[394] = {}; // cos term for each participant for n = 2 
	Float_t e2_sin[394] = {}; // sin term for each participant for n = 2
	
	Float_t e2_cos_sum[10000]; // sum of cos terms for each participant for n = 2
	Float_t e2_sin_sum[10000]; // sum of sin terms for each participant for n = 2
	
	Float_t e2_cos_sum_avg[10000]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision
	Float_t e2_sin_sum_avg[10000]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision
	
	Float_t e2_num[10000]; // Numerator of eccentricity expression for n = 2
	Float_t e2[10000]; //Eccentricity for each event for n = 2
	Float_t e2_2[10000]; // Eccentricity^2  for each event for n = 2
	
	
	
	Float_t e3_cos[394] = {}; // cos term for each participant for n = 3 
	Float_t e3_sin[394] = {}; // sin term for each participant for n = 3
	
	Float_t e3_cos_sum[10000]; // sum of cos terms for each participant for n = 3 
	Float_t e3_sin_sum[10000]; // sum of sin terms for each participant for n = 3
	
	Float_t e3_cos_sum_avg[10000]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision
	Float_t e3_sin_sum_avg[10000]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision
	
	
	Float_t e3_num[10000]; // Numerator of eccentricity expression for n = 3
	Float_t e3[10000]; //Eccentricity for each event for n = 3
	Float_t e3_2[10000]; // Eccentricity^2  for each event for n = 3
	
	Float_t phi_part[394] = {}; //Azimuthal angle of participating nucleons
	Int_t N[10000];	// For counting the number of participants.
	Int_t cent = 0;
	
	Float_t sum_e2_2[5] = {};
	Float_t sum_e3_2[5] = {};
	Int_t counter[5] = {};
	Float_t RMS_2[5] = {};
	Float_t RMS_3[5] = {};
	Float_t centr[5] = {1, 2, 3, 4, 5};		
	
	for(Int_t i = 0; i<entries; i++)
		{	
			tree->GetEntry(i);
			Array[i] = Imp;
			MArray[i] = Mult;
			N[i] = 0;
			e2_2[i] = 0;
			e3_2[i] = 0;
			r_2_sum[i] = 0;
			e2_cos_sum[i] = 0;
			e2_sin_sum[i] = 0;
			e3_cos_sum[i] = 0;
			e3_sin_sum[i] = 0;
			if(Array[i]<=1.862)cent = 1;
			else if(Array[i]>1.862 && Array[i]<=2.45)cent = 2;
			else if(Array[i]>2.45 && Array[i]<=3.038)cent = 3;
			else if(Array[i]>3.038 && Array[i]<=3.234)cent = 4;
			else if(Array[i]>3.234 && Array[i]<=3.626)cent = 5;
			for(Int_t j = 0; j<394; j++){
				
				//cout << Nx_arr[i] << "  " << Ny_arr[i] << "  " << Nz_arr[i] << "  " << Stat_arr[i] << endl;
				
				if(Stat_arr[j] > 0){
				
				
					
					phi_part[j] = TMath::ATan2(Ny_arr[j],Nx_arr[j]);
					
					N[i]+=1; //counts the number of participating nucleons for each collision
					
					r_2[j] = ((TMath::Power(Nx_arr[j],2)) + (TMath::Power(Ny_arr[j],2)));
					
					r_2_sum[i] = r_2_sum[i] + r_2[j];
									
					e2_cos[j]= r_2[j]*(TMath::Cos((2*phi_part[j])));
					e2_sin[j]= r_2[j]*(TMath::Sin((2*phi_part[j])));
				
					e2_cos_sum[i] = e2_cos_sum[i] + e2_cos[j];
					e2_sin_sum[i] = e2_sin_sum[i] + e2_sin[j];
					
					e3_cos[j]= r_2[j]*(TMath::Cos((3*phi_part[j])));
					e3_sin[j]= r_2[j]*(TMath::Sin((3*phi_part[j])));
					
					e3_cos_sum[i] = e3_cos_sum[i] + e3_cos[j];
					e3_sin_sum[i] = e3_sin_sum[i] + e3_sin[j];
					 	
				}
		}
		
		r_2_avg[i] = r_2_sum[i]/N[i];
		e2_cos_sum_avg[i] = e2_cos_sum[i]/N[i];
		e2_sin_sum_avg[i] = e2_sin_sum[i]/N[i];
		
 		
 		e2_num[i] = TMath::Sqrt((TMath::Power(e2_cos_sum_avg[i],2)) + (TMath::Power(e2_sin_sum_avg[i],2))); //Numerator term for eccentricity 
		e2[i] = e2_num[i]/r_2_avg[i];
		e2_2[i] = TMath::Power(e2[i],2);
		
		e3_cos_sum_avg[i] = e3_cos_sum[i]/N[i];
		e3_sin_sum_avg[i] = e3_sin_sum[i]/N[i];
				
		
		e3_num[i] = TMath::Sqrt((TMath::Power(e3_cos_sum_avg[i],2)) + (TMath::Power(e3_sin_sum_avg[i],2))); //Numerator term for eccentricity 
		e3[i] = e3_num[i]/r_2_avg[i];
		e3_2[i] = TMath::Power(e3[i],2);
		
		
				
		if(Array[i]<=1.862){
			sum_e2_2[0] = sum_e2_2[0] + e2_2[i];	// 0 - 1% 
			sum_e3_2[0] = sum_e3_2[0] + e3_2[i];
			counter[0]+=1;
			
		}
		else if(Array[i]>1.862 && Array[i]<=2.45){
			sum_e2_2[1] = sum_e2_2[1] + e2_2[i];	// 1 - 2%
			sum_e3_2[1] = sum_e3_2[1] + e3_2[i];
			counter[1]+=1;
		}
		else if(Array[i]>2.45 && Array[i]<=3.038){
			sum_e2_2[2] = sum_e2_2[2] + e2_2[i];	// 2 - 3%
			sum_e3_2[2] = sum_e3_2[2] + e3_2[i];
			counter[2]+=1;
		}
		else if(Array[i]>3.038 && Array[i]<=3.234){
			sum_e2_2[3] = sum_e2_2[3] + e2_2[i];	// 3 - 4%
			sum_e3_2[3] = sum_e3_2[3] + e3_2[i];
			counter[3]+=1;
		}
		else if(Array[i]>3.234 && Array[i]<=3.626){
			sum_e2_2[4] = sum_e2_2[4] + e2_2[i];	// 4 - 5%
			sum_e3_2[4] = sum_e3_2[4] + e3_2[i];
			counter[4]+=1;
		}
		
	}
	
  	
	for(Int_t l = 0; l<5; l++){
		RMS_2[l] = TMath::Sqrt((sum_e2_2[l]/counter[l]));
		RMS_3[l] = TMath::Sqrt((sum_e3_2[l]/counter[l]));
	}
	
	TCanvas *c16 = new TCanvas();
	c16->cd();
	
	TGraph* gr1 = new TGraph(5, centr, RMS_2);
	gr1->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 2 S = 0;Centrality %(Imp);RMS Eccentricity n = 2");
	gr1->SetLineColor(2);
	gr1->SetMarkerColor(4);
	gr1->SetMarkerStyle(21);
	gr1->Draw("ACP");
	
	TCanvas *c17 = new TCanvas();
	c17->cd();
	
	TGraph* gr2 = new TGraph(5, centr, RMS_3);
	gr2->SetTitle("RMS Eccentricity vs Centrality % (Imp b) n = 3 S = 0;Centrality %(Imp);RMS Eccentricity n = 3");
	gr2->SetLineColor(3);
	gr2->SetMarkerColor(5);
	gr2->SetMarkerStyle(22);
	gr2->Draw("ACP");
	
		
	/*TCanvas *c22 = new TCanvas();
	auto MG1 = new TMultiGraph();
	MG1->SetTitle("RMS Eccentricity vs Centrality % (Imp b) for S = 0;Centrality % (Imp); RMS Eccentricity");
	MG1->Add(gr1, "cp");
	MG1->Add(gr2, "cp");
	MG1->Draw("acp");
	c22->BuildLegend();*/
	
			
	/*TFile* f1 = new TFile("Draft0.root", "CREATE");
	gr1->Write("Graph1");
	gr2->Write("Graph2");
	MG1->Write("GS0");
	hist1->Write();*/
	
	//input->Close();
}