// 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 theta_res_90deg(void)
{
	gROOT->Reset();

	int max_events = 2000000;

	TChain *fine = chain_fine("event");
	TChain *course = chain_course("event");
	
	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	
	TH1D *thetares_fine = new TH1D("thetares_fine", "", 100, -30.0, 30.0);
	TH1D *thetares_course = new TH1D("thetares_course", "", 100, -30.0, 30.0);
	
	TCut cut1("E.thrown.E()>0.1 && E.thrown.E()<0.15 && abs(E.thrown.Theta()*57.3-90.0)<10.0");
	
	fine->Project("thetares_fine","(E.recon.p.Theta()-E.thrown.Theta())*57.3",cut1,"",max_events);
	course->Project("thetares_course","(E.recon.p.Theta()-E.thrown.Theta())*57.3",cut1,"",max_events);
	
	TF1 *voigt = new TF1("voigt","[0]*TMath::Voigt(x-[1], [2], [3])", -30.0, 30.0);
	voigt->SetParName(0, "Amp");
	voigt->SetParName(1, "x0");
	voigt->SetParName(2, "sigma");
	voigt->SetParName(3, "gamma");

	voigt->SetParameter(0, thetares_fine->GetMaximum());
	voigt->SetParameter(1, 0.0);
	voigt->SetParameter(2, 2.0);
	voigt->SetParameter(3, 2.0);
	thetares_fine->Fit(voigt);
	thetares_course->Fit(voigt);
	
	thetares_fine->SetLineColor(kRed);
	thetares_fine->SetLineWidth(2);
	thetares_fine->GetFunction("voigt")->SetLineColor(thetares_fine->GetLineColor());
	thetares_course->SetLineColor(kBlue);
	thetares_course->SetLineWidth(2);
	thetares_course->GetFunction("voigt")->SetLineColor(thetares_course->GetLineColor());

	thetares_fine->SetTitle("BCAL Polar angle Resolution");
	thetares_fine->SetXTitle("#theta_{rec}-#theta_{gen}(deg.)");
	thetares_fine->SetYTitle("Events/0.2 deg.");
	thetares_fine->SetStats(0);
	
	thetares_fine->Draw();
	thetares_course->Draw("same");
	
	//TH2D *axes = new TH2D("axes","BCAL Reconstruction Efficiency", 100, 0.0, 110.0, 100, 0.0, 1.05);
	//thetares_fine->SetStats(0);
	//thetares_fine->SetXTitle("#theta_{rec}-#theta_{gen} (degrees)");
	//axes->SetYTitle("Events/0.1 deg.");
	//axes->Draw();
	
	TLegend *leg = new TLegend(0.564,0.680, 0.883, 0.872);
	leg->SetFillColor(kWhite);
	leg->AddEntry(thetares_fine,"Fine Segmentation");
	leg->AddEntry(thetares_course,"Course Segmentation");
	leg->Draw();

	char events_lab[256]="";
	sprintf(events_lab,"%d events   No cuts on primary conversion point.", fine->GetEntries(cut1));
	
	// Make stats box
	double x1 = -29.0;
	double x2 = x1 + 6.0;
	double y1 = 0.4*thetares_fine->GetMaximum();
	double y2 = 0.95*thetares_fine->GetMaximum();
	DrawStatsLabels( x1, y1, x2, y2);
	x1 = x2;
	x2 = x1 + 6.0;
	DrawStats("fine", thetares_fine->GetFunction("voigt"), x1, y1, x2, y2);
	x1 = x2;
	x2 = x1 + 8.0;
	DrawStats("course", thetares_course->GetFunction("voigt"), x1, y1, x2, y2);

	
	StandardLabels1D(thetares_fine, "100MeV #leq E #leq 150MeV  80^{o}#leq#theta#leq110^{o}", "Single #gamma", events_lab);
	
	c1->SaveAs("theta_res_90deg.png");
	c1->SaveAs("theta_res_90deg.pdf");
}

double FindWidth(TF1 *fun, double conf_level)
{
	// Integrate out to 1 million
	double x0 = fun->GetParameter(1);
	double norm = fun->Integral(x0, 1000000.0);
	double x1 = 0.0;
	double x2 = 25.0;
	double x3 = 1000.0;
	double I1=fun->Integral(0.0,x1)/norm;
	double I3=fun->Integral(0.0,x3)/norm;
	for(int i=0; i<100; i++){
		double I2=fun->Integral(x0, x2)/norm;
		
		if(I2>conf_level){
			I3 = I2;
			x3 = x2;
		}else{
			I1 = I2;
			x1 = x2;
		}
		if(fabs(I1-conf_level)<0.0001){
			x2=x1;
			break;
		}else if(fabs(I2-conf_level)<0.0001){
			break;
		}else if(fabs(I3-conf_level)<0.0001){
			x2 = x3;
			break;
		}
		
		x2 = (x1+x3)/2.0;
		
	}
	
	return x2-x0;
}

void DrawStatsLabels(double x1, double y1, double x2, double y2)
{
	TPaveStats *pav = new TPaveStats(x1,y1,x2,y2);
	pav->AddText("  ");
	pav->AddText("#sigma");
	pav->AddText("#Gamma");
	pav->AddText("95%CL");
	pav->Draw();
}


void DrawStats(const char *lab, TF1 *f, double x1, double y1, double x2, double y2)
{
	char str[256];

	TPaveStats *pav = new TPaveStats(x1,y1,x2,y2);
	pav->AddText(lab);
	
	double width = FindWidth(f, 0.95);

	sprintf(str, "%3.2f", f->GetParameter(2));
	pav->AddText(str);
	sprintf(str, "%3.2f", f->GetParameter(3));
	pav->AddText(str);
	sprintf(str, "%3.2f", width);
	pav->AddText(str);
	pav->Draw();
}