// WARNING: This script seems to only work the first 
// time it is run. Subsequent invocations cause empty
// histos. Restart ROOT before re-running!


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

void E_res_vs_theta(void)
{
	gROOT->Reset();

	int max_events = 1000000;

	TChain *fine = chain_fine("event");
	TChain *course = chain_course("event");
	
	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	
	TH2D *Eres_fine = new TH2D("Eres_fine", "", 50, 10.0, 110.0, 80, -2.0, 2.0);
	TH2D *Eres_course = new TH2D("Eres_course", "", 50, 10.0, 110.0, 80, -2.0, 2.0);
	
	TCut cut1("");
	
	fine->Project("Eres_fine","(E.recon.p.E()-E.thrown.E())/E.thrown.E():E.thrown.Theta()*57.3",cut1,"",max_events);
	course->Project("Eres_course","(E.recon.p.E()-E.thrown.E())/E.thrown.E():E.thrown.Theta()*57.3",cut1,"",max_events);
	
	Eres_fine->FitSlicesY();
	Eres_course->FitSlicesY();
	
	TH1D *Eres_fine_2 = (TH1D*)gROOT->FindObject("Eres_fine_2");
	TH1D *Eres_course_2 = (TH1D*)gROOT->FindObject("Eres_course_2");
	
	Eres_fine_2->SetLineColor(kRed);
	Eres_fine_2->SetLineWidth(2);
	Eres_course_2->SetLineColor(kBlue);
	Eres_course_2->SetLineWidth(2);

	TH2D *axes = new TH2D("axes","BCAL Energy resolution", 100, 10.0, 110.0, 100, 0.0, 1.0);
	axes->SetStats(0);
	axes->SetYTitle("#sigma((E_{rec}-E_{gen})/E_{gen}) (deg.)");
	axes->SetXTitle("#theta_{#gamma} (GeV)");
	axes->Draw();

	Eres_fine_2->Draw("same");
	Eres_course_2->Draw("same");
	
	TLegend *leg = new TLegend(0.564,0.680, 0.883, 0.872);
	leg->SetFillColor(kWhite);
	leg->AddEntry(Eres_fine_2,"Fine Segmentation");
	leg->AddEntry(Eres_course_2,"Course Segmentation");
	leg->Draw();
	
	char str[256];
	sprintf(str,"%d events   No cuts on primary conversion point.", fine->GetEntries(cut1));
	
	StandardLabels(axes, "0GeV #leq E #leq 4GeV  10^{o}#leq#theta#leq110^{o}", "Single #gamma", str, "Sigmas from Gaussian Fits");
	
	c1->SaveAs("E_res_vs_theta.png");
	c1->SaveAs("E_res_vs_theta.pdf");
}