void A1()
{
	
	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
	 
	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 eccen_cos[12000]; // cos term for each participant 
	Float_t eccen_sin[12000]; // sin term for each participant
	
	Float_t eccen_cos_sum[12000]; // sum of cos terms for each participant 
	Float_t eccen_sin_sum[12000]; // sum of sin terms for each participant 
	
	Float_t eccen_cos_sum_avg[12000]; //eccen_cos_sum divided by # of participapting nucleons for that particular collision
	Float_t eccen_sin_sum_avg[12000]; //eccen_sin_sum divided by # of participapting nucleons for that particular collision
	
	Float_t eccen_num[12000]; // Numerator of eccentricity expression
	Float_t eccen[12000]; //Eccentricity 
	
	
	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
		}
			 
	}	
	
	
	TCanvas *c10 = new TCanvas();
	
	TH1F *myhist = new TH1F("myhist", "Impact Parameter b", 100, 0., 0.); // auto binning
	for (Int_t i=0; i<entries; i++){
		//leaf->GetBranch()->GetEntry(i);	
		//myhist->Fill(leaf->GetValue());
		tree->GetEntry(i);
		myhist->Fill(Imp);
		//cout << i << " : " << Imp << "  " << Array[i] << endl;
	}
	//c10->cd();
	//myhist->Draw();
	
	Int_t TBins = myhist->GetNbinsX();
	//cout << TBins << endl; // 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 = 0;
	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];
	
	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 >= 950 && sum<1000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
			
		}
		if(sum >= 1950 && sum<2000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);			
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 2900 && sum<3000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 3900 && sum<4000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 4900 && sum<5000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 5900 && sum<6050){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 6900 && sum<7000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 7900 && sum<=8050){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		/*if(sum >= 8900 && sum<9000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}
		if(sum >= 9900 && sum<10000){
			binnumber = i;
			bc = myhist->GetXaxis()->GetBinCenter(binnumber);
			bl = myhist->GetXaxis()->GetBinLowEdge(binnumber);
			cout << binnumber << "  " << bc << "  "  << bl << endl;
		}*/
		}
	
	TCanvas *c9 = new TCanvas();
	TCanvas *c11 = new TCanvas();
	
	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", 100, 0, 80);
	
		
	Int_t N[12000];	// For counting the number of participants

	 
	Int_t n = 2;
	Int_t cent = 0;
	
	
	for(Int_t i = 0; i<entries; i++)
		{	
			tree->GetEntry(i);
			N[i] = 0;
			Array[i] = Imp;
			if(Array[i]<=4.704)cent = 10;
			else if(Array[i]>4.704 && Array[i]<=6.664)cent = 20;
			else if(Array[i]>6.664 && Array[i]<=8.232)cent = 30;
			else if(Array[i]>8.232 && Array[i]<=9.604)cent = 40;
			else if(Array[i]>9.604 && Array[i]<=10.78)cent = 50;
			else if(Array[i]>10.78 && Array[i]<=11.76)cent = 60;
			else if(Array[i]>11.76 && Array[i]<=12.544)cent = 70;
			else if(Array[i]>12.544 && Array[i]<=13.524)cent = 80;
			else if(Array[i]>13.524)cent = 90;
			r_2_sum[i] = 0;
			eccen_cos_sum[i] = 0;
			eccen_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
				
					eccen_cos[j]= r_2[j] * TMath::Cos(n * phi_part[j]);
					eccen_sin[j]= r_2[j] * TMath::Sin(n * phi_part[j]);
				
					eccen_cos_sum[i] = eccen_cos_sum[i] + eccen_cos[j];
					eccen_sin_sum[i] = eccen_sin_sum[i] + eccen_sin[j];
				}
		}
		
		r_2_avg[i] = r_2_sum[i]/N[i];
		eccen_cos_sum_avg[i] = eccen_cos_sum[i]/N[i];
		eccen_sin_sum_avg[i] = eccen_sin_sum[i]/N[i];
		eccen_num[i] = TMath::Sqrt((TMath::Power(eccen_cos_sum_avg[i],2)) + TMath::Power(eccen_sin_sum_avg[i],2)); //Numerator term for eccentricity 
		eccen[i] = eccen_num[i]/r_2_avg[i];
		//cout << eccen[i] << endl;
		Array[i] = Imp;
		hprof1->Fill(Array[i],eccen[i]);
		hprof2->Fill(cent,eccen[i]);
	}
	
	c9->cd();
	hprof1->GetXaxis()->SetTitle("Impact Parameter b (in fm)");
	hprof1->GetYaxis()->SetTitle("Eccentricity for n = 2 for S = 0");
	hprof1->Draw();
	c9->BuildLegend();
	
	c11->cd();
	hprof2->GetXaxis()->SetTitle("Centrality %");
	hprof2->GetYaxis()->SetTitle("Eccentricity for n = 2 for S = 0");
	hprof2->Draw();
	c11->BuildLegend();
	
	
	//TFile* f1 = new TFile("2_S0.root", "CREATE");
	//hprof1->Write();
	

	
	//input->Close();
}