#include <iostream>
#include <fstream>
#include <vector>
#include <map>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TRandom3.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TParticle.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TStyle.h"
#include "TNtuple.h"
#include "TLatex.h"
#include "TGaxis.h"
#include "TGraphAsymmErrors.h"
#include "TEventList.h"
#include "TLegend.h"

std::map<int, TString> sel45;
std::map<int, TString> sel55;
TString channels[8];

// *** setup  style
void setStyle()
{
	gStyle->SetPadTopMargin(0.13);
	gStyle->SetPadBottomMargin(0.25);
	gStyle->SetPadLeftMargin(0.15);
	gStyle->SetPadRightMargin(0.02);

	gStyle->SetLabelSize(0.062,"X");
	gStyle->SetLabelSize(0.052,"Y");

/*	gStyle->SetStatY(0.87);
	gStyle->SetStatX(0.986);
	gStyle->SetStatW(0.4);
	gStyle->SetStatH(0.12);
	gStyle->SetOptStat(0); // only entries*/
	
	gStyle->SetTitleH(0.12);
	gStyle->SetTitleX(0.1);
	gStyle->SetTitleW(0.9);
	
	gStyle->SetTitleOffset(1.1,"x");
	gStyle->SetTitleXSize(0.1);
	
	gStyle->SetNdivisions(505);
}

void config_histo(TH1 *h)
{
	h->SetYTitle("t [d]");
	h->GetXaxis()->LabelsOption("v");
	h->SetStats(0);
	h->SetMarkerSize(1.5);
	h->SetMarkerStyle(20);
	h->SetMarkerColor(4);
}

void config_graph(TGraph *h, TString name, int col=4)
{
	h->SetName(name);
	h->SetMarkerSize(1.5);
	h->SetMarkerStyle(20);
	h->SetMarkerColor(col);
}

void init_channels()
{
	channels[0]  = "#eta_{c}(KK#pi^{0}) 2#pi";
	channels[1]  = "#eta_{c}(KK#pi^{0}) 2#pi^{0}";
	channels[2]  = "#eta_{c}(K_{S}K#pi) 2#pi";
	channels[3]  = "#eta_{c}(K_{S}K#pi) 2#pi^{0}";
	
	channels[4]  = "J/#psi(e^{+}e^{-}) 2#pi";
	channels[5]  = "J/#psi(e^{+}e^{-}) 2#pi^{0}";
	channels[6]  = "J/#psi(#mu^{+}#mu^{-}) 2#pi";
	channels[7]  = "J/#psi(#mu^{+}#mu^{-}) 2#pi^{0}";
	
/*	channels[0]  = "#eta_{c}(1) 2#pi";
	channels[1]  = "#eta_{c}(1) 2#pi^{0}";
	channels[2]  = "#eta_{c}(2) 2#pi";
	channels[3]  = "#eta_{c}(2) 2#pi^{0}";
	
	channels[4]  = "J/#psi(1) 2#pi";
	channels[5]  = "J/#psi(1) 2#pi^{0}";
	channels[6]  = "J/#psi(2) 2#pi";
	channels[7]  = "J/#psi(2) 2#pi^{0}";*/
	
/*	channels[0]  = "#eta_{c} #pi^{+} #pi^{-}; #eta_{c} -> K^{+} K^{-} #pi^{0}";
	channels[1]  = "#eta_{c} #pi^{0} #pi^{0}; #eta_{c} -> K^{+} K^{-} #pi^{0}";
	channels[2]  = "#eta_{c} #pi^{+} #pi^{-}; #eta_{c} -> K_{S} K^{#pm} #pi^{#mp}";
	channels[3]  = "#eta_{c} #pi^{0} #pi^{0}; #eta_{c} -> K_{S} K^{#pm} #pi^{#mp}";
	
	channels[4]  = "J/#psi #pi^{+} #pi^{-}; J/#psi -> e^{+} e^{-}";
	channels[5]  = "J/#psi #pi^{0} #pi^{0}; J/#psi -> e^{+} e^{-}";
	channels[6]  = "J/#psi #pi^{+} #pi^{-}; J/#psi -> #mu^{+} #mu^{-}";
	channels[7]  = "J/#psi #pi^{0} #pi^{0}; J/#psi -> #mu^{+} #mu^{-}";*/
	
}

void init_cuts()
{
	sel45[0]   = "";
	sel45[1]   = "";	
	sel45[10]  = "";
	sel45[11]  = "";	
	sel45[20]  = "abs(fpbarp_d0m-3.096)<0.036&&pbarp_d0d0pide>0.5&&pbarp_d0d1pide>0.5";
	sel45[21]  = "abs(fpbarp_d0m-3.096)<0.036&&pbarp_d0d0pide>0.5&&pbarp_d0d1pide>0.5";	
	sel45[30]  = "abs(fpbarp_d0m-3.096)<0.028&&pbarp_d0d0pidmu>0.5&&pbarp_d0d1pidmu>0.5";
	sel45[31]  = "abs(fpbarp_d0m-3.096)<0.028&&pbarp_d0d0pidmu>0.5&&pbarp_d0d1pidmu>0.5";	
	
	sel55[0]   = "abs(fpbarp_d0m-2.985)<0.075&&pbarp_d0d0pidk>0.5&&pbarp_d0d1pidk>0.5&&fchi2<15&&esptmax>0.9&&abs(pbarp_d0d2m-0.1342)<0.011";
	sel55[1]   = "abs(fpbarp_d0m-2.985)<0.08&&pbarp_d0d0pidk>0.5&&pbarp_d0d1pidk>0.5&&fchi2<15&&abs(pbarp_d0d2m-0.1342)<0.011&&esptmax>0.8";	
	sel55[10]  = "";
	sel55[11]  = "";	
	sel55[20]  = "abs(fpbarp_d0m-3.096)<0.06&&pbarp_d0d0pide>0.25&&pbarp_d0d1pide>0.25";
	sel55[21]  = "abs(fpbarp_d0m-3.096)<0.06&&pbarp_d0d0pide>0.25&&pbarp_d0d1pide>0.25";	
	sel55[30]  = "abs(fpbarp_d0m-3.096)<0.04&&pbarp_d0d0pidmu>0.5&&pbarp_d0d1pidmu>0.5";
	sel55[31]  = "abs(fpbarp_d0m-3.096)<0.04&&pbarp_d0d0pidmu>0.25&&pbarp_d0d1pidmu>0.25";	
	
	
}

double signi(double t, double effs, double effb, double fBR, double L, double sigs_nb, double sigb_mb=60.)
{
	double val = sqrt(L*t)*(sigs_nb*1e-33*effs*fBR)/sqrt(sigs_nb*sigs_nb*1e-66*effs*effs*fBR*fBR + sigb_mb*sigb_mb*1e-54*effb*effb);
	return val;
}

double t_until_signi(double Ns, double effs, double effb, double fBR, double L, double sigs_nb, double sigb_mb=60.)
{
	double val = Ns*Ns/L*(sigs_nb*1e-33*effs*fBR+sigb_mb*1e-27*effb)/(sigs_nb*sigs_nb*1e-66*effs*effs*fBR*fBR);
	return val/86400.;
}

double comptimes(double &dtl, double &dth, double Ns, double effs, double effb, double fBR, double L, double sigs_nb, double sigb_mb=60.)
{
	double t = t_until_signi(Ns, effs, effb, fBR, L, sigs_nb, sigb_mb);
	dtl      = t - t_until_signi(Ns, effs, effb, fBR, L, sigs_nb*2., sigb_mb);
	dth      = t_until_signi(Ns, effs, effb, fBR, L, sigs_nb*0.5, sigb_mb) - t;
	
	return t;
}

void calc_fom_gr(double sigs_nb=1., double Ns = 5., double sigb_mb=60.)
{
	setStyle();
	init_channels();
	// dmode = 0 (etac -> K K pi0), 1 (etac -> KS K pi), 2 (J/psi -> e+ e-), 3 (J/psi -> mu+ mu-)
	// rmode = 0 (pi+ pi+), 1 (pi0 pi0)
	
	// mode =          		00		01		10		11		20		21		30		31
	double eff_s_45[8] = { 0.23,    0.17,   0.172,   0.14,    0.316,  0.252,  0.375,  0.307  };
	double eff_b_45[8] = { 2e-6,  4.8e-7, 4.8e-5,  1.3e-5,    1e-8,   1e-8,   7e-8,   4e-8  };

	double eff_s_55[8] = { 0.23,   0.174,  0.174,   0.14,    0.318,  0.262,  0.359,   0.353  };
	double eff_b_55[8] = { 4.5e-6, 1.5e-6,  1.8e-5,  2.8e-6,    1e-8,   1e-8,   5e-8,    5e-8  };
	
	double fBR[8]      = { 0.035,  0.035,  0.035,  0.035,   0.058,  0.058,  0.058,   0.058 };
	
	TH1F *fom_hlum45 = new TH1F("fom_hlum45",TString::Format("Time for %0.0f#sigma: #sigma_{s} = %4.1fnb, #sqrt{s} = 4.5GeV",Ns,sigs_nb),8,0,8);
/*	TH1F *fom_llum45 = new TH1F("fom_llum45",TString::Format("Time for %0.0f#sigma: L = %s/cm^{2}s, #sigma_{s} = %4.1fnb, #sqrt{s} = 4.5GeV",Ns,"1·10^{31}",sigs_nb),8,0,8);*/
	
	TH1F *fom_hlum55 = new TH1F("fom_hlum55",TString::Format("Time for %0.0f#sigma: #sigma_{s} = %4.1fnb, #sqrt{s} = 5.5GeV",Ns,sigs_nb),8,0,8);
/*	TH1F *fom_llum55 = new TH1F("fom_llum55",TString::Format("Time for %0.0f#sigma: L = %s/cm^{2}s, #sigma_{s} = %4.1fnb, #sqrt{s} = 5.5GeV",Ns,"1·10^{31}",sigs_nb),8,0,8);*/
	
	TGraphAsymmErrors *gh45=new TGraphAsymmErrors(8);
	TGraphAsymmErrors *gl45=new TGraphAsymmErrors(8);
	TGraphAsymmErrors *gh55=new TGraphAsymmErrors(8);
	TGraphAsymmErrors *gl55=new TGraphAsymmErrors(8);
	
	config_graph(gh45, "gh45");
	config_graph(gl45, "gl45", 6);
	config_graph(gh55, "gh55");
	config_graph(gl55, "gl55",6);
	
	double max45=0, max55=0;
	
	for (int i=0;i<8;++i)
	{		
		double t, dtl, dth;
		// sqrts = 4.5
		t  = comptimes(dtl, dth, Ns, eff_s_45[i], eff_b_45[i], fBR[i], 2e32, sigs_nb, sigb_mb);
		gh45->SetPoint(i,0.35+i,t); gh45->SetPointError(i,0,0,dtl,dth);
				
		t  = comptimes(dtl, dth, Ns, eff_s_45[i], eff_b_45[i], fBR[i], 1e31, sigs_nb, sigb_mb);
		gl45->SetPoint(i,0.65+i,t); gl45->SetPointError(i,0,0,dtl,dth);
		
		if ((t+dth)>max45) max45 = t+dth;
		
		// sqrts = 5.5
		t  = comptimes(dtl, dth, Ns, eff_s_55[i], eff_b_55[i], fBR[i], 2e32, sigs_nb, sigb_mb);
		gh55->SetPoint(i,0.35+i,t); gh55->SetPointError(i,0,0,dtl,dth);
		
		t  = comptimes(dtl, dth, Ns, eff_s_55[i], eff_b_55[i], fBR[i], 1e31, sigs_nb, sigb_mb);
		gl55->SetPoint(i,0.65+i,t); gl55->SetPointError(i,0,0,dtl,dth);
		
		if ((t+dth)>max55) max55 = t+dth;		
		
		
		fom_hlum45->GetXaxis()->SetBinLabel(i+1,channels[i]);
/*		fom_llum45->GetXaxis()->SetBinLabel(i+1,channels[i]);*/
		fom_hlum55->GetXaxis()->SetBinLabel(i+1,channels[i]);
/*		fom_llum55->GetXaxis()->SetBinLabel(i+1,channels[i]);*/
		
/*		fom_hlum45->SetBinContent(i+1,t_until_signi(Ns, eff_s_45[i], eff_b_45[i], fBR[i], 2e32, sigs_nb, sigb_mb));
		fom_llum45->SetBinContent(i+1,t_until_signi(Ns, eff_s_45[i], eff_b_45[i], fBR[i], 1e31, sigs_nb, sigb_mb));
		
		fom_hlum55->SetBinContent(i+1,t_until_signi(Ns, eff_s_55[i], eff_b_55[i], fBR[i], 2e32, sigs_nb, sigb_mb));
		fom_llum55->SetBinContent(i+1,t_until_signi(Ns, eff_s_55[i], eff_b_55[i], fBR[i], 1e31, sigs_nb, sigb_mb));*/
	}

	config_histo(fom_hlum45);
	//config_histo(fom_llum45);
	config_histo(fom_hlum55);
	//config_histo(fom_llum55);

	fom_hlum45->SetMaximum(max45*2);
	fom_hlum45->SetMinimum(0.005);
	//fom_llum45->SetMinimum(0.01);
	fom_hlum55->SetMaximum(max55*2);
	fom_hlum55->SetMinimum(0.005);
	//fom_llum55->SetMinimum(0.01);

	TCanvas *c1=new TCanvas("c1","c1",10,10,1200,600);
	c1->Divide(2,1);
	
	int limit1 = 30;
	int limit2 = 360;
	
	TLine lgood, lbad;
	TLatex latgood, latbad;
	TBox boxg, boxy, boxr;
	boxg.SetFillStyle(3003);boxg.SetFillColor(kGreen);
	boxy.SetFillStyle(3003);boxy.SetFillColor(kOrange);
	boxr.SetFillStyle(3003);boxr.SetFillColor(kRed);

	lgood.SetLineColor(kGreen+3);lgood.SetLineStyle(2);lgood.SetLineWidth(2);
	lbad.SetLineColor(2);lbad.SetLineStyle(2);lbad.SetLineWidth(2);

	latgood.SetTextColor(kGreen+3);latgood.SetTextSize(0.035);
	latbad.SetTextColor(2);latbad.SetTextSize(0.035);
	
	c1->cd(1); gPad->SetLogy();
	fom_hlum45->Draw("P");
	lgood.DrawLine(0,limit1,8,limit1);lbad.DrawLine(0,limit2,8,limit2);
	boxg.DrawBox(0,0.005,8,30); boxy.DrawBox(0,30,8,360); boxr.DrawBox(0,360,8,fom_hlum45->GetMaximum());
	gh45->Draw("Psame");
	gl45->Draw("Psame");
	TLegend *leg1=new TLegend(0.75,0.75,0.98,0.87);
	leg1->AddEntry("gh45","L = 2·10^{32}","p");
	leg1->AddEntry("gl45","L = 1·10^{31}","p");
	leg1->Draw();
	latgood.DrawLatex(2,30*1.5,"1 month");
	latbad.DrawLatex(2,365*1.5,"1 year");
	
	//fom_hlum45->Draw("Psame");
	
	
/*	c1->cd(2); gPad->SetLogy();
	fom_llum45->Draw("P");
	lgood.DrawLine(0,limit1,8,limit1);lbad.DrawLine(0,limit2,8,limit2);
	boxg.DrawBox(0,0.01,8,30); boxy.DrawBox(0,30,8,360); boxr.DrawBox(0,360,8,fom_llum45->GetMaximum()*2);
	gl45->Draw("Psame");
	fom_llum45->Draw("Psame");
*/
	
	c1->cd(2); gPad->SetLogy();
	fom_hlum55->Draw("P");
	lgood.DrawLine(0,limit1,8,limit1);lbad.DrawLine(0,limit2,8,limit2);
	boxg.DrawBox(0,0.005,8,30); boxy.DrawBox(0,30,8,360); boxr.DrawBox(0,360,8,fom_hlum55->GetMaximum());
	gh55->Draw("Psame");
	gl55->Draw("Psame");
	TLegend *leg2=new TLegend(0.75,0.75,0.98,0.87);
	leg2->AddEntry("gh55","L = 2·10^{32}","p");
	leg2->AddEntry("gl55","L = 1·10^{31}","p");
	leg2->Draw();
	latgood.DrawLatex(2,30*1.5,"1 month");
	latbad.DrawLatex(2,365*1.5,"1 year");
	//fom_hlum55->Draw("Psame");
/*
	c1->cd(4); gPad->SetLogy();
	fom_llum55->Draw("P");
	lgood.DrawLine(0,limit1,8,limit1);lbad.DrawLine(0,limit2,8,limit2);
	boxg.DrawBox(0,0.01,8,30); boxy.DrawBox(0,30,8,360); boxr.DrawBox(0,360,8,fom_llum55->GetMaximum()*2);
	gl55->Draw("Psame");
	fom_llum55->Draw("Psame");
	*/
}