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

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

	TCanvas *c1 = new TCanvas();
	c1->SetTicks();
	c1->SetGrid();
	
	char title[256];
	sprintf(title, "BCAL #theta Res. (%s)", who);
	TH2D *axes = new TH2D("axes",title, 100, 0.0, 3.0, 100, 0.0, 6.0);
	axes->SetStats(0);
	axes->SetXTitle("Energy (GeV)");
	axes->SetYTitle("#theta Resolution(degrees)");
	axes->Draw();
		
	// Legend
	TLegend *leg = new TLegend(0.526, 0.119+0.40, 0.879, 0.462+0.40);
	leg->SetFillColor(kWhite);

	
	// Make plot for several different energy bins
	double thetas[]={12.0, 16.0, 20.0, 45.0, 90.0};
	int colors[] = {kRed,kBlue,kGreen,kMagenta,kBlack};
	int styles[] = {20, 21, 22, 23, 29};
	int Nthetas = 5;
	//for(int i=0; i<1; i++){
	for(int i=0; i<Nthetas; i++){
	
		// Project out single theta bin
		int Nbins = 4;
		int xbin = theta_diff->GetXaxis()->FindBin(thetas[i]);
		int xbinstart = xbin - (Nbins-1)/2;
		int xbinend = xbinstart + Nbins-1;
		if(xbinstart<1)xbinstart=1;
		if(xbinend>theta_diff->GetNbinsX())xbinend=theta_diff->GetNbinsX();
		double ystart = theta_diff->GetXaxis()->GetBinLowEdge(xbinstart);
		double yend = theta_diff->GetXaxis()->GetBinLowEdge(xbinend+1);
		
		theta_diff->GetXaxis()->SetRange(xbinstart, xbinend);
		TH2D *theta_diff_zy = theta_diff->Project3D("zy");
		//theta_diff_zy->FitSlicesY();
		//TH1D *theta_diff_zy_2 = (TH1D*)gROOT->FindObject("theta_diff_zy_2");
		TH1D *theta_diff_zy_2 = MyFitSlicesY(theta_diff_zy);

		char hname[256];
		sprintf(hname, "theta%05d_px", thetas[i]*100.0);
		TH1D *theta_res = (TH1D*)theta_diff_zy_2->Clone(hname);

		theta_res->SetLineColor(colors[i]);
		theta_res->SetMarkerColor(colors[i]);
		theta_res->SetMarkerStyle(styles[i]);
	
		theta_res->Draw("PCLsame");

		char lab[256];
		double theta = thetas[i];
		double Lfrac = (65.0-17.0+65.0/tan(theta/57.3))/390.0;
		sprintf(lab,"%2.0f^{o}- %2.0f^{o} (z=%4.02fL)", ystart, yend, Lfrac);
		leg->AddEntry(theta_res, lab);
	}

	leg->Draw();
	
	StandardLabels(axes, "20M events   No cuts on primary conversion point.", "Single #gamma", who);

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

TH1D* MyFitSlicesY(TH2D *h)
{
	// Fit slices, giving better starting values for the fit
	// so it will converge on the correct peak. This is needed
	// since for the low energy theta_diff spectra being fit to
	// a single Gaussian, there was a small shoulder that the fit 
	// occasionally converged on using TH2D::FitSlices.

cout<<__LINE__<<endl; // this line prevents the macro from crashing !!?!?

	int NbinsY = h->GetNbinsY();
	int NbinsX = h->GetNbinsX();
	
	const char *basename = h->GetName();
	char hname[256];
	sprintf(hname, "%s_2", basename);
	
	// Make 1D histo able to hold sigmas
	TH1D *h_sigma = h->ProjectionX(hname);
	h_sigma->Reset();
	
	TF1 *fun = new TF1("fun", "[0]*TMath::Gaus(x, [1], [2])");
	
	// Loop over Y bins
	for(int ixbin=NbinsX; ixbin>0; ixbin--){
		TH1D *h_py = h->ProjectionY("_py", ixbin, ixbin);
		
		int bin_max = h_py->GetMaximumBin();
		double xmax = h_py->GetBinCenter(bin_max);
		
		if(ixbin==NbinsX){
			fun->SetParameter(0, h_py->GetMaximum());
			fun->SetParameter(1, xmax);
			fun->SetParameter(2, 0.1);
		}

		h_py->Fit(fun, "0Q");
		//c2->Update();
		
		double sigma = fun->GetParameter(2);
		//cout<<"E="<<h_sigma->GetBinCenter(ixbin)<<"  sigma="<<sigma<<endl;
		h_sigma->SetBinContent(ixbin, fabs(sigma));
	}

	cout<<"Fit sigmas for \""<<h_sigma->GetName()<<"\""<<endl;
	
	return h_sigma;
}