#include <TGraph.h>
#include "constants.h"
double mass[]  = {70, 80, 90, 100, 110, 120, 130, 140, 150};
double br_LEP[]= {0.012,  0.014,  0.02,  0.03,  0.1,  0.8,  100,  100,  100};
//double br[]       = {0.81 , 0.7 ,0.4 , 0.18, 0.06, 0.028, 0.01, 0.006, 0.002};
double br[]       = {0.81 , 0.7 ,0.4 , 0.18, 0.07, 0.032, 0.012, 0.0068, 0.0038};
double br_sm[]       = {0.001, 0.0013, 0.0017, 0.002 ,0.0026, 0.003, 0.0027, 0.0025, 0.0019};

double sigma_lo[] = {0.8, 0.5,0.30,0.10, 0.03, 0.02,0.01, 0.002,0.001};
double sigma[]    = {1.47, 1.04, 0.75, 0.55, 0.42, 0.32, 0.25, 0.19, 0.15};
double err_sigma[] = {0.039, 0.029, 0.022, 0.016, 0.013, 0.010, 0.008, 0.008, 0.006};

double run1[] = {0.46, 0.44, 0.37, 0.35, 0.34, 0.33,0.33,0.32,0.32};

double theory[9], theory_lo[9], err_theory[9], br_run1[9], br_run2[9], br_run2_exp[9];
double mass2[] = {70,80,90,100,110,120,130.5,140.5,150.5};

plot_limit(){
  
  // Read Expected and Observed Limits
  double expected[9];
  double observed[9];

  double sm_limit_obs[9];
  double sm_limit_exp[9];
  double sm_limit_exp_b[9];

  ifstream f;
  f.open("limit.txt",ios::in);
  string line;
  getline (f,line);
  cout<<line<<endl;
  for(int i=0;i<9;++i){
    f>>mass[i]>>observed[i]>>expected[i];
    observed[i];
    expected[i];
  }
  f.close(); 

  for(int i = 0;i<9;++i)
    {
      theory_lo[i]  = sigma_lo[i] *br[i];
      theory[i]     = sigma[i]    *br[i];  
      
      br_run1[i] = run1[i]    /sigma_lo[i];
      br_run2[i] = observed[i]/sigma[i];
      br_run2_exp[i] = expected[i]/sigma[i];

      if (br_run1[i]>1) br_run1[i] =100;
      if (br_run2[i]>1) br_run2[i] =100;
      sm_limit_obs[i]  = observed[i]/(br_sm[i]*sigma[i]);
      sm_limit_exp [i] = expected[i]/(br_sm[i]*sigma[i]);
      sm_limit_exp_b [i] = expected[i]/(1.41*br_sm[i]*sigma[i]);
      printf("obs: = %.3g, exp = %.3g\n",sm_limit_obs[i], sm_limit_exp_b[i]);
     }
  for(int i = 0;i<9;++i)
       printf("%.3g : %.3g\n",mass[i],br_LEP[i]/br_run2[i]);
  
  
  for(int i=0;i<9;++i)
    printf(" %.3g \n ",br_run1[i]);
  
  TGraph * g_th_lo = new TGraph(9,mass, theory_lo);
  //   float t[2];
  //   g_th_lo->GetX(2,t[0],t[1]);
  //   printf("%.3g %.3g %.3g\n", theory_lo[0], t[0],t[1]);
  //  exit(0);

  TGraph * g_th = new TGraph(9,mass, theory);
  TGraph * g_run1 = new TGraph(9,mass,run1);
  TGraph * g_exp = new TGraph(9,mass,expected);
  TGraph * g_obs = new TGraph(9,mass,observed);

  TGraph * g_br_lep = new TGraph(9, mass2, br_LEP);
  TGraph * g_br_run1 = new TGraph(9, mass, br_run1);
  TGraph * g_br_run2 = new TGraph(9, mass, br_run2);
  TGraph * g_br_run2_exp = new TGraph(9, mass, br_run2_exp);


  TGraph * g_br = new TGraph(9, mass, br);

  TGraph * g_sm_obs = new TGraph(9, mass, sm_limit_obs);
  TGraph * g_sm_exp = new TGraph(9, mass, sm_limit_exp);
  TGraph * g_sm_exp_b = new TGraph(9, mass, sm_limit_exp_b);
  
  TLegend *leg = new TLegend(0.14,0.186,0.492,0.456);
  leg->AddEntry(g_th_lo,"B*#sigma, LO");
  leg->AddEntry(g_th,"B*#sigma, NLO");
  leg->AddEntry(g_run1,"B*#sigma^{95 %CL}_{Run I}");
  leg->AddEntry(g_exp,"B*#sigma^{95 %CL}_{Run II} expected");
  leg->AddEntry(g_obs,"B*#sigma^{95 %CL}_{Run II} observed");
  leg->SetFillColor(kWhite);

  Color_t white=10;
  gStyle->SetCanvasColor(white);
  
  TCanvas *c = new TCanvas("c","", 800, 600);
  c->SetGridy();
  c->SetLogy();
/*

  g_th_lo->GetYaxis()->SetRangeUser(0.000001,2.5);
  g_th_lo->SetTitle("");
  g_th_lo->GetXaxis()->SetTitle("m_{h_{f}} (GeV)");
  g_th_lo->GetYaxis()->SetTitle("#sigma_{h_{f}} B(h_{f}#rightarrow #gamma#gamma), pb");
  g_th_lo->SetLineColor(kRed);
  g_th_lo->SetLineWidth(2);
  g_th_lo->SetLineStyle(2);
  g_th_lo->SetFillColor(kWhite);
  g_th_lo->Draw("AC");
  
  g_th->SetTitle("");
  g_th->GetXaxis()->SetTitle("m_{h_{f}} (GeV)");
  g_th->SetLineColor(kRed);
  g_th->SetLineWidth(2);
  g_th->SetFillColor(kWhite);
  g_th->Draw("CSAME");
  
  g_run1->SetMarkerStyle(21);
  g_run1->SetMarkerColor(kBlack);
  g_run1->SetFillColor(kWhite);
  g_run1->SetMarkerSize(1.1);
  g_run1->Draw("Psame");
  
  g_exp->SetMarkerStyle(23);
  g_exp->SetMarkerColor(kGreen);
  g_exp->SetFillColor(kWhite);
  g_exp->SetMarkerSize(1.1);
  g_exp->Draw("Psame");
  
  g_obs->SetMarkerStyle(22);
  g_obs->SetMarkerColor(4);
  g_obs->SetFillColor(kWhite);
  g_obs->SetMarkerSize(1.1);
  
  g_obs->Draw("Psame");
  
  leg->Draw();

  label("D\351, 1.1 fb^{ -1}");

 //   c->Print("figures/limits_graphs.gif");
//   c->Print("figures/limits_graphs.pdf");
//   c->Print("figures/limits_graphs.root");

  //  c->Print("figures/fig_4.gif");
  //  c->Print("figures/fig_4.pdf");
  //  c->Print("figures/fig_4.root");


*/
  // Plot B's
  //  TLegend *leg = new TLegend(0.63,0.33,0.98,0.60);
  TLegend *leg = new TLegend(0.7,0.38,0.87,0.648);
  leg->AddEntry(g_br_lep,"LEP");
  leg->AddEntry(g_br_run1,"D\351, Run I");
  leg->AddEntry(g_br_run2,"observed");
  leg->AddEntry(g_br_run2_exp,"expected");
  leg->AddEntry(g_br,"Benchmark");
  
  leg->SetFillColor(kWhite);
  leg->SetTextSize(label_font_size);
  
  Color_t white=10;
  gStyle->SetCanvasColor(white);
  
  TCanvas *c = new TCanvas("c","", 800, 600);
  //  c->SetGridy();
  
  c->SetLogy();


  
  // Format fonts
  g_br_lep->GetXaxis()->SetLabelSize(label_font_size);
  g_br_lep->GetYaxis()->SetLabelSize(label_font_size);
  
  g_br_lep->GetXaxis()->SetTitleSize(label_font_size);
  g_br_lep->GetYaxis()->SetTitleSize(label_font_size);
  
  g_br_lep->SetTitle("");
  g_br_lep->GetXaxis()->SetTitle("m_{h_{f}} (GeV)");
  
  g_br_lep->GetXaxis()->CenterTitle();
  g_br_lep->GetYaxis()->CenterTitle();
  g_br_lep->GetXaxis()->SetTitleOffset(0.9);
  g_br_lep->GetYaxis()->SetTitle("B(h_{f}#rightarrow #gamma#gamma)");
  g_br_lep->SetMarkerColor(kRed);
  g_br_lep->SetMarkerStyle(20);
  //g_br_lep->SetFillColor(kWhite);
  g_br_lep->SetFillColor(kRed);
  g_br_lep->SetFillStyle(3004);
  g_br_lep->SetLineWidth(8/*2*/002);
  g_br_lep->SetLineColor(kRed);

  g_br_lep->GetYaxis()->SetRangeUser(0.001,2);
  g_br_lep->SetMarkerSize(1.3);


// h1->SetContour(2);
// h1->SetContourLevel(0,level1);  //value for your first level
// h1->SetContourLevel(1,1e6);     //non-existing high level
// h1->Draw("cont1");
// h2->SetContour(2);
// h2->SetContourLevel(0,0);       //non existing low level
// h2->SetContourLevel(1,level2);  //value for your second level 
// h2->Draw("cont1 same");
//g_br_lep->SetContour(2);
//g_br_lep->SetContourLevel(0,1);  //value for your first level
//g_br_lep->SetContourLevel(1,1e6);	  //non-existing high level


//g_br_lep->Draw("APcont2");
//g_br_lep->Draw("APsurf");

  g_br_lep->Draw("APC
");
 
// 	TMultiGraph *mg = new TMultiGraph();
//  
// 	mg->Add(g_br_lep);
// 	//mg->Add(g_br_run1);
// 	//mg->Add(g_br_run2);
// 	//mg->Add(g_br_run2_exp);
// 	//mg->Add(g_br);
// 	mg->Draw("APC");
 
  g_br_run1->SetMarkerStyle(21);
  g_br_run1->SetMarkerColor(kBlack);
  g_br_run1->SetFillColor(kWhite);
  g_br_run1->SetMarkerSize(1.3);
  g_br_run1->SetFillColor(kBlack);
  g_br_run1->SetFillStyle(3001);
  g_br_run1->SetLineWidth(6000);
  g_br_run1->SetLineColor(kBlack);
  g_br_run1->Draw("Plsame");
 
 
  g_br_run2->SetMarkerStyle(23);
  g_br_run2->SetMarkerColor(kBlue);
  g_br_run2->SetFillColor(kWhite);
  g_br_run2->SetMarkerSize(1.3);
  g_br_run2->SetFillColor(kBlue);
  g_br_run2->SetFillStyle(3002);
  g_br_run2->SetLineWidth(8/*2*/002);
  g_br_run2->SetLineColor(kBlue);


  g_br_run2->Draw("PCsame");
  
  g_br_run2_exp->SetMarkerStyle(26);
  g_br_run2_exp->SetMarkerColor(kBlue);
  g_br_run2_exp->SetFillColor(kWhite);
  g_br_run2_exp->SetMarkerSize(1.3);
  g_br_run2_exp->Draw("Psame");
  
  //  g_br->SetMarkerStyle(24);
  g_br->SetLineColor(kGreen);
  g_br->SetFillColor(kWhite);
  g_br->SetLineWidth(2);
  g_br->Draw("PCsame");



// TMultiGraph *mg = new TMultiGraph();
// 
// mg->Add(g_br_lep);
// //mg->Add(g_br_run1);
// //mg->Add(g_br_run2);
// //mg->Add(g_br_run2_exp);
// //mg->Add(g_br);
// mg->Draw("APCSame");

  //leg->SetShadowColor(0);
  leg->SetLineColor(0);
  leg->SetTextSize(0.04);
  leg->Draw();

 






// 
// 
//    TMultiGraph *mg = new TMultiGraph();
//    mg->SetTitle("Exclusion graphs");
// 
//    const Int_t n = 35;
//    Double_t x1[n], x2[n], x3[n], y1[n], y2[n], y3[n];
//    for (Int_t i=0;i<n;i++) {
//      x1[i]  = i*0.1;
//      x2[i]  = x1[i];
//      x3[i]  = x1[i]+.5;
//      y1[i] = 10*sin(x1[i]);
//      y2[i] = 10*cos(x1[i]);
//      y3[i] = 10*sin(x1[i])-2;
//    }
// 
//    gr1 = new TGraph(n,x1,y1);
//    gr1->SetLineColor(2);
//    gr1->SetLineWidth(1504);
//    gr1->SetFillStyle(3005);
// 
//    gr2 = new TGraph(n,x2,y2);
//    gr2->SetLineColor(4);
//    //gr2->SetLineWidth(-2002);
//    gr2->SetLineWidth(-4002);
// 
//    gr2->SetFillStyle(3004);
//    gr2->SetFillColor(9);
// 
//    gr3 = new TGraph(n,x3,y3);
//    gr3->SetLineColor(5);
//    gr3->SetLineWidth(-802);
//    gr3->SetFillStyle(3002);
//    gr3->SetFillColor(2);
// 
//    mg->Add(gr1);
//    mg->Add(gr2);
//    mg->Add(gr3);
//    mg->Draw("AC");


 
  TPaveText *pt  =  new TPaveText(0.18,0.20,.47,.28,"NDC");
  pt->SetTextSize(label_font_size);
  pt->AddText("D\351, 1.1 fb^{ -1}");
  pt->SetFillColor(kWhite);
  pt->Draw("same");

   TPave *pave = new TPave(62.112,1.0,158.028,1.99,0,"br");
   pave->SetFillColor(0);
   pave->Draw("same");
  
  //	label("D\351, 1.1 fb^{ -1}");
  
  //  c->SetLogy();
	
//   c->Print("figures/br_limits_graphs.gif");
//   c->Print("figures/br_limits_graphs.pdf");
//   c->Print("figures/br_limits_graphs.root"); 

  c->Print("figures/fig_5.gif");
  c->Print("figures/fig_5.pdf");
  c->Print("figures/fig_5.root");
  c->Print("figures/fig_5.eps");

//   // Produce limit over sm expectation
//   delete leg;
//   
//   TLegend *leg = new TLegend(0.23,0.40,0.55,0.67);
//   leg->AddEntry(g_sm_obs,"observed");
//   leg->AddEntry(g_sm_exp,"expected");
//   leg->AddEntry(g_sm_exp_b,"Run2a+b extrapolation");
//   leg->SetFillColor(kWhite);
// 
//   Color_t white=10;
//   gStyle->SetCanvasColor(white);
//   delete c;
	
//   TCanvas *c = new TCanvas("c","", 800, 600);
// 
//   g_sm_exp->SetTitle("");
//   g_sm_exp->GetXaxis()->SetTitle("m_{h^{SM}} (GeV)");
//   
//   g_sm_exp->GetXaxis()->CenterTitle();
//   g_sm_exp->GetXaxis()->SetTitleOffset(1.2);
//   g_sm_exp->GetYaxis()->SetTitle("(B(h_{SM}#rightarrow #gamma#gamma)#times #sigma_{SM})/(SM prediction)");
//   g_sm_exp->SetLineWidth(3);
//   g_sm_exp->SetLineStyle(2);
//   g_sm_exp->SetLineColor(kRed);
// 	g_sm_exp->SetFillColor(kWhite);
//   g_sm_exp->GetYaxis()->SetRangeUser(40,200);
// 	 g_sm_exp->GetXaxis()->SetRangeUser(95,155);
//   g_sm_exp->Draw("AL");
// 	
//   g_sm_exp_b->GetXaxis()->CenterTitle();
//   g_sm_exp_b->GetXaxis()->SetTitleOffset(1.2);
// 	g_sm_exp_b->SetLineWidth(3);
//   g_sm_exp_b->SetLineStyle(2);
//   g_sm_exp_b->SetLineColor(kBlue);
//   g_sm_exp_b->SetFillColor(kWhite);
//   g_sm_exp_b->GetYaxis()->SetRangeUser(0,200);
//   g_sm_exp_b->Draw("Lsame");
// 	
//   g_sm_obs->GetXaxis()->CenterTitle();
//   g_sm_obs->GetXaxis()->SetTitleOffset(1.2);
//   g_sm_obs->SetLineWidth(3);
//   g_sm_obs->SetLineColor(kBlack);
//   g_sm_obs->SetFillColor(kWhite);
//   g_sm_obs->GetYaxis()->SetRangeUser(0,200);
//   g_sm_obs->SetMarkerSize(1.3);
//   g_sm_obs->Draw("Lsame");
//   leg->Draw();
//   label("VV#rightarrow h, V#rightarrowhV, D#oslash Run IIa, 1.1 fb^{ -1}"); 
}