#include <TROOT.h>
#include <TH1.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TColor.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TLatex.h>
#include <TProfile.h>

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

void GetPoints(TH2D *h, TGraphErrors* &gres, TGraphErrors* &gmean);

#include "Angle.h"

void HDGeant_leakage(void)
{
	TFile *f = new TFile("leakage2.root");
	TTree *leaks = (TTree*)gROOT->FindObject("leaks");

	TH2D *dE_vs_E =new TH2D("dE_vs_E", "HDGeant missing energy, Absolute", 50, 0.001, 2.0, 200, 0.0, 100.0);
	TH2D *dE_vs_E_over_E =new TH2D("dE_vs_E_over_E", "HDGeant missing energy, Relative", 50, 0.001, 2.0, 200, 0.0, 0.1);
	
	leaks->Project("dE_vs_E", "1000.0*(Egen-(E_in_bcal+E_out_bcal+E_leave_mother)):Egen");
	leaks->Project("dE_vs_E_over_E", "(Egen-(E_in_bcal+E_out_bcal+E_leave_mother))/Egen:Egen");
	
	dE_vs_E->SetStats(0);
	dE_vs_E->SetXTitle("Energy, generated (GeV)");
	dE_vs_E->SetYTitle("E_{gen}-#Sigma dE_{step} (MeV)");

	dE_vs_E_over_E->SetStats(0);
	dE_vs_E_over_E->SetXTitle("Energy, generated (GeV)");
	dE_vs_E_over_E->SetYTitle("(E_{gen}-#Sigma dE_{step})/E_{gen}");	
	
	dE_vs_E->SetMarkerStyle(8);
	dE_vs_E->SetMarkerSize(0.5);
	dE_vs_E->SetMarkerColor(kMagenta);
	dE_vs_E_over_E->SetMarkerStyle(8);
	dE_vs_E_over_E->SetMarkerSize(0.5);
	dE_vs_E_over_E->SetMarkerColor(kMagenta);
	
	TCanvas *c1 = new TCanvas("c1", "", 800, 800);
	c1->Divide(1,2);
	
	c1->cd(1);
	gPad->SetTicks();
	gPad->SetGrid();
	dE_vs_E->Draw();

	c1->cd(2);
	gPad->SetTicks();
	gPad->SetGrid();
	dE_vs_E_over_E->Draw();

	//StandardLabels(axes, Scheme(), "", AngleStr("#theta_{#gamma}="));
	
	c1->SaveAs("HDGeant_leakage.pdf");
	c1->SaveAs("HDGeant_leakage.png");
}

//------------------
// GetPoints
//------------------
void GetPoints(TH2D *h, TGraphErrors* &gres, TGraphErrors* &gmean)
{
	int min_entries = 1000;

	// Project onto the x-axis so we can find the limits for
	// for the projection on the y-axis based on number of events
	TH1D *h_px = h->ProjectionX("_px", 1);

	// Loop over (uneven) bins
	int Npoints = 0;
	vector<double> x;
	vector<double> y;
	vector<double> y_mean;
	vector<double> xerr;
	vector<double> yerr;
	vector<double> y_mean_err;
	int Nbins = h->GetNbinsX();
	for(int start_bin=1; start_bin<=Nbins; ){

		// Find limits to give us min_entries entries
		int end_bin = start_bin;
		while(h_px->Integral(start_bin, end_bin)<min_entries && end_bin<Nbins){
			++end_bin;
		}
		
		// If only a few events are left, go ahead and include them here.
		if(end_bin<Nbins && h_px->Integral(end_bin, Nbins)<min_entries)end_bin=Nbins;

		// Project onto y-axis
		cout<<"projecting bins "<<start_bin<<" to "<<end_bin<<"  Nentries="<<h_px->Integral(start_bin, end_bin)<<endl;
		TH1D *h_py = h->ProjectionY("_py",start_bin, end_bin);

		if(h_py->Integral()>=20){
		
			// Fit to Gaussian
			//h_py->Fit("gaus","0Q");
			//c1->Update();

			// Get fit results
			double mean = h_py->GetBinCenter(h_py->GetMaximumBin());
			double sigma = h_py->GetRMS();
			double mean_err = sigma/sqrt(h_py->Integral());
			double sigma_err = 0.0;

			// Find average fADC value for hits in this bin range
			double fADC = 0.0;
			double norm = 0.0;
			for(int bin=start_bin; bin<=end_bin; bin++){
				double weight = h_px->GetBinContent(bin);
				fADC += h_px->GetBinCenter(bin)*weight;
				norm += weight;
			}
			fADC /= norm;
		
			// Add floor term to error obtained from Gaussian fit
			double epsilon = 0.001; // GeV
			mean_err = sqrt(mean_err*mean_err + epsilon*epsilon);
			//sigma_err = sqrt(sigma_err*sigma_err + epsilon*epsilon);

			// Error on reconstructed energy is difficult to estimate here.
			// It is definitely worse than the final algorithm will yield.
			// The fact that the "reconstructed" energy extends out to 2.5GeV
			// when the generated data set only went to 2GeV says it is 
			// considerably worse than 5.5%/sqrt(E). For lack of a better
			// idea, I use the RMS of the fADC values about the mean
			// in the current range.
			double fADC_rms = 0.0;
			double norm_rms = 0.0;
			for(int bin=start_bin; bin<=end_bin; bin++){
				double weight = h_px->GetBinContent(bin);
				fADC_rms += pow(h_px->GetBinCenter(bin)-fADC, 2.0)*weight;
				norm_rms += weight;
			}
			fADC_rms /= norm_rms;
			fADC_rms = sqrt(fADC_rms);

			x.push_back(fADC);
			y.push_back(sigma);
			y_mean.push_back(mean);
			xerr.push_back(fADC_rms);
			yerr.push_back(sigma_err);
			y_mean_err.push_back(sigma_err);
			Npoints++;
		}

		// Increment for next iteration
		start_bin=end_bin+1;
	}
	
	gres = new TGraphErrors(Npoints, &x[0], &y[0], &xerr[0], &yerr[0]);
	gres->SetMarkerColor(kCyan);
	gres->SetLineColor(kCyan);
	gres->SetLineWidth(2);

	gmean = new TGraphErrors(Npoints, &x[0], &y_mean[0], &xerr[0], &y_mean_err[0]);
	gmean->SetMarkerColor(kBlue);
	gmean->SetMarkerStyle(22);
	gmean->SetLineColor(kBlue);
	
}