#include "StandardLabels.C"

void Eres_vs_E(void)
{
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGridy();

	// Make list of energies
	vector<double> energies;
	energies.push_back(0.100);
	energies.push_back(0.250);
	energies.push_back(0.500);
	energies.push_back(0.750);
	energies.push_back(1.000);
	energies.push_back(1.250);
	energies.push_back(1.500);
	energies.push_back(1.750);
	energies.push_back(2.000);
	energies.push_back(4.000);
	
	// Make function for energy resolution function
	TF1 *fun = new TF1("fun", "sqrt([0]*[0]/x + [1]*[1] +[2]*[2]/x/x)");
	fun->SetParName(0, "stochastic");
	fun->SetParName(1, "floor");
	fun->SetParName(2, "noise");
	
	// Get Data points
	vector<string> labels;
	TGraphErrors *degrees12 = MakeTGraph(energies, 12.0, labels, kRed);
	TGraphErrors *degrees15 = MakeTGraph(energies, 15.0, labels, kGreen);
	TGraphErrors *degrees20 = MakeTGraph(energies, 20.0, labels, kBlue);
	TGraphErrors *degrees30 = MakeTGraph(energies, 30.0, labels, kMagenta);

	TH2D *axes = new TH2D("axes", "BCAL Energy resolution from M.C.", 100, 0.0, 4.250, 100, 0.0, 0.60);
	axes->SetStats(0);
	axes->SetXTitle("Incident Particle Energy (GeV)");
	axes->SetYTitle("#sigma_{#deltaE}/E");
	axes->Draw();

	degrees12->Draw("PE");
	degrees15->Draw("PE");
	degrees20->Draw("PE");
	degrees30->Draw("PE");
	
	// Make legend
	TLegend *leg = new TLegend(0.428, 0.414, 0.891, 0.888);
	leg->SetFillColor(kWhite);
	leg->AddEntry(degrees12, labels[0].c_str());
	leg->AddEntry(degrees15, labels[1].c_str());
	leg->AddEntry(degrees20, labels[2].c_str());
	leg->AddEntry(degrees30, labels[3].c_str());
	leg->Draw();
	
	StandardLabels(axes);

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

TGraphErrors* MakeTGraph(vector<double> &energies, double theta, vector<string> &labels, int color)
{
	vector<double> Eres;
	vector<double> Eres_err;

	for(unsigned int i=0; i<energies.size(); i++){
		char fname[256];
		sprintf(fname, "../%1.3fGeV/%ddegrees/hd_root.root", energies[i], (int)theta);
		//sprintf(fname, "../%1.3fGeV/%ddegrees/hd_root_no_darkhits.root", energies[i], (int)theta);
		cout<<"Opening \""<<fname<<"\""<<endl;
		TFile *f = new TFile(fname);
		TH2D *Ethrown_vs_Erecon = (TH2D*)gROOT->FindObject("Ethrown_vs_Erecon");
		
		TH1D *Ethrown_vs_Erecon_px = Ethrown_vs_Erecon->ProjectionX();
		Ethrown_vs_Erecon_px->Fit("gaus", "", "", energies[i]-1.0, energies[i]+2.0);
		TF1 *g = Ethrown_vs_Erecon_px->GetFunction("gaus");
		
		double mean = g->GetParameter(1);
		double sigma = g->GetParameter(2);
		double sigma_err = g->GetParError(2);
		
		// Convert to GeV
		double fac = energies[i]/mean;
		double sigma_GeV = sigma*fac;
		double sigma_err_GeV = sigma_err*fac;
		double sigma_fraction = sigma_GeV/energies[i];
		double sigma_err_fraction = sigma_err_GeV/energies[i];
		
		cout<<"mean="<<mean<<endl;
		cout<<"sigma="<<sigma<<" ("<<sigma_GeV<<"GeV, "<<sigma_fraction*100.0<<"%)"<<endl;
		cout<<"sigma_err="<<sigma_err<<endl;

		Eres.push_back(sigma_fraction);
		Eres_err.push_back(sigma_err_fraction);
	}
	
	TGraphErrors *tgraph = new TGraphErrors(energies.size(), &energies[0], &Eres[0], 0, &Eres_err[0]);

	tgraph->SetMarkerStyle(20);
	tgraph->SetMarkerColor(color);
	tgraph->SetLineColor(color);
	tgraph->SetLineWidth(2);
	tgraph->SetFillColor(kWhite);
	
	// Set starting parameters for fit
	TF1 *fun = (TF1*)gROOT->FindObject("fun");
	fun->SetParameter(0, 0.042);  // Stochastic term
	fun->SetParameter(1, 0.013);  // Floor term
	fun->SetParameter(2, 0.000);  // Noise term
	
	//fun->FixParameter(0, 0.042);
	fun->FixParameter(1, 0.013);
	//fun->FixParameter(2, 0.0);
	
	tgraph->Fit("fun", "", "", theta>12.0 ? 0.050:0.150, 4.50);
	//tgraph->Fit("fun", "", "", 0.050, 4.50);
	TF1 *myfun = tgraph->GetFunction("fun");
	myfun->SetLineColor(color);
	myfun->SetLineStyle(1);
	myfun->SetLineWidth(2);
	
	// Create Latex string representing this function
	char str[256];
	sprintf(str, "%d^{o}     #frac{%2.1f%%}{#sqrt{E}} #oplus %2.1f%% #oplus #frac{%2.1f%%}{E}"
		, (int)theta
		, fabs(myfun->GetParameter(0))*100.0
		, fabs(myfun->GetParameter(1))*100.0
		, fabs(myfun->GetParameter(2))*100.0);
	labels.push_back(str);

	// Append chi2/dof result to output file
	cout<<"Chisq="<<myfun->GetChisquare()<<" NDF="<<myfun->GetNDF()<<endl;
	ofstream ofs("chisq.dat", ios_base::app);
	ofs<<theta<<" "<<myfun->GetChisquare()/myfun->GetNDF()<<endl;

	return tgraph;
}