#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 <TMath.h>
#include <TFitResult.h>

#include <unistd.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

#include "StandardLabels.C"

#include "E_calibrate.cc"

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

void GetPoints(TH2D *h, TGraphErrors* &gres, TGraphErrors* &gmean);
Double_t asymmetric_gaus(Double_t *xptr, Double_t *par);

//------------------------------------
// Copied from JEventProcessor_bcal_timing.cc
typedef struct{
	int event;
	int layer;
	int sector;
	int fADC_up;
	int fADC_dn;
	float Etot;
	float geometric_mean;
	float tup;
	float tdn;
	float tup_corrected;
	float tdn_corrected;
	float theta_thrown;
	float E_thrown;
}HIT_t;
//------------------------------------

//-----------------
// Eres_vs_E
//-----------------
void Eres_vs_E(double fADC_thresh_MeV=1.0)
{

	TFile *f = new TFile("hd_root.root");
	f->cd();
	
	TTree *tree = (TTree*)gROOT->FindObject("tree");
	HIT_t hit;
	HIT_t *hitptr = &hit;
	TBranch *branch = tree->GetBranch("T");
	branch->SetAddress(hitptr);

	TH2D *h = new TH2D("h", "", 200, 0.150, 0.650, 200, -0.5, 0.5);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();

	// Loop over all events and fill histogram. We do this
	// rather than TTree::Project since we need to call the 
	// BCAL_tdiff_res for each entry
	int Nentries = branch->GetEntries();
	double E_thrown = 0.0;
	double E_tot = 0.0;
	int last_event=1;
	for(int i=1; i<=Nentries; i++){
		tree->GetEntry(i);
		
		// Fill histogram on event boundaries and reset for next event
		if(hit.event!=last_event){
		
			E_tot = BCAL_fADC2E(E_tot); // convert to GeV

			h->Fill(E_thrown, E_tot - E_thrown);
			
			E_thrown = 0.0;
			E_tot = 0.0;
			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>MaxSummedLayers())continue; //

		// Apply fADC threshold
		double mV_per_MeV = 5.0833; // see bcal_timing/DBCALTimeSpectrum.cc::Convolute
		       mV_per_MeV *= 0.346/0.668; // 5.0833 value uses new fiber light output (145 photons/MeV)
				                            // scale back down to fibers used for beam test (75 photons/MeV)
		       mV_per_MeV /= 5.0; // This scales down the second stage TDC amplification so we are in units of the fADC leg
		double LSB_pC_fADC = 2.0/50.0*4.0E3/4096.0; // see top of bcal_timing/DBCALTimeSpectrum.cc
		double SiPM_pulse_integral_pC = 1200.0;
		double SiPM_pulse_peak_mV = 2293.0;
		double pC_per_MeV = mV_per_MeV * (SiPM_pulse_integral_pC/SiPM_pulse_peak_mV);
		double LSB_MeV_fADC = LSB_pC_fADC/pC_per_MeV;
		double fADC_threshold = fADC_thresh_MeV/LSB_MeV_fADC; // threshold for 3.0 attenuated MeV in fADC counts
		       fADC_threshold = floor(0.5 + fADC_threshold); // integerize
		if(hit.fADC_up<fADC_threshold)continue;
		if(hit.fADC_dn<fADC_threshold)continue;
		//if(hit.tup_corrected<-100)continue;
		//if(hit.tdn_corrected<-100)continue;
		
		// Accumulate info for this event
		double geomn = hit.geometric_mean;
		double E = geomn;
		E_tot += E;
		E_thrown = hit.E_thrown;
	}
	
	TGraphErrors *gres;
	TGraphErrors *gmean;
	GetPoints(h, gres, gmean);
	//return;

	// Fit to find sigma of delta E
	TF1 *fun = new TF1("fun", "sqrt(pow([0]*sqrt(x),2.0) + pow([1]*x,2.0) + pow([2],2.0) + pow([3]*x*log(x),2.0))");
	fun->SetParameter(0, 0.00);
	fun->SetParameter(1, 0.0);
	fun->SetParameter(2, 0.0);
	fun->SetParameter(3, 0.00);
	gres->Fit(fun, "0SE", "", 0.150, 0.650);

	// Now, convert function and data to DeltaE/E fo plotting
	// I'm not sure if this is necessary, but it seems most
	// resolution plots are of sigmaE/E. The fitting to find
	// sigmaE though seems like it should be done to DeltaE
	for(int i=0; i<gres->GetN(); i++){
		double E, DeltaE;
		gres->GetPoint(i, E, DeltaE);
		DeltaE /= E;
		gres->SetPoint(i, E, DeltaE);
	}

	TF1 *fun2 = new TF1("fun2", "sqrt(pow([0]/sqrt(x),2.0) + pow([1],2.0) + pow([2]/x,2.0) + pow([3]*log(x),2.0))", 0.0, 2.0);
	fun2->SetParameter(0, fabs(fun->GetParameter(0)));
	fun2->SetParameter(1, fabs(fun->GetParameter(1)));
	fun2->SetParameter(2, fabs(fun->GetParameter(2)));
	fun2->SetParameter(3, fabs(fun->GetParameter(3)));
	//gres->Fit(fun2);

	// From Blake's analysis
	TF1 *nim = new TF1("nim", "sqrt(pow([0]/sqrt(x),2.0) + pow([1],2.0))", 0.0, 2.0);
	nim->SetParameter(0, 0.054);
	nim->SetParameter(1, 0.023);
	nim->SetLineColor(kRed);
	nim->SetLineStyle(2);
	nim->SetLineWidth(2);

	
	double ymin = 0.064; // match to Blake's plot
	double lab_x = 0.32;  // for resolution label
	double lab_y = 0.072; //  "          "
 	if(fun2->Eval(0.3)<0.06){
		ymin = 0.0;
		lab_x = 0.37;
		lab_y = 0.03;
	}

	TH2D *axes = new TH2D("axes","", 100, 0.120, 0.685, 100, ymin, 0.145);
	axes->SetStats(0);
	axes->SetYTitle("#sigma_{E}/E");
	axes->SetXTitle("Generated Energy (GeV)");
	axes->Draw();

	gres->Draw("Psame");
	fun2->Draw("same");
	nim->Draw("same");
	
	char str[256];
	sprintf(str, "#frac{#sigma_{E}}{E} = #frac{(%3.1f#pm%3.1f)%%}{#sqrt{E}} #oplus (%3.1f#pm%3.1f)%% #oplus #frac{(%3.1f#pm%3.1f)%%}{E} #oplus (%3.1f#pm%3.1f)%%#timesln(E)"
		, fabs(fun->GetParameter(0))*100.0, fun->GetParError(0)*100.0
		, fabs(fun->GetParameter(1))*100.0, fun->GetParError(1)*100.0
		, fabs(fun->GetParameter(2))*100.0, fun->GetParError(2)*100.0
		, fabs(fun->GetParameter(3))*100.0, fun->GetParError(3)*100.0);
	TLatex *lab = new TLatex(0.13, lab_y, str);
	lab->SetTextSize(0.035);
	lab->SetTextAlign(12);
	lab->Draw();
	
	StandardLabels(axes, Scheme(), "", AngleStr("#theta_{#gamma}="));
	
	c1->SaveAs("Eres_vs_E.pdf");
	c1->SaveAs("Eres_vs_E.png");

	// Write parameter to a file so they can be read in and used
	// to automatically generate the LaTeX table later
	ofstream ofs("Eres_parms.dat");
	ofs<<Scheme()<<endl;
	ofs<<Angle()<<endl;
	ofs<<fun->GetParameter(0)*100.0<<" "<<fun->GetParError(0)*100.0<<endl;
	ofs<<fun->GetParameter(1)*100.0<<" "<<fun->GetParError(1)*100.0<<endl;
	ofs<<fun->GetParameter(2)*100.0<<" "<<fun->GetParError(2)*100.0<<endl;
	ofs<<fun->GetParameter(3)*100.0<<" "<<fun->GetParError(3)*100.0<<endl;
	ofs.close();
}

//------------------
// 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);
	
	// Define an asymmetric Gaussian function
	TF1 *f = new TF1("agaus", asymmetric_gaus, 0.0, 2.0, 4); 

	// 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){

#if 0		
			// Fit to Gaussian
			int status = h_py->Fit("gaus","");
			if(status!=0){
				start_bin++;
				continue;
			}
			TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
			c1->Update();
			//c1->SaveAs("symmetric_gauss_fit.png");
			//if(start_bin>30)return;

			// Get fit results
			double mean = h_py->GetFunction("gaus")->GetParameter(1);
			double sigma = h_py->GetFunction("gaus")->GetParameter(2);
			double mean_err = h_py->GetFunction("gaus")->GetParError(1);
			double sigma_err = h_py->GetFunction("gaus")->GetParError(2);
#else
			// Fit to Asymmetric Gaussian
			f->SetParameter(0, h_py->GetMaximum());
			f->SetParameter(1, h_py->GetBinCenter(h_py->GetMaximumBin()));
			f->SetParameter(2, 0.03);
			f->SetParameter(3, 0.006);
			//h_py->Fit(f,"0Q");
			Int_t status = h_py->Fit(f,"Q");
			double sigma1 = fabs(h_py->GetFunction("agaus")->GetParameter(2));
			double sigma2 = fabs(h_py->GetFunction("agaus")->GetParameter(3));
			if(status==0 && (sigma1>1.0 || sigma2>1.0))status = -1;
			if(status != 0){
				// Fit failed
				cout<<"BAD FIT! Skipping start_bin and trying again ...."<<endl;
				//end_bin = start_bin;
				start_bin++;
				continue;
			}
			// Fit succeeded
			
			char hname[256];
			sprintf(hname, "asym%02d", Npoints);
			h_py->Clone(hname);

			TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
			if(c1)c1->Update();
			//c1->SaveAs("asymmetric_gauss_fit.png");
			//sleep(1);
			//if(start_bin>30)return;

			// Get fit results
			double mean = h_py->GetFunction("agaus")->GetParameter(1);
			double mean_err = h_py->GetFunction("agaus")->GetParError(1);

			// Variance of asymmetric Gaussian function about center
			// of individual Gaussians is given by
			// 
			//
			//                 (s1^3 + s2^3)
			//   sigma_tot^2 = ----------------------------
			//                 (s1   +  s2)
			//
			//
			//  where s1 = sigma_1  and s2 = sigma_2
			//
			//
			// The uncertainty on the variance is derived by taking
			// derivatives of the above. It ends up being:
			//
			//
			//                   3[(s1^2)(s1_err) + (s2^2)(s2_err)] - (sigma_tot^2)[s1_err+s2_err]
			//  sigma_tot_err = --------------------------------------------------------------------
			//                                        2sigma_tot(s1 + s2)
			//
			//
			//
			// The variance about the mean of the asymmetric distribution
			// is given by
			//
			//  sigma_mu^2 = sigma_tot^2 - mu^2
			//
			//
			//  The mean "mu" is given by
			//
			//                        (s2^2 - s1^2)
			//     mu = sqrt(2/pi)* -----------------
			//                          (s1 + s2)
			//
			//
			cout<<"sigma1="<<sigma1<<"  sigma2="<<sigma2<<endl;
			double mu = sqrt(2.0/TMath::Pi())*(sigma2*sigma2 - sigma1*sigma1)/(sigma2 + sigma1);
				    //mu -= mean;
			double sigma = sqrt((pow(sigma1, 3.0) + pow(sigma2, 3.0))/(sigma1 + sigma2)); // (verified via macro that this is correct)
				    sigma = sqrt(sigma*sigma - mu*mu);
			double sigma1_err = h_py->GetFunction("agaus")->GetParError(2);
			double sigma2_err = h_py->GetFunction("agaus")->GetParError(2);
			double sigma_err = 3.0*(pow(sigma1,2.0)*sigma1_err + pow(sigma2,2.0)*sigma2_err) - pow(sigma,2.0)*(sigma1_err+sigma2_err);
				    sigma_err /= 2.0*sigma*(sigma1 + sigma2);
#endif
			// 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(kMagenta);
	gres->SetLineColor(kMagenta);
	gres->SetLineWidth(2);

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

//-----------------
// asymmetric_gaus
//-----------------
Double_t asymmetric_gaus(Double_t *xptr, Double_t *par)
{
	double x = xptr[0];
	double amp = par[0];
	double mean = par[1];
	double sigma = x<mean ? par[2]:par[3];

	return amp*exp(-0.5*pow((x-mean)/sigma, 2.0));
}