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

// CDF = "Cumulative Distribution Function" (i.e. integral fraction)

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

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

	int max_events = 2000000;

	TChain *fine = chain_fine("event");
	TChain *course = chain_course("event");
	
	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	
	TH2D *thetares_fine = new TH2D("thetares_fine", "", 50, 0.0, 2.0, 120, -30.0, 30.0);
	TH2D *thetares_course = new TH2D("thetares_course", "", 50, 0.0, 2.0, 120, -30.0, 30.0);
	
	TCut cut1("abs(E.recon.p.E()-E.thrown.E())/E.thrown.E()<0.2");
	
	fine->Project("thetares_fine","(E.recon.p.Theta()-E.thrown.Theta())*57.3:E.thrown.E()",cut1,"",max_events);
	course->Project("thetares_course","(E.recon.p.Theta()-E.thrown.Theta())*57.3:E.thrown.E()",cut1,"",max_events);
	
	TH1D *thetares_fine_2 = Make95CLHisto(thetares_fine);
	TH1D *thetares_course_2 = Make95CLHisto(thetares_course);
	
	thetares_fine_2->SetMarkerColor(kRed);
	thetares_fine_2->SetLineColor(thetares_fine_2->GetMarkerColor());
	thetares_fine_2->SetMarkerStyle(21);
	thetares_course_2->SetMarkerColor(kBlue);
	thetares_course_2->SetLineColor(thetares_course_2->GetMarkerColor());
	thetares_course_2->SetMarkerStyle(22);

	TH2D *axes = new TH2D("axes","BCAL Polar Angle Resolution", 100, 0.0, 2.0, 100, 0.0, 50.0);
	axes->SetStats(0);
	axes->SetYTitle("95% CDF of Voigt fit to (#theta_{rec}-#theta_{gen}) (deg.)");
	axes->SetXTitle("E_{#gamma} (GeV)");
	axes->Draw();

	thetares_fine_2->Draw("Psame");
	thetares_course_2->Draw("Psame");
	
	TLegend *leg = new TLegend(0.564,0.680, 0.883, 0.872);
	leg->SetFillColor(kWhite);
	leg->AddEntry(thetares_fine_2,"Fine Segmentation");
	leg->AddEntry(thetares_course_2,"Course Segmentation");
	leg->Draw();
	
	char str[256];
	sprintf(str,"2M events   No cuts on primary conversion point.");
	
	StandardLabels(axes, "0GeV #leq E #leq 4MeV  10^{o}#leq#theta#leq110^{o}", "Single #gamma", "2M events   No cuts on primary conversion point","Cut on |E_{recon}-E_{gen}|/E_{gen}<0.2");
	
	c1->SaveAs("theta_res_95CDF_vs_E.png");
	c1->SaveAs("theta_res_95CDF_vs_E.pdf");
}

TH1D* Make95CLHisto(TH2D *h)
{
	char hname[256];
	sprintf(hname, "%s_2", h->GetName());
	
	TAxis *xaxis = h->GetXaxis();
	TAxis *yaxis = h->GetYaxis();
	Int_t Nbins = xaxis->GetNbins();
	TH1D *h_2 = new TH1D(hname, "", Nbins, xaxis->GetXmin(), xaxis->GetXmax());
	TH1D *tmp_py = new TH1D("tmp_py", "", yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax());

	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");
	
	for(int bin=1; bin<=Nbins; bin++){
		tmp_py->Reset();
		h->ProjectionY("tmp_py", bin, bin);
		

		voigt->SetParameter(0, tmp_py->GetMaximum());
		voigt->SetParameter(1, 0.0);
		voigt->SetParameter(2, 0.3);
		voigt->SetParameter(3, 5.0);

		double width = 0.0;
		if(h->GetEntries()>10.0)
		{
			tmp_py->Fit(voigt);
			width = FindWidth(tmp_py->GetFunction("voigt"), 0.95);
		}

		h_2->SetBinContent(bin, width);
		h_2->SetBinError(bin, 0.0);
	}
	
	return h_2;
}

double FindWidth(TF1 *fun, double conf_level)
{
	if(!fun)return 0.0;

	// 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(x0,x1)/norm;
	double I3=fun->Integral(x0,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;
}