// This macro is to be used on a root file produced with the
// event.C selector.


#include "chain_fine.h"
#include "chain_course.h"
#include "StandardLabels.C"

void Eres_vs_theta(const char *who="fine")
{
	gROOT->Reset();
	//gStyle->SetPadRightMargin(0.15);

	char fname[256];
	sprintf(fname,"event_%s.root", who);
	TFile *f = new TFile(fname);
	TH3D *E_diff = (TH3D*)gROOT->FindObject("E_diff");	

	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	c1->SetGrid();
	
	char title[256];
	sprintf(title, "BCAL Energy Res. (%s)", who);
	TH2D *axes = new TH2D("axes",title, 100, 10.0, 110.0, 100, 0.0, 60.0);
	axes->SetStats(0);
	axes->SetXTitle("#theta (degrees)");
	axes->SetYTitle("Energy Resolution(%)");
	axes->Draw();
		
	// Legend
	TLegend *leg = new TLegend(0.526, 0.119, 0.879, 0.462);
	leg->SetFillColor(kWhite);

	
	// Make plot for several different energy bins
	double energies[]={0.020, 0.060, 0.200, 0.500, 1.0};
	int colors[] = {kRed,kBlue,kGreen,kMagenta,kBlack};
	int styles[] = {20, 21, 22, 23, 29};
	int Nenergies = 5;
	for(int i=0; i<Nenergies; i++){
	
		// Project out single energy bin
		int Nbins = 4;
		int ybin = E_diff->GetYaxis()->FindBin(energies[i]);
		int ybinstart = ybin - (Nbins-1)/2;
		int ybinend = ybinstart + Nbins-1;
		double ystart = E_diff->GetYaxis()->GetBinLowEdge(ybinstart);
		double yend = E_diff->GetYaxis()->GetBinLowEdge(ybinend+1);
		
		E_diff->GetYaxis()->SetRange(ybinstart, ybinend);
		TH2D *E_diff_zx = E_diff->Project3D("zx");
		E_diff_zx->FitSlicesY();
		TH1D *E_diff_zx_2 = (TH1D*)gROOT->FindObject("E_diff_zx_2");

		char hname[256];
		sprintf(hname, "E%05d_px", energies[i]*100.0);
		TH1D *Eres = (TH1D*)E_diff_zx_2->Clone(hname);

		Eres->Scale(100.0); // convert to percent

		Eres->SetLineColor(colors[i]);
		Eres->SetMarkerColor(colors[i]);
		Eres->SetMarkerStyle(styles[i]);
	
		Eres->Draw("PCLsame");

		char lab[256];
		double E = energies[i];
		double E_binwidth = E_diff->GetYaxis()->GetBinWidth(ybin);
		const char *units = "GeV";
		sprintf(lab,"%5.3f-%5.3f %s", ystart, yend, units);
		if(E<1.0){
			E *= 1000.0;
			E_binwidth *= 1000.0;
			units = "MeV";
			sprintf(lab,"%d-%d %s", (int)(1000.0*ystart), (int)(1000.0*yend), units);
		}
		leg->AddEntry(Eres, lab);
		
		double y = Eres->GetBinContent(Eres->GetXaxis()->FindBin(55.0)) + 3.0 - 1.75*energies[i];
		TLatex *label = new TLatex(55.0, y, lab);
		label->SetTextAlign(22);
		label->SetTextSize(0.035);
		label->SetTextColor(colors[i]);
		label->Draw();
	}

	//leg->Draw();
	
	StandardLabels(axes, "20M events   No cuts on primary conversion point.", "Single #gamma");

	sprintf(fname,"Eres_vs_theta_%s.png", who);
	c1->SaveAs(fname);
	sprintf(fname,"Eres_vs_theta_%s.pdf", who);
	c1->SaveAs(fname);
}