#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 "TEventList.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.22);
	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(2.);
	h->SetMarkerStyle(20);
	h->SetMarkerColor(4);
}

void init_channels()
{
	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]   = "";
	sel55[1]   = "";	
	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;
}

void calc_fom(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.99,   0.99,   0.174,   0.14,    0.318,  0.262,  0.359,   0.353  };
	double eff_b_55[8] = { 1e-13,  1e-13,  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: L = %s/cm^{2}s, #sigma_{s} = %4.1fnb, #sqrt{s} = 4.5GeV",Ns,"2·10^{32}",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: L = %s/cm^{2}s, #sigma_{s} = %4.1fnb, #sqrt{s} = 5.5GeV",Ns,"2·10^{32}",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);
	
	
	for (int i=0;i<8;++i)
	{
		fom_hlum45->GetXaxis()->SetBinLabel(i+1,channels[i]);
		fom_llum45->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)/86400);
		fom_llum45->SetBinContent(i+1,t_until_signi(Ns, eff_s_45[i], eff_b_45[i], fBR[i], 1e31, sigs_nb, sigb_mb)/86400);
		
		fom_hlum55->GetXaxis()->SetBinLabel(i+1,channels[i]);
		fom_llum55->GetXaxis()->SetBinLabel(i+1,channels[i]);
		fom_hlum55->SetBinContent(i+1,t_until_signi(Ns, eff_s_55[i], eff_b_55[i], fBR[i], 2e32, sigs_nb, sigb_mb)/86400);
		fom_llum55->SetBinContent(i+1,t_until_signi(Ns, eff_s_55[i], eff_b_55[i], fBR[i], 1e31, sigs_nb, sigb_mb)/86400);
	}

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

	fom_hlum45->SetMinimum(0.01);
	fom_llum45->SetMinimum(0.01);
	fom_hlum55->SetMinimum(0.01);
	fom_llum55->SetMinimum(0.01);

	TCanvas *c1=new TCanvas("c1","c1",10,10,800,600);
	c1->Divide(2,2);
	
	int limit1 = 30;
	int limit2 = 360;
	
	TLine lgood, lbad;
	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);

	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.01,8,30); boxy.DrawBox(0,30,8,360); boxr.DrawBox(0,360,8,fom_hlum45->GetMaximum()*2);
	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);
	fom_llum45->Draw("Psame");
	
	c1->cd(3); gPad->SetLogy();
	fom_hlum55->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_hlum55->GetMaximum()*2);
	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);
	fom_llum55->Draw("Psame");
	
}