void A3()
{
	
	TFile *input = new TFile("AuAu200GeV_SM_NT3_DecaysOff_0_20fm_file0.root", "read");
	
	TTree *tree = (TTree*)input->Get("tr");
	
	Int_t Event, Mult;
	
	Int_t Stat[12000]; //Stat = 0 means Spectators Stat>0 means participants
	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
 	
	Float_t Imp; //Impact parameter 
	Float_t Array[12000]; //Impact parameter array
	
	Int_t MArray[12000];	// Multiplicity array
	
	Float_t phi_part[12000]; //Azimuthal angle of participating nucleons
	
	Float_t r_2[12000]; // r^2 = Nx^2 +Ny^2 +Nz^2 It is for each multiplicity. Index runs over j. r_2[j]
	Float_t r_2_sum[12000]; // 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[12000]; // <r^2> .	It is for each collision. Index runs over i. r_2_avg[i]
	
	
	
	Float_t e2_cos[12000]; // cos term for each participant for n = 2 
	Float_t e2_sin[12000]; // sin term for each participant for n = 2
	
	Float_t e2_cos_sum[12000]; // sum of cos terms for each participant for n = 2
	Float_t e2_sin_sum[12000]; // sum of sin terms for each participant for n = 2
	
	Float_t e2_cos_sum_avg[12000]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision
	Float_t e2_sin_sum_avg[12000]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision
	
	Float_t e2_num[12000]; // Numerator of eccentricity expression for n = 2
	Float_t e2[12000]; //Eccentricity for each event for n = 2
	Float_t e2_2[12000]; // Eccentricity^2  for each event for n = 2
	
	
	
	Float_t e3_cos[12000]; // cos term for each participant for n = 3 
	Float_t e3_sin[12000]; // sin term for each participant for n = 3
	
	Float_t e3_cos_sum[12000]; // sum of cos terms for each participant for n = 3 
	Float_t e3_sin_sum[12000]; // sum of sin terms for each participant for n = 3
	
	Float_t e3_cos_sum_avg[12000]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision
	Float_t e3_sin_sum_avg[12000]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision
	
	Float_t e3_num[12000]; // Numerator of eccentricity expression for n = 3
	Float_t e3[12000]; //Eccentricity for each event for n = 3
	Float_t e3_2[12000]; // Eccentricity^2  for each event for n = 3
	
	
	Float_t Nx[12000];
	Float_t Ny[12000];
	Float_t Nz[12000];
	
	
	Float_t Px[12000];
	Float_t Py[12000];
	Float_t Pz[12000];
	
	tree->SetBranchAddress("Event", &Event);
	tree->SetBranchAddress("Mult", &Mult);
	tree->SetBranchAddress("Imp", &Imp);
	
	tree->SetBranchAddress("PID",PID);
	tree->SetBranchAddress("Mass",Mass);
	
	tree->SetBranchAddress("Nx",Nx);
	tree->SetBranchAddress("Ny",Ny);
	tree->SetBranchAddress("Nz",Nz);
	
	
	tree->SetBranchAddress("Px",Px);
	tree->SetBranchAddress("Py",Py);
	tree->SetBranchAddress("Pz",Pz);
	
	Int_t entries = tree->GetEntries();
	

	/* Rapidity (y), Pseudorapidity (eta), Transverse momentum (Pt), Azimuthal angle (phi) */ 
		
		
	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
		}
			 
	}
	
	/*
		
	
		
	for(Int_t i = 0; i<entries; i++)
	{	
		tree->GetEntry(i);
		K[i] = 0;
		MArray[i] = Mult;
		
		for(Int_t j = 0; j<Mult; j++){
			if(Stat[j] > 0){
				K[i]+=1;	//counts the number of participating nucleons for each collision
			}
		}
		cout << Mult << "  " << MArray[i] << "  " << K[i] << endl;
		hprof5->Fill(K[i],MArray[i]);
	}
	
	
	
	*/
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
		
	
	/* 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 
	
	//Double_t F = myhist->GetXaxis()->GetBinLowEdge(3);
	//Double_t G = myhist->GetXaxis()->GetBinCenter(35);
	
	//cout << F << endl;
	//cout << G << endl;
	
	
	Int_t sum;
	Int_t binnumber;
	Float_t bc = 0; // bin center value (Impact parameter) of a bin
	Float_t bl = 0; // bin low  value Impact parameter of a bin
	
	Double_t binentries[12000];
	Int_t A[8] = {1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000};
	
	for(Int_t j = 0; j<8; j++){
		sum = 0;
		for(Int_t i = 0; i<TBins; i++){
		binentries[i] = 0;
		binnumber = 0;
		binentries[i] = myhist->GetBinContent(i); // Gives the number of enteries in each bin
		//cout << binentries[i] << endl;
		sum = sum + binentries[i];
		if(sum>A[j]){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			//bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  " << endl;
			break;}
			
		}
		}
		
	
	TCanvas *c10 = new TCanvas();
	TCanvas *c11 = new TCanvas();
	TCanvas *c12 = new TCanvas();
	TCanvas *c13 = new TCanvas();
	TCanvas *c18 = new TCanvas();
	
	TProfile *hprof5 = new TProfile("hprof5", "npart vs Mult", 100, 0, 10000);
	
	TProfile *hprof1 = new TProfile("hprof1", "Eccentricity vs b for n = 2 for S = 0", 40, 0, 20);
	TProfile *hprof2 = new TProfile("hprof2", "Eccentricity vs Centrality for n = 2 for S = 0", 40, 0, 90);
	
	TProfile *hprof3 = new TProfile("hprof3", "Eccentricity vs b for n = 3 for S = 0", 40, 0, 20);
	TProfile *hprof4 = new TProfile("hprof4", "Eccentricity vs Centrality for n = 3 for S = 0", 40, 0, 90);
		
	Int_t N[12000];	// For counting the number of participants
	
	//Int_t n = 2;
	
	Int_t cent = 0;
	Float_t sum_e2_2[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t sum_e3_2[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Int_t counter[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t RMS_2[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t RMS_3[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t centr[9] = {10, 20, 30, 40, 50, 60, 70, 80, 90};
	
	Float_t E23avg[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};	// <e2.e3>
	Float_t E23[12000];	// e2[i]*e3[i] for each collision
	Float_t sum_E23[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	
	Float_t sum_X[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t sum_X_avg[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	
	Float_t sum_Y[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	Float_t sum_Y_avg[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	
	Float_t W[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
	for(Int_t i = 0; i<entries; i++)
		{	
			tree->GetEntry(i);
			N[i] = 0;
			e2_2[i] = 0;
			e3_2[i] = 0;
			Array[i] = Imp;
			MArray[i] = Mult;
			if(Array[i]<=4.998)cent = 10;
			else if(Array[i]>4.998 && Array[i]<=6.958)cent = 20;
			else if(Array[i]>6.958 && Array[i]<=8.526)cent = 30;
			else if(Array[i]>8.526 && Array[i]<=9.898)cent = 40;
			else if(Array[i]>9.898 && Array[i]<=11.074)cent = 50;
			else if(Array[i]>11.074 && Array[i]<=11.858)cent = 60;
			else if(Array[i]>11.858 && Array[i]<=12.838)cent = 70;
			else if(Array[i]>12.838 && Array[i]<=13.622)cent = 80;
			else if(Array[i]>13.622)cent = 90;
			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;
			for(Int_t j = 0; j<Mult; j++){
						
				if(Stat[j] > 0){
				
					phi_part[j] = TMath::ATan2(Py[j],Px[j]);
					r_2[j] = (TMath::Power(Nx[j],2)) + (TMath::Power(Ny[j],2));
					r_2_sum[i] = r_2_sum[i] + r_2[j];
					N[i]+=1; //counts the number of participating nucleons for each collision
				
					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];
					
				
				}
		}
		cout << Mult << "  " << MArray[i]<< "  " << N[i] << endl;
		hprof5->Fill(N[i],MArray[i]);
		
		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);
		E23[i] = e2[i] * e3[i];	// e2*e3 for each collision
		
		//cout << eccen[i] << "  " << eccen_2[i] << endl;
		hprof1->Fill(Array[i],e2[i]);
		hprof2->Fill(cent,e2[i]);
		hprof3->Fill(Array[i],e3[i]);
		hprof4->Fill(cent,e3[i]);
		
		if(Array[i]<=4.998){
			sum_e2_2[0] = sum_e2_2[0] + e2_2[i];	// 0 - 10% 
			sum_e3_2[0] = sum_e3_2[0] + e3_2[i];
			sum_E23[0] = sum_E23[0] + E23[i];
			sum_X[0] = sum_X[0] + e2[i];
			sum_Y[0] = sum_Y[0] + e3[i];
			counter[0]+=1;
			
		}
		else if(Array[i]>4.998 && Array[i]<=6.958){
			sum_e2_2[1] = sum_e2_2[1] + e2_2[i];	// 10 - 20%
			sum_e3_2[1] = sum_e3_2[1] + e3_2[i];
			sum_E23[1] = sum_E23[1] + E23[i];
			sum_X[1] = sum_X[1] + e2[i];
			sum_Y[1] = sum_Y[1] + e3[i];
			counter[1]+=1;
		}
		else if(Array[i]>6.958 && Array[i]<=8.526){
			sum_e2_2[2] = sum_e2_2[2] + e2_2[i];	// 20 - 30%
			sum_e3_2[2] = sum_e3_2[2] + e3_2[i];
			sum_E23[2] = sum_E23[2] + E23[i];
			sum_X[2] = sum_X[2] + e2[i];
			sum_Y[2] = sum_Y[2] + e3[i];
			counter[2]+=1;
		}
		else if(Array[i]>8.526 && Array[i]<=9.898){
			sum_e2_2[3] = sum_e2_2[3] + e2_2[i];	// 30 - 40%
			sum_e3_2[3] = sum_e3_2[3] + e3_2[i];
			sum_E23[3] = sum_E23[3] + E23[i];
			sum_X[3] = sum_X[3] + e2[i];
			sum_Y[3] = sum_Y[3] + e3[i];
			counter[3]+=1;
		}
		else if(Array[i]>9.898 && Array[i]<=11.074){
			sum_e2_2[4] = sum_e2_2[4] + e2_2[i];	// 40 - 50%
			sum_e3_2[4] = sum_e3_2[4] + e3_2[i];
			sum_E23[4] = sum_E23[4] + E23[i];
			sum_X[4] = sum_X[4] + e2[i];
			sum_Y[4] = sum_Y[4] + e3[i];
			counter[4]+=1;
		}
		else if(Array[i]>11.074 && Array[i]<=11.858){
			sum_e2_2[5] = sum_e2_2[5] + e2_2[i];	// 50 - 60%
			sum_e3_2[5] = sum_e3_2[5] + e3_2[i];
			sum_E23[5] = sum_E23[5] + E23[i];
			sum_X[5] = sum_X[5] + e2[i];
			sum_Y[5] = sum_Y[5] + e3[i];
			counter[5]+=1;
		}
		else if(Array[i]>11.858 && Array[i]<=12.838){
			sum_e2_2[6] = sum_e2_2[6] + e2_2[i];	// 60 - 70%
			sum_e3_2[6] = sum_e3_2[6] + e3_2[i];
			sum_E23[6] = sum_E23[6] + E23[i];
			sum_X[6] = sum_X[6] + e2[i];
			sum_Y[6] = sum_Y[6] + e3[i];
			counter[6]+=1;
		}
		else if(Array[i]>12.838 && Array[i]<=13.622){
			sum_e2_2[7] = sum_e2_2[7] + e2_2[i];	// 70 - 80%
			sum_e3_2[7] = sum_e3_2[7] + e3_2[i];
			sum_E23[7] = sum_E23[7] + E23[i];
			sum_X[7] = sum_X[7] + e2[i];
			sum_Y[7] = sum_Y[7] + e3[i];
			counter[7]+=1;
		}
		else if(Array[i]>13.622){
			sum_e2_2[8] = sum_e2_2[8] + e2_2[i];	// > 80%
			sum_e3_2[8] = sum_e3_2[8] + e3_2[i];
			sum_E23[8] = sum_E23[8] + E23[i];
			sum_X[8] = sum_X[8] + e2[i];
			sum_Y[8] = sum_Y[8] + e3[i];
			counter[8]+=1;
		}
		
	
	}
	
	c18->cd();
	hprof5->GetXaxis()->SetTitle("Multiplicity");
	hprof5->GetYaxis()->SetTitle("npart (# of participating nucleons)");
	hprof5->Draw();
	c18->BuildLegend();
	
	TCanvas *c19 = new TCanvas();
	TGraph* gr5 = new TGraph(10000, N, MArray);
	c19->cd();
	gr5->SetTitle("npart vs Mult S = 0;npart;Mult");
	gr5->SetLineColor(6);
	gr5->SetMarkerColor(8);
	gr5->SetMarkerStyle(24);
	gr5->Draw("ALP");
	
	
	//cout << "Sum of Eccentricity^2" << "  " << "No. of events" << "  " << "RMS" << endl;
	
	/*for(Int_t l = 0; l<9; l++){
		RMS_2[l] = TMath::Sqrt((sum_e2_2[l]/counter[l]));
		RMS_3[l] = TMath::Sqrt((sum_e3_2[l]/counter[l]));
		E23avg[l] = sum_E23[l]/counter[l];
		sum_X_avg[l] = sum_X[l]/counter[l];
		sum_Y_avg[l] = sum_Y[l]/counter[l];
		W[l] = sum_X_avg[l]*sum_Y_avg[l];
		//cout << sum_e2_2[l] << "  " << counter[l] << "  " << curl[l] << endl;
	}
	
	
	c10->cd();
	hprof1->GetXaxis()->SetTitle("Impact Parameter b (in fm)");
	hprof1->GetYaxis()->SetTitle("Eccentricity for n = 2 for S = 0");
	hprof1->Draw();
	c10->BuildLegend();
	
	
	c11->cd();
	hprof2->GetXaxis()->SetTitle("Centrality %");
	hprof2->GetYaxis()->SetTitle("Eccentricity for n = 2 for S = 0");
	hprof2->Draw();
	c11->BuildLegend();
	
	c12->cd();
	hprof3->GetXaxis()->SetTitle("Impact Parameter b (in fm)");
	hprof3->GetYaxis()->SetTitle("Eccentricity for n = 3 for S = 0");
	hprof3->Draw();
	c12->BuildLegend();

	c13->cd();
	hprof4->GetXaxis()->SetTitle("Centrality %");
	hprof4->GetYaxis()->SetTitle("Eccentricity for n = 3 for S = 0");
	hprof4->Draw();
	c13->BuildLegend();
	
	
	TCanvas *c14 = new TCanvas();
	c14->cd();
	
	TGraph* gr1 = new TGraph(9, centr, RMS_2);
	gr1->SetTitle("RMS Eccentricity vs Centrality n = 2 S = 0;Centrality %;RMS Eccentricity");
	gr1->SetLineColor(2);
	gr1->SetMarkerColor(4);
	gr1->SetMarkerStyle(21);
	gr1->Draw("ALP");
	
	TCanvas *c15 = new TCanvas();
	c15->cd();
	
	TGraph* gr2 = new TGraph(9, centr, RMS_3);
	gr2->SetTitle("RMS Eccentricity vs Centrality n = 3 S = 0;Centrality %;RMS Eccentricity");
	gr2->SetLineColor(3);
	gr2->SetMarkerColor(5);
	gr2->SetMarkerStyle(22);
	gr2->Draw("ALP");
	
	TCanvas *c16 = new TCanvas();
	c16->cd();
		
	TGraph* gr3 = new TGraph(9, centr, E23avg);
	gr3->SetTitle(" <e2.e3> S = 0;Centrality %;RMS Eccentricity");
	gr3->SetLineColor(4);
	gr3->SetMarkerColor(6);
	gr3->SetMarkerStyle(23);
	gr3->Draw("ALP");
	
	
	TCanvas *c17 = new TCanvas();
	c17->cd();
		
	TGraph* gr4 = new TGraph(9, centr, W);
	gr4->SetTitle(" <e2>.<e3> S = 0;Centrality %;RMS Eccentricity");
	gr4->SetLineColor(6);
	gr4->SetMarkerColor(8);
	gr4->SetMarkerStyle(24);
	gr4->Draw("ALP");*/
	
	/*TFile* f1 = new TFile("2_S0.root", "CREATE");
	hprof1->Write();	
	hprof2->Write();
	gr1->Write("Graph1");*/
	
	/*TFile* f2 = new TFile("S0.root", "CREATE");
	gr3->Write("Graph3");
	gr4->Write("Graph4");*/	
	//input->Close();
}