// 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 eff_vs_theta(const char *who="fine")
{
	gROOT->Reset();

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

	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	c1->SetGrid();
	
	char title[256];
	sprintf(title, "BCAL Recon. Efficiency (%s)", who);
	TH2D *axes = new TH2D("axes",title, 100, 10.0, 110.0, 100, 0.0, 1.05);
	axes->SetStats(0);
	axes->SetXTitle("#theta (degrees)");
	axes->SetYTitle("Efficiency");
	axes->Draw();

	// Fix cut at 4 sigma
	TAxis *zaxis = eff_numerator->GetZaxis();
	int max_bin = zaxis->FindBin(4.0);
	zaxis->SetRange(1, max_bin);
	TH2D *eff = (TH2D*)eff_numerator->Project3D("yx");
	eff->Divide(eff_denominator);
	
	// 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.100, 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 ybin = eff->GetYaxis()->FindBin(energies[i]);
		char hname[256];
		sprintf(hname, "E%4.2f_px", energies[i]);
		TH1D *my_eff = (TH1D*)eff->ProjectionX(hname,ybin, ybin);
		my_eff->SetLineColor(colors[i]);
		my_eff->SetMarkerColor(colors[i]);
		my_eff->SetMarkerStyle(styles[i]);
	
		my_eff->Draw("PCLsame");

		char lab[256];
		double E = energies[i];
		double E_binwidth = eff->GetYaxis()->GetBinWidth(ybin);
		const char *units = "GeV";
		sprintf(lab,"%5.3f-%5.3f %s", E, E + E_binwidth, units);
		if(E<1.0){
			E *= 1000.0;
			E_binwidth *= 1000.0;
			units = "MeV";
			sprintf(lab,"%d-%d %s", (int)E, (int)(E + E_binwidth), units);
		}
		leg->AddEntry(my_eff, lab);
	}

	leg->Draw();

	char lab[256];
	sprintf(lab,"Cut on E_{recon} within 4#sigma of 5.5%% #times #sqrt{E_{thrown}}");
	StandardLabels(axes, "20M events   No cuts on primary conversion point.", "Single #gamma", "", lab);

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