#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TGaxis.h"
#include "TRandom.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TGraph.h"
#include "TSystem.h"
#include "TTree.h"
#include "TTreePlayer.h"
#include "TF1.h"
#include "TCut.h"
#include "TPaletteAxis.h"
#include "TFitResult.h"
#include "TMath.h"

double errort=1;
double tmin=-5000;
double tpeakmin=-2500;
double tpeakmax=2500;
double tmax=5000;
double errorTotArea, errorBkgAreaSx, errorBkgAreaDx, errorPeakArea;



void DrawVerticalLine(Double_t x)
{
   TLine l;
   l.SetLineColor(kRed);
   l.SetLineStyle(10);
   l.SetLineWidth(1.);
   Double_t lm = gPad->GetLeftMargin();
   Double_t rm = 1.-gPad->GetRightMargin();
   Double_t tm = 1.-gPad->GetTopMargin();
   Double_t bm = gPad->GetBottomMargin();
   Double_t xndc = (rm-lm)*((x-gPad->GetUxmin())/(gPad->GetUxmax()-gPad->GetUxmin()))+lm;
   l.DrawLineNDC(xndc,bm,xndc,tm);
}

void DrawHorizontalLine(Double_t x)
{
   TLine l;
   l.SetLineColor(kRed);
   l.SetLineStyle(10);
   l.SetLineWidth(1.);
   Double_t lm = gPad->GetLeftMargin();
   Double_t rm = 1.-gPad->GetRightMargin();
   Double_t tm = 1.-gPad->GetTopMargin();
   Double_t bm = gPad->GetBottomMargin();
   //Double_t yndc = (rm-lm)*((x-gPad->GetUymin())/(gPad->GetUymax()-gPad->GetUymin()))+lm;
   Double_t yndc = (tm-bm)*((x-gPad->GetUymin())/(gPad->GetUymax()-gPad->GetUymin()))+bm;
   l.DrawLineNDC(lm,yndc,rm,yndc);
}


void deltatbkgest()
{

	float fillcolor=0.35;
	gStyle->SetFitFormat("8.6g");
	gStyle->SetOptStat(111110);
	ofstream results;
	results.open("ResultsBackgroundRun0610.txt",  ios::out);
	TFile *fin = TFile::Open("run0610.root");
	 	fin->ls ();
		TGaxis::SetMaxDigits(3);
 		if (fin == 0) {
	    	printf("Error: cannot open the file!\n");
	   	}
   		else {
			
			TTree *te=0;
			fin->GetObject("tree",te);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			te->SetLineColor(kBlue);
			TH1::StatOverflows(true);
			TCanvas *c01 = new TCanvas("c01","c01",1280,1024);
			c01->SetLogy();
			TString string = TString::Format("t1-t0 >> htemp(1000, -5000, 5000)");
			TCut coincidence = TString::Format("abs(t1-t0)<5000").Data();
			TCut pulser0 = TString::Format("e2==0").Data();
		    TCut pulser1 = TString::Format("e3==0").Data();
		    TCut GePD = TString::Format("e0>0").Data();
		    TCut GeDD = TString::Format("e1>0").Data();
			te->Draw(string,coincidence && GePD && GeDD && pulser0 && pulser1);
 			te->GetHistogram()->SetTitle("");
 			te->SetScanField(0);
  	 		TH1F *hgraph= (TH1F*)gPad->GetPrimitive("htemp"); 
  			hgraph->GetXaxis()->SetTitle("#DeltaT (ns)"); 
			hgraph->GetYaxis()->SetTitle("Counts");
			hgraph->GetYaxis()->SetTitleSize(40);
  			hgraph->GetYaxis()->SetTitleFont(43);
   			hgraph->GetYaxis()->SetTitleOffset(1.1);
   			hgraph->GetYaxis()->SetLabelFont(43); 
   			hgraph->GetYaxis()->SetLabelSize(40);
   			hgraph->GetXaxis()->SetTitleSize(40);
   			hgraph->GetXaxis()->SetTitleFont(43);
   			hgraph->GetXaxis()->SetTitleOffset(1.1);
   			hgraph->GetXaxis()->SetLabelFont(43); 
   			hgraph->GetXaxis()->SetLabelSize(40);
   			hgraph->SetFillStyle(3005);
 			hgraph->SetFillColorAlpha(kBlue, fillcolor);
   			hgraph->SetStats(0);
	   		hgraph->Draw();
	   		hgraph->SetName("");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c01->Update();
			gPad->Modified();
   			/*TLegend* leg = new TLegend(0.25, 0.7, .35, .75);
   			leg->SetHeader("Legend");
  			leg->SetNColumns(1);
   			leg->AddEntry(hgraph, "MC", "l");
   			leg->Draw();*/ 
   			c01->Update();
   			gPad->Modified();
			gPad->Update();
			/*TPaveStats *stats = (TPaveStats*)hgraph->GetListOfFunctions()->FindObject("stats");
			stats->SetTextColor(kBlue);
			stats->SetX1NDC(0.80); stats->SetX2NDC(0.98);
			stats->SetY1NDC(0.83); stats->SetY2NDC(0.98);
			//stats->AddText(TString::Format("Text");
			stats->DrawClone();*/
			gPad->Update();
			DrawVerticalLine(-2500.);
   			DrawVerticalLine(2500.);
   			DrawHorizontalLine(10);
   			TBox *box1 = new TBox(0.5,0.5,0.8,0.8);
   			box1->SetFillColor(19);
   			box1->Draw();
   			TPaveText *t=new TPaveText(0.42,0.15,0.57,0.3,"brNDC");
 			t->AddText("Random coincidences");
 			t->SetLineColor(0);
 			t->SetFillStyle(0);
  			t->Draw();
   			c01->Update();
   			gPad->Modified();
   			TPaveText *t2=new TPaveText(0.42,0.45,0.57,0.6,"brNDC");
 			t2->AddText("True coincidences");
 			t2->SetLineColor(0);
 			t2->SetFillStyle(0);
  			t2->Draw();
  			c01->Update();
   			gPad->Modified();
			c01->Print("DeltaT[-5us,5us]run0610delay_bkgestimation.png");
			int BinMin = hgraph->GetXaxis()-> FindBin(tmin);
   			int BinMax = hgraph->GetXaxis()-> FindBin(tmax);
			int BinPeakMin = hgraph->GetXaxis()-> FindBin(tpeakmin);
   			int BinPeakMax = hgraph->GetXaxis()-> FindBin(tpeakmax);
			double TotArea= hgraph->IntegralAndError(BinMin,BinMax,errorTotArea,""); 
			double BkgAreaSx= hgraph->IntegralAndError(BinMin,BinPeakMin-1,errorBkgAreaSx,""); 
			double BkgAreaDx= hgraph->IntegralAndError(BinPeakMax+1,BinMax,errorBkgAreaDx,""); 
			double BkgArea= BkgAreaSx+BkgAreaDx;
			double errorBkgArea=(errorBkgAreaSx*errorBkgAreaSx)+(errorBkgAreaDx*errorBkgAreaDx);
			errorBkgArea=TMath::Sqrt(errorBkgArea);
			double PeakArea= hgraph->IntegralAndError(BinPeakMin,BinPeakMax,errorPeakArea,""); 
			double texternalSx=abs(tpeakmin-tmin);
			double errortexternalSx= ((((tpeakmin-tmin)/abs(tpeakmin-tmin))*errort)*(((tpeakmin-tmin)/abs(tpeakmin-tmin))*errort))+((((tmin-tpeakmin)/abs(tpeakmin-tmin))*errort)*(((tmin-tpeakmin)/abs(tpeakmin-tmin))*errort));
			errortexternalSx=TMath::Sqrt(errortexternalSx);
			double texternalDx=abs(tmax-tpeakmax);
			double errortexternalDx=((((tmax-tpeakmax)/abs(tmax-tpeakmax))*errort)*(((tmax-tpeakmax)/abs(tmax-tpeakmax))*errort))+((((tpeakmax-tmax)/abs(tmax-tpeakmax))*errort)*(((tpeakmax-tmax)/abs(tmax-tpeakmax))*errort));
			errortexternalDx= TMath::Sqrt(errortexternalDx);
			double texternal = texternalDx + texternalSx;
			double errortexternal = (errortexternalSx*errortexternalSx)+(errortexternalDx*errortexternalDx);
			errortexternal = TMath::Sqrt(errortexternal);
			double tpeak=abs(tpeakmax-tpeakmin);
			double errortpeak= ((((tpeakmax-tpeakmin)/abs(tpeakmax-tpeakmin))*errort)*((tpeakmax-tpeakmin)/abs(tpeakmax-tpeakmin))*errort)+ ((((tpeakmin-tpeakmax)/abs(tpeakmax-tpeakmin))*errort)*((tpeakmin-tpeakmax)/abs(tpeakmax-tpeakmin))*errort);
			errortpeak = TMath::Sqrt(errortpeak);
			double PeakBkgArea= (BkgArea/texternal)*tpeak;
			double errorPeakBkgArea= (((tpeak/texternal)*errorBkgArea)*((tpeak/texternal)*errorBkgArea))+(((BkgArea*tpeak)/((texternal)*(texternal)))*errortexternal)*((((BkgArea*tpeak)/((texternal)*(texternal)))*errortexternal))+((BkgArea/texternal*errortpeak)*(BkgArea/texternal*errortpeak));
			errorPeakBkgArea= TMath::Sqrt(errorPeakBkgArea);
			double PeakBkgPercentage= 100*PeakBkgArea/PeakArea;
			double errorPeakBkgPercentage= ((errorPeakBkgArea/PeakBkgArea)*(errorPeakBkgArea/PeakBkgArea))+((errorPeakArea/PeakArea)*(errorPeakArea/PeakArea));
			errorPeakBkgPercentage=100*TMath::Sqrt(errorPeakBkgPercentage);
			double TrueEventsArea= PeakArea-PeakBkgArea;
			double errorTrueEventsArea= (errorPeakArea*errorPeakArea)+(errorPeakBkgArea*errorPeakBkgArea);
			errorTrueEventsArea=TMath::Sqrt(errorTrueEventsArea);
			results << "****************Background Estimation run0610 ****************"<<  endl;
    		results << "Time min= " << tmin << " +/- " << errort <<   endl;
    		results << "Time peak min= " << tpeakmin << " +/- " << errort  <<  endl;
    		results << "Time peak max= " << tpeakmax << " +/- " << errort  <<  endl;
    		results << "Time max= " << tmax << " +/- " << errort <<   endl;
			results << "Time background external Sx= " << texternalSx << " +/- " << errortexternalSx <<  endl;
			results << "Time background external Dx= " << texternalDx << " +/- " << errortexternalDx <<  endl;
			results << "Time background external = " << texternal << " +/- " << errortexternal <<  endl;
    		results << "Time peak= " << tpeak << " +/- " << errortpeak <<  endl;
			results << "Total Area= " << TotArea << " +/- " << errorTotArea <<  endl;
    		results << "Background Area Sx= " << BkgAreaSx << " +/- " << errorBkgAreaSx <<  endl;
    		results << "Background Area Dx= " << BkgAreaDx << " +/- " << errorBkgAreaDx  << endl;
    		results << "Background Area= " << BkgArea << " +/- " << errorBkgArea  << endl;
    		results << "Peak Area= " << PeakArea << " +/- " << errorPeakArea  << endl;
    		results << "Peak Background Area= " << PeakBkgArea  << " +/- " << errorPeakBkgArea  << endl;
    		results << "Peak Background Area (%)= 100*PeakBkgArea/PeakArea=  " << PeakBkgPercentage << " +/- " << errorPeakBkgPercentage  << endl;
    		results << "True Events area= " << TrueEventsArea << " +/- " << errorTrueEventsArea  << endl;
			delete c01;	
	}
}