#include <TCanvas.h>
#include <TGraph.h>

#include <fstream>
#include <string>
using namespace std;

#include "StandardLabels.C"

void GetParms(double thresh, double &p0, double &p1, double &p2, double &p3);

#include "Angle.h"
#include "Geom.h"

//--------------
// Eres_vs_thresh_2006BeamTest
//--------------
void Eres_vs_thresh_2006BeamTest(void)
{
	double thresholds[7] = {0.1, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0};
	double sigs[7];
	int N = 0;
	for(; N<7; N++){
		double p0, p1, p2, p3;
		GetParms(thresholds[N], p0, p1, p2, p3);
		
		double E = 0.4; // GeV
		double sigma2 = pow(p0*sqrt(E), 2.0) + pow(p1*E, 2.0) + pow(p2, 2.0) + pow(p3*E*log(E), 2.0);
		double sigma_E_over_E = sqrt(sigma2)/E;
		
		sigs[N] = sigma_E_over_E; // already in %
		cout<<"thresh="<<thresholds[N]<<"  sig="<<sigs[N]<<endl;
	}
	
	TGraph *gr = new TGraph(N, thresholds, sigs);
	gr->SetMarkerStyle(22);
	gr->SetMarkerColor(kRed);
	gr->SetMarkerSize(1.5);
	
	gr->Fit("pol1","","", 1.0, 5.0);
	gr->GetFunction("pol1")->SetRange(0.0, 6.0);
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetGrid();
	c1->SetTicks();

	TH2D *axes = new TH2D("axes","Energy Resolution vs. fADC threshold", 100, 0.0, 6.0, 100, 8.5, 10.5);
	axes->SetStats(0);
	axes->SetXTitle("fADC Threshold (MeV)");
	axes->SetYTitle("#sigma_{E}/E at 400MeV (%)");
	axes->Draw();
	gr->Draw("Psame");
	
	TLatex *lab = new TLatex(1.5, 9.7, "fit to x#geq1");
	lab->Draw();

	StandardLabels(axes, Scheme(), AngleStr("#theta_{#gamma}="));

	c1->SaveAs("Eres_vs_thresh_2006BeamTest.png");
	c1->SaveAs("Eres_vs_thresh_2006BeamTest.pdf");
}

//--------------
// GetParms
//--------------
void GetParms(double thresh, double &p0, double &p1, double &p2, double &p3)
{
	p0 = 0.0;
	p1 = 0.0;

	char fname[256];
	sprintf(fname, "%03.1f/Eres_parms.dat", thresh);
	
	ifstream ifs(fname);
	if(!ifs.is_open()){
		cout<<"Couldn't open \""<<fname<<"\"!"<<endl;
		return;
	}

	string scheme;
	double angle;
	double p0_err, p1_err, p2_err, p3_err;

	ifs >> scheme;
	ifs >> angle;
	ifs >> p0 >> p0_err;
	ifs >> p1 >> p1_err;
	ifs >> p2 >> p2_err;
	ifs >> p3 >> p3_err;
}