#include "../macros/StandardLabels.C"

void Erel_90degrees_clean(const char *who="fine")
{
	//gROOT->Reset();
	
	char fname[256];
	sprintf(fname,"hd_root_%s.root", who);
	TFile *f = new TFile(fname);
	TTree *event = (TTree*)gROOT->FindObject("event");
	
	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	c1->SetGrid();

	double xmin = -1.0; // MeV
	double xmax = +0.4; // MeV


	// Legend
	TLegend *leg = new TLegend(0.171, 0.458, 0.420, 0.801);
	leg->SetFillColor(kWhite);

	// Make plot for several different energy bins
	double energies[]={0.180, 0.250, 0.450, 0.650};
	int colors[] = {kRed,kBlue,kGreen,kMagenta,kBlack};
	int styles[] = {20, 21, 22, 23, 29};
	int Nenergies = 4;
	TH1D* hists[100];
	double ymax = 0.0;
	for(int i=0; i<Nenergies; i++){

		char hname[256];
		sprintf(hname, "Ediff_%03d", (int)(1000.0*energies[i]));
		TH1D *Ediff = new TH1D(hname, "", 100, xmin, xmax);

		Ediff->SetLineColor(colors[i]);
		Ediff->SetLineWidth(2);
		Ediff->SetMarkerColor(colors[i]);
		Ediff->SetMarkerStyle(styles[i]);
	
		char cut[256];
		sprintf(cut, "A.recon.p.P()>0.0 && A.recon.p.P()<10.0 && abs(A.thrown.E()-%f)<0.02", energies[i]);
		event->Project(hname, "(A.recon.p.E() - A.thrown.E())/A.thrown.E()",cut);
		
		hists[i] = Ediff;
		if(Ediff->GetMaximum()>ymax)ymax = Ediff->GetMaximum();
		
		char lab[256];
		sprintf(lab, "%d MeV", (int)(1000.0*energies[i]));
		leg->AddEntry(Ediff, lab);
	}
	
	TH2D *axes = new TH2D("axes", "BCAL Relative Energy Difference for #theta=90^{o}", 100, xmin, xmax, 100, 0.0, ymax*1.05);
	axes->SetStats(0);
	axes->SetXTitle("(E_{recon}-E_{gen})/E_{gen}");
	axes->Draw();
	
	for(unsigned int j=0; j<Nenergies; j++)hists[j]->Draw("same");

	leg->Draw();
	
	StandardLabels(axes, "1M events   No cuts on primary conversion point.", "Single #gamma at 90^{o}", who);

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