#include <TROOT.h>
#include <TTree.h>
#include <TH2.h>
#include <TGraphErrors.h>
#include <TFile.h>
#include <TF1.h>
#include <TCanvas.h>

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


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

#include "StandardLabels.C"


//-----------------
// enoise_compare
//-----------------
void enoise_compare(void)
{

	TFile *f = new TFile("electronic_noise.root");
	f->cd();
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	// Small signal
	TH1D *before = (TH1D*)gROOT->FindObject("elec_up_before_lay2_sec1");
	TH1D *after  = (TH1D*)gROOT->FindObject("elec_up_after_lay2_sec1");
	
	TH2D *axes = new TH2D("axes","", 100, 0.0, 130.0, 100, -5.0, 1.1*before->GetMaximum());
	axes->SetStats(0);
	axes->SetXTitle("Time (ns)");
	axes->SetYTitle("Amplitude (mV)");
	axes->Draw();
	
	after->SetLineColor(kRed);
	after->SetLineWidth(3);
	before->SetLineWidth(2);

	after->Draw("same");
	before->Draw("same");
	before->Draw("same");
	
	TLine *lin = new TLine(0.0, 44.7, 130.0, 44.7);
	lin->SetLineColor(kMagenta);
	lin->SetLineWidth(2);
	lin->SetLineStyle(2);
	lin->Draw();
	
	TLatex *lab = new TLatex(20.0, 45.0, "44.7mV");
	lab->SetTextColor(kMagenta);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(21);
	lab->Draw();
	
	StandardLabels(axes, AngleStr("#theta_{#gamma}=", Scheme()), "Upstream - small signal", "electronic noise: 1GHz   1mV   2mV/ns");

	c1->SaveAs("enoise_compare_small.pdf");
	c1->SaveAs("enoise_compare_small.png");

	// Large signal
	TH2D *axes = new TH2D("axes","", 100, 0.0, 130.0, 100, -5.0, 3500.0);
	axes->SetStats(0);
	axes->SetXTitle("Time (ns)");
	axes->SetYTitle("Amplitude (mV)");
	axes->Draw();
	
	TH1D *before = (TH1D*)gROOT->FindObject("elec_dn_before_lay2_sec2");
	TH1D *after  = (TH1D*)gROOT->FindObject("elec_dn_after_lay2_sec2");
	
	after->SetLineColor(kRed);
	after->SetLineWidth(3);
	before->SetLineWidth(2);

	after->Draw("same");
	before->Draw("same");
	
	TLine *lin = new TLine(0.0, 44.7, 130.0, 44.7);
	lin->SetLineColor(kMagenta);
	lin->SetLineWidth(2);
	lin->SetLineStyle(2);
	lin->Draw();
	
	TLatex *lab = new TLatex(40.0, 60.0, "44.7mV");
	lab->SetTextColor(kMagenta);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(21);
	lab->Draw();
	
	StandardLabels(axes, AngleStr("#theta_{#gamma}=", Scheme()), "Downstream - large signal", "electronic noise: 1GHz   1mV   2mV/ns");

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


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

	// 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> yerr;
	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->GetFunction("gaus")->GetParameter(1);
			double sigma = h_py->GetFunction("gaus")->GetParameter(2);
			double sigma_err = h_py->GetFunction("gaus")->GetParError(2);

			// 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.050; // ns
			sigma_err = sqrt(sigma_err*sigma_err + epsilon*epsilon);

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

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

	gmean = new TGraph(Npoints, &x[0], &y_mean[0]);
	gmean->SetMarkerColor(kRed);
	gmean->SetMarkerStyle(22);
	gmean->SetLineColor(kRed);
	
	// Fit the graph. If there are less than 8 points, just fit to 
	// a flat line.
	TF1 *fun;
	if(Npoints<4){
		fun  = new TF1("fun", "[0]");
	}else{
		fun = new TF1("fun", "[0] + [1]/(pow(x,[2])+[3])");
		fun->SetParameter(0, 0.037);
		fun->SetParameter(1, 1000.0);
		fun->SetParameter(2, 1.0);
		fun->SetParameter(3, -3.0);
		
		fun->SetParLimits(2, 0.0, 10.0);
	}
	
	fun->SetLineColor(kBlue);
	gres->Fit(fun);

}