#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 "TPad.h"
#include "TCut.h"
#include "TFractionFitter.h"


float offx=1.2;
float offy=1.0;
float margr=0.08;
float w=3;
float margl=0.12;
double Emin=1200;
double Emax=5000;
int BinSumMin0;
int BinSumMax0;

void mcdata()
{
	gStyle->SetOptStat(00000);
	gStyle->SetOptFit(0);
	TH1::StatOverflows(true);
	TFile *fRea = TFile::Open("run0611.root");
    if (fRea == 0) {
      // if we cannot open the file, print an error message and return immediatly
    	printf("Error: cannot open Rea.root!\n");
      	return;
    }
    fRea->ls();
    TFile *fBIB = TFile::Open("run0735.root");
    if (fBIB == 0) {
      // if we cannot open the file, print an error message and return immediatly
    	printf("Error: cannot open BIB.root!\n");
      	return;
    }
   	fBIB->ls();
   	TFile *fMC97 = TFile::Open("SimLuna.1.root");
    if (fMC97 == 0) {
      // if we cannot open the file, print an error message and return immediatly
    	printf("Error: cannot open MC 6997 .root!\n");
      	return;
    }
    fMC97->ls();
   TFile *fMC98 = TFile::Open("SimLuna98.1.root");
    if (fMC98 == 0) {
      // if we cannot open the file, print an error message and return immediatly
    	printf("Error: cannot open MC 6998.root!\n");
      	return;
    }
   	fMC98->ls();
	TTree *tRea=0;
   	TTree *tBIB=0;
   	TTree *tMC97=0;
   	TTree *tMC98=0;
   	fRea->GetObject("tree",tRea);
   	fBIB->GetObject("tree",tBIB);
   	fMC97->GetObject("Tree1",tMC97);
   	fMC98->GetObject("Tree1",tMC98);
	if (tRea == nullptr ) {
      printf("Error: cannot get Tree Rea!\n");
    	return;
   	}
   	if (tBIB == nullptr ) {
      printf("Error: cannot get Tree BIB!\n");
    	return;
   	}
   	if (tMC97 == nullptr ) {
      printf("Error: cannot get Tree MC 970!\n");
    	return;
   	}
   	if (tMC98 == nullptr ) {
      printf("Error: cannot get Tree MC 980!\n");
    	return;
   	}
   	tRea->SetLineColor(kBlue);
	tBIB->SetLineColor(kRed);
	tMC97->SetLineColor(kMagenta);
	tMC98->SetLineColor(kGreen);
    TCanvas *c01 = new TCanvas("c01","c01",1280,1024);
	TPad *pad0a = new TPad("pad0a", "pad0a", 0.0, 0.3, 1.0, 1.0);
	pad0a->Draw();             // Draw the upper pad: pad0a
   	pad0a->cd();
   	pad0a->SetLogy();
    gPad->SetLeftMargin(margl);
   	gPad->SetRightMargin(margr);
	gROOT->cd(); // newly created histograms should go here
	TString hStringRea0 = TString::Format("e0 >> htempRea0(8000, 0., 8000.)"); 
	TString hStringBIB0 = TString::Format("e0 >> htempBIB0(8000, 0., 8000.)"); 
	TString hStringMC970 = TString::Format("Edep[0] >> htempMC970(8000, 0., 8000.)"); 
	TString hStringMC980 = TString::Format("Edep[0] >> htempMC980(8000, 0., 8000.)"); 
	TCut cut0a = TString::Format("e0>0").Data();
	TCut cut0b = TString::Format("e2==0").Data();
	TCut cut0c = TString::Format("Edep[0]>0").Data();
	tRea->Draw(hStringRea0,cut0a && cut0b,"HIST");
	TH1F *hRea0 = (TH1F*)gPad->GetPrimitive("htempRea0");
	hRea0->SetOption("HIST"); 
	hRea0->Scale(1./28.51842226);
	c01->Modified(); 
	c01->Update(); 
   	hRea0->SetTitle("; Energy (keV); Counts");
	c01->Modified(); 
	c01->Update(); 
   	hRea0->SetFillColor(kBlue);
   	hRea0->SetFillStyle(3005);
   	hRea0->SetMaximum(1.e4);
   	hRea0->SetMinimum(1.e-5);
   	tBIB->Draw(hStringBIB0,cut0a && cut0b,"HIST SAME");	
	TH1F *hBIB0 = (TH1F*)gPad->GetPrimitive("htempBIB0");
	hBIB0->SetOption("HIST"); 
	hBIB0->Scale(1./19.28593866);
    c01->Modified(); 
	c01->Update();  
    hBIB0->SetFillColor(kRed);
    hBIB0->SetFillStyle(3005);
    hBIB0->SetTitle("; Energy (keV); Counts");
    hBIB0->GetXaxis()->SetTitleOffset(offx);
    hBIB0->GetYaxis()->SetTitleOffset(offy);
    hBIB0->SetStats(0);
    hBIB0->SetMaximum(1.e4);
    hBIB0->SetMinimum(1.e-5);
   	c01->Update();
   	tMC97->Draw(hStringMC970,cut0c, "HIST SAME");	
	TH1F *hMC970 = (TH1F*)gPad->GetPrimitive("htempMC970");
	hMC970->SetOption("HIST"); 
	//hMC970->Scale(1.);
    c01->Modified(); 
	c01->Update();  
    hMC970->SetFillColor(kMagenta);
    hMC970->SetFillStyle(3005);
    hMC970->SetTitle("; Energy (keV); Counts");
    hMC970->GetXaxis()->SetTitleOffset(offx);
    hMC970->GetYaxis()->SetTitleOffset(offy);
    hMC970->SetStats(0);
    hMC970->SetMaximum(1.e4);
    hMC970->SetMinimum(1.e-5);
   	c01->Update();
   	tMC98->Draw(hStringMC980,cut0c, "HIST SAME");	
	TH1F *hMC980 = (TH1F*)gPad->GetPrimitive("htempMC980");
	hMC980->SetOption("HIST"); 
	//hMC980->Scale(1.);
    c01->Modified(); 
	c01->Update();  
    hMC980->SetFillColor(kGreen);
    hMC980->SetFillStyle(3005);
    hMC980->SetTitle("; Energy (keV); Counts");
    hMC980->GetXaxis()->SetTitleOffset(offx);
    hMC980->GetYaxis()->SetTitleOffset(offy);
    hMC980->SetStats(0);
    hMC980->SetMaximum(1.e4);
    hMC980->SetMinimum(1.e-5);
   	c01->Update();
	TH1F *hSum0 = (TH1F *)hMC970->Clone("MCSum0");
    if (!(hSum0->GetSumw2N() > 0))hSum0->Sumw2(kTRUE);
	hSum0->Add(hMC980,+1);
	hSum0->Add(hBIB0,+1);
	//hMC980->Scale(1.);
	hSum0-> SetLineColor(kOrange);
	hSum0->SetFillStyle(3005);
	hSum0->GetXaxis()->SetRangeUser(0,8000);
	hSum0->SetFillColor(kOrange);
	hSum0->SetTitle("; Energy (keV); Difference");
	hSum0->SetOption("HIST"); 
//	BinSumMin0 = hSum0->GetXaxis()-> FindBin(Emin);
  // 	BinSumMax0 = hSum0->GetXaxis()-> FindBin(Emax);
   //	hSum0->Scale(1./hSum0->Integral(BinSumMin0,BinSumMax0), "width");
	hSum0->Draw("HIST SAME");
   	c01->Update();
    TLegend* leg0 = new TLegend(0.35, 0.7, .5, .8);
   	leg0->SetHeader("");
  	leg0->SetNColumns(1);
   	leg0->AddEntry(hRea0, "Rea/Q_{Rea}: Run0611");
   	leg0->AddEntry(hBIB0, "BIB/Q_{BIB}: Run0735");
   	leg0->AddEntry(hMC970, "MC 6997");
   	leg0->AddEntry(hMC980, "MC 6998");
   	leg0->AddEntry(hSum0, "MC 6997 + MC 6998 +  BIB/Q_{BIB}");
   	leg0->SetBorderSize(0);
   	leg0->Draw(); //if wanted
    c01->cd();
    c01->Update();
    TPad *pad0b = new TPad("pad0b", "pad0b", 0.0, 0.0, 1.0, 0.3);
    pad0b->SetTopMargin(0.0);
    pad0b->SetGridx(); // vertical grid
    pad0b->SetLogy();
    pad0b->Draw();
    pad0b->cd();       // pad0b becomes the current pad
    TH1F *diff0 = (TH1F *)hRea0->Clone("diff0");
    if (!(diff0->GetSumw2N() > 0))diff0->Sumw2(kTRUE);
	diff0->Add(hSum0,-1);
	diff0-> SetLineColor(kGreen);
	diff0->GetXaxis()->SetRangeUser(0,8000);
	diff0->SetFillColor(kGreen);
	diff0->SetTitle("; Energy (keV); Difference");
	diff0->SetOption("HIST"); 
	diff0->Draw("EP");
//	diff0->Draw();
	TLegend* leg0r = new TLegend(0.45, 0.8, .6, .9);
   	leg0r->SetHeader("");
    leg0r->SetNColumns(1);
   	leg0r->AddEntry(diff0, "Rea/Q_{Rea} - (MC 6997 + MC 6998 + BIB/Q_{BIB})");
   	leg0r->SetBorderSize(0);
   	leg0r->Draw(); //if wanted
    gPad->SetLeftMargin(margl);
    gPad->SetRightMargin(margr);
    gPad->SetBottomMargin(0.15);
	c01->Print("Simulazione/Rea0611-MC6998-MC6997-BIB0735_GePD.png"); 
	
	TCanvas *c02 = new TCanvas("c01","c01",1280,1024);
	TPad *pad1a = new TPad("pad1a", "pad1a", 0.0, 0.3, 1.0, 1.0);
	pad1a->Draw();             // Draw the upper pad: pad0a
   	pad1a->cd();
   	pad1a->SetLogy();
    gPad->SetLeftMargin(margl);
   	gPad->SetRightMargin(margr);
   	hRea0->Draw("HIST");
   	hSum0->Draw("HIST SAME");
   	TLegend* leg1 = new TLegend(0.35, 0.7, .5, .8);
   	leg1->SetHeader("");
  	leg1->SetNColumns(1);
   	leg1->AddEntry(hRea0, "Rea/Q_{Rea}: Run0611");
   	leg1->AddEntry(hSum0, "MC 6997 + MC 6998 +  BIB/Q_{BIB}");
   	leg1->SetBorderSize(0);
   	leg1->Draw(); //if wanted
   	c02->cd();
    c02->Update();
   	TPad *pad1b = new TPad("pad1b", "pad1b", 0.0, 0.0, 1.0, 0.3);
    pad1b->SetTopMargin(0.0);
    pad1b->SetGridx(); // vertical grid
    pad1b->SetLogy();
    pad1b->Draw();
    pad1b->cd();       // pad1b becomes the current pad
    diff0->Draw("EP");
    TLegend* leg1r = new TLegend(0.45, 0.8, .6, .9);
   	leg1r->SetHeader("");
    leg1r->SetNColumns(1);
   	leg1r->AddEntry(diff0, "Rea/Q_{Rea} - (MC 6997 + MC 6998 + BIB/Q_{BIB})");
   	leg1r->SetBorderSize(0);
   	leg1r->Draw(); //if wanted
    gPad->SetLeftMargin(margl);
    gPad->SetRightMargin(margr);
    gPad->SetBottomMargin(0.15);
    c02->Print("Simulazione/Rea0611-MC6998-MC6997-BIB0735-JustSum_GePD.png"); 
    
    TCanvas *c03 = new TCanvas("c01","c01",1280,1024);
    gPad->SetLeftMargin(margl);
   	gPad->SetRightMargin(margr);
   	 TObjArray *mcbib = new TObjArray(3);        // MC and BIB histograms are put in this array
   	 mcbib->Add(hMC970);
   	 mcbib->Add(hMC980);
   	 mcbib->Add(hBIB0);
   	 TFractionFitter* fit = new TFractionFitter(hRea0, mcbib); // initialise
   fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
   fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
   Int_t status = fit->Fit();               // perform the fit
   std::cout << "fit status: " << status << std::endl;
   if (status == 0) {                       // check on fit status
     TH1F* result = (TH1F*) fit->GetPlot();
     hRea0->Draw("Ep");
     result->Draw("same");
 }
     c03->Print("Simulazione/Rea0611-MC6998-MC6997-BIB0735-Fitter_GePD.png"); 
    
    delete c01;
 	delete c02;
 	delete c03;
 
}