TH1D *background;
TH1D *signal;
TH1D *data;
TH1D *L;
TH1D *L95;

TCanvas *c1;
TGaxis *axis;
//gStyle->SetOptStat(kFALSE);
Float_t rightmax;
Float_t scale;

const int sigsteps=300;
const double maxsig = 5.;//25.; //in pb
const double minsig = 0.1;
const double masscut=0.;
//const double masscut=350.; //Z'(m=500)
//const double masscut=600.; //Z'(m=750)
//const double masscut=800.; //Z'(m=1000)
//const double masscut=900.; //Z'(m=1250)
//const double masscut=1100.;//Z'(m=1500)
//const double masscut=1600.;//Z'(m=2000)

void ReadMe();

void InitializeMe(){
  L= new TH1D("L","L",sigsteps,minsig,maxsig); 
  L95= new TH1D("L95","L95",sigsteps,minsig,maxsig);
  L->SetTitle("");
  L->GetYaxis()->SetTitleSize(0.05);
  L->GetYaxis()->SetTitleOffset(1.40);
  L->GetYaxis()->SetTitleColor(kBlack);
  L->GetYaxis()->SetTitle("L"); 
  L->GetXaxis()->SetTitleSize(0.05);
  L->GetXaxis()->SetTitleOffset(1.18);
  L->GetXaxis()->SetTitleColor(kBlack);
  L->GetXaxis()->SetTitle("#sigma(PP #rightarrow Z') x BR(Z'#rightarrow t#bar{t}) [pb]"); 
}
void ReadMe(){

  TFile f("MC_Sig_500GeV_Bck.root"); 
  //TFile f("MC_Sig_750GeV_Bck.root"); 
  //TFile f("MC_Sig_1000GeV_Bck.root");
  //TFile f("MC_Sig_1250GeV_Bck.root");
  //TFile f("MC_Sig_1500GeV_Bck.root");
  //TFile f("MC_Sig_2000GeV_Bck.root");

  background = new TH1D(*((TH1D*)f.Get("h_all")));

  signal = new TH1D(*((TH1D*)f.Get("zprime_m500")));
  //signal = new TH1D(*((TH1D*)f.Get("zprime_m750")));
  //signal = new TH1D(*((TH1D*)f.Get("zprime_m1000")));
  //signal = new TH1D(*((TH1D*)f.Get("zprime_m1250")));
  //signal = new TH1D(*((TH1D*)f.Get("zprime_m1500")));
  //signal = new TH1D(*((TH1D*)f.Get("zprime_m2000")));

  signal->Scale(1./50.); //signal corresponds to 50 pb crossection
  data = new TH1D(*((TH1D*)f.Get("data")));
  background->SetDirectory(gROOT);
  signal->SetDirectory(gROOT);
  data->SetDirectory(gROOT);
  f.Close();
  InitializeMe();
}


void LikelihoodMe(double xsec){
  double l=1.;
  for(int i=signal->FindBin(masscut); i<= signal->FindBin(2500.);i++){
    double s=signal->GetBinContent(i)*xsec;
    double b=background->GetBinContent(i);
    double data=data->GetBinContent(i);
    //l*=TMath::Poisson(data,s+b);//observed limit
    l*=TMath::Poisson(b,s+b);   //expected limit
  }
  int sigbin=L->FindBin(xsec);
  L->SetBinContent(sigbin,l);
}

void DoMe(){
  ReadMe();
  for(int i=1;i<=sigsteps;i++){
    LikelihoodMe(L->GetBinCenter(i));
  }
  double ul95=0;
  for(int i=1;L->Integral(1,i)/L->Integral(1,sigsteps)<0.95&&i<=sigsteps;i++){
    L95->SetBinContent(i,L->GetBinContent(i));
    ul95=L95->GetBinCenter(i);
  }



  double _margin_top    = 0.070;  double _margin_right  = 0.055;
  double _margin_left   = 0.185;  double _margin_bottom = 0.160;
  
  c1 = new TCanvas("c1","transparent pad",200,10,400,400);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  c1->SetHighLightColor(1);
  c1->SetFillColor(kWhite);
  c1->SetBorderMode(0);
  c1->SetBorderSize(0);
  c1->SetFrameBorderMode(0);
  c1->SetFrameBorderSize(0);
  c1->Range(0,0,1,1);
  c1->SetTopMargin(_margin_top);
  c1->SetRightMargin(_margin_right);
  c1->SetLeftMargin(_margin_left);
  c1->SetBottomMargin(_margin_bottom);
  //c1->SetLogy();

  L95->SetFillStyle(3004);
  L95->SetFillColor(2);//expected limit
  //L95->SetFillColor(4);//observed limit
  L->Draw();  
  L95->Draw("SAME");      
  //c1->Print("Likelihood500.pdf");
  
 TPaveText *tpt = new TPaveText(.72,.80,.93,.91,"brNDC");
   
  tpt->SetName("tpavetext");
  //tpt->SetBorderSize(0);
  tpt->SetFillColor(10);
  tpt->SetShadowColor(10);
  tpt->SetTextAlign(22);
  tpt->SetTextFont(42);
  tpt->SetTextColor(kBlack);
  tpt->SetTextSize(0.037);

  //observed Limit     
  //tpt->AddText("M_{z'}=0.5 TeV        ");
  //tpt->AddText("UL(95%)=14 pb");
  //tpt->AddText("M_{z'}=0.75 TeV        ");
  //tpt->AddText("UL(95%)=4.27 pb");
  //tpt->AddText("M_{z'}=1 TeV               ");
  //tpt->AddText("UL(95%)=1.27 pb");
  //tpt->AddText("M_{z'}=1.25 TeV         ");
  //tpt->AddText("UL(95%)=0.99 pb");
  //tpt->AddText("M_{z'}=1.5 TeV          ");
  //tpt->AddText("UL(95%)=0.86 pb");
  //tpt->AddText("M_{z'}=2 TeV              ");
  //tpt->AddText("UL(95%)=0.75 pb");

  //expected Limit     
  tpt->AddText("M_{z'}=0.5 TeV          ");
  tpt->AddText("UL(95%)=5.46 pb");
  //tpt->AddText("M_{z'}=0.75 TeV        ");
  //tpt->AddText("UL(95%)=3.15 pb");
  //tpt->AddText("M_{z'}=1 TeV               ");
  //tpt->AddText("UL(95%)=2.28 pb");
  //tpt->AddText("M_{z'}=1.25 TeV         ");
  //tpt->AddText("UL(95%)=1.77 pb");
  //tpt->AddText("M_{z'}=1.5 TeV          ");
  //tpt->AddText("UL(95%)=1.60 pb");
  //tpt->AddText("M_{z'}=2 TeV              ");
  //tpt->AddText("UL(95%)=1.22 pb");
  
  tpt->Draw();
 
  std::cout<<"95 % UL:"<<ul95<<std::endl;
}