#include "Math/Interpolator.h"
#include <iostream>
#include "AtlasStyle.C"
void BRHaalimitSetHigh_low_mass() {

  TString files_lm[]={"./LowMass_a0a0.txt"};
  TString files_hm[]={"./HighMass_a0a0.txt"};
  TString files[]={"./type-1.dat"};
  // TString s[]={"3"};
   TString canv[]={"can3"};
  
  gSystem->Load("AtlasStyle.C");
  SetAtlasStyle();
  
  

  const unsigned int  Npoints_lm = 28;
  const unsigned int  Npoints_hm = 91;
  const unsigned int  Npoints = 673;
  const unsigned int  Npoints_max = 1000; 
  
  
  Double_t mzd[Npoints_max];
  Double_t Brall[Npoints_max];
  
  
  
  Double_t ggFSigma[Npoints_max];
  Double_t fac[Npoints_max];
  
  
  
  double mass[Npoints_max];
  double obs[Npoints_max];
  double Exp[Npoints_max];
  double onesiglow[Npoints_max];
  double onesighigh[Npoints_max];
  double twosiglow[Npoints_max];
  double twosighigh[Npoints_max];
  
  
  for(Int_t i=0;i<1;i++){
    
    TCanvas *can = new TCanvas(canv[i], canv[i],66,52,800,600);
    can->Range(0,0,1,1);
    can->SetFillColor(0);
    can->SetBorderMode(0);
    can->SetBorderSize(2);
    can->SetFrameBorderMode(0);
    
    //Plot the observed limits
    TGraphErrors *gre = new TGraphErrors(Npoints_max);
    gre->SetName("CLs_observed");
    gre->SetTitle("Observed CLs");
    gre->SetFillColor(1);
    gre->SetLineColor(1);
    gre->SetLineWidth(2);
    gre->SetMarkerStyle(20);
    
    // cout<<Exp[0]*(ggFSigma[0])<<endl;
    
    
    TH1F *Graph_CLs_observed1 = new TH1F("Graph_CLs_observed1","",100,1,60);
    Graph_CLs_observed1->SetMinimum(0.0e-4);
    Graph_CLs_observed1->SetMaximum(1000);
    Graph_CLs_observed1->SetDirectory(0);
    Graph_CLs_observed1->SetStats(0);
    TGaxis *yaxis = (TGaxis*)Graph_CLs_observed1->GetYaxis();
    yaxis->SetMaxDigits(2);
    
    Int_t ci;   // for color index setting
    Int_t cii;
    ci = TColor::GetColor("#000099");
    Graph_CLs_observed1->SetLineColor(ci);
    Graph_CLs_observed1->GetXaxis()->SetTitle("m_{Z_{d}} [GeV]");
    Graph_CLs_observed1->GetXaxis()->SetLabelFont(42);
    Graph_CLs_observed1->GetXaxis()->SetLabelSize(0.035);
    Graph_CLs_observed1->GetXaxis()->SetTitleSize(0.035);
    Graph_CLs_observed1->GetXaxis()->SetTitleFont(42);
    Graph_CLs_observed1->GetYaxis()->SetTitle("Upper Bound on BR [H #rightarrow Z_{d}Z_{d}]");
    Graph_CLs_observed1->GetYaxis()->SetLabelFont(42);
    Graph_CLs_observed1->GetYaxis()->SetLabelSize(0.035);
    Graph_CLs_observed1->GetYaxis()->SetTitleSize(0.035);
    Graph_CLs_observed1->GetYaxis()->SetTitleFont(42);
    Graph_CLs_observed1->GetYaxis()->SetTitleOffset(1.85);
    Graph_CLs_observed1->GetZaxis()->SetLabelFont(42);
    Graph_CLs_observed1->GetZaxis()->SetLabelSize(0.035);
    Graph_CLs_observed1->GetZaxis()->SetTitleSize(0.035);
    Graph_CLs_observed1->GetZaxis()->SetTitleFont(42);
    gre->SetHistogram(Graph_CLs_observed1);
    
    gre->Draw("al");

     //Plot for the expected limits and the pm 1 and 2 sigma fluctuations
    TMultiGraph *multigraph = new TMultiGraph();
    multigraph->SetName("Plot_expected");
    multigraph->SetTitle("Expected Plot");
   
    TString tmps;
    ifstream fin_lm(files_lm[i]);
    ifstream fin(files[i]);
    
    for (Int_t k=0; k<Npoints; k++)
      {
	fin>>tmps;mzd[k]=tmps.IsFloat()?tmps.Atof():0;fin>>tmps; Brall[k]=tmps.IsFloat()?tmps.Atof():0;
	
      }
    
    cout<<Brall[672]<<endl;
    ROOT::Math::Interpolator inter(673, ROOT::Math::Interpolation::kCSPLINE);
    inter.SetData(673, mzd, Brall);
    
    double invCrossSection = (1.0 / 48.59)*0.001;
    
    // fin>>tmps;fin>>tmps;fin>>tmps;fin>>tmps;fin>>tmps;fin>>tmps;fin>>tmps;//read first line.
    
    for (Int_t j=0; j<Npoints_lm; j++)
      {
      
      fin_lm>>tmps;mass[j]=tmps.IsFloat()?tmps.Atof():0;fin_lm>>tmps; obs[j]=tmps.IsFloat()?tmps.Atof():0; fin_lm>>tmps; Exp[j]=tmps.IsFloat()?tmps.Atof():0; fin_lm>>tmps; onesiglow[j]=tmps.IsFloat()?tmps.Atof():0; fin_lm>>tmps; onesighigh[j]=tmps.IsFloat()?tmps.Atof():0; fin_lm>>tmps; twosiglow[j]=tmps.IsFloat()?tmps.Atof():0; fin_lm>>tmps; twosighigh[j]=tmps.IsFloat()?tmps.Atof():0;

      double interBr = inter.Eval( mass[j] );
      fac[j] =  invCrossSection/ (interBr*interBr);
      ggFSigma[j] = fac[j]; 
      
      }

   TGraphAsymmErrors *grae = new TGraphAsymmErrors(Npoints_lm);
    grae->SetName("sig2");
    grae->SetTitle("Expected CLs #pm 2 #sigma");
    
    ci = TColor::GetColor("#ffff00");
    grae->SetFillColor(ci);
    
    for ( Int_t j=0; j<Npoints_lm; j++) {
      grae->SetPoint(j,mass[j],Exp[j]*(ggFSigma[j]));
      grae->SetPointError(j,0,0,Exp[j]*(ggFSigma[j])-twosiglow[j]*(ggFSigma[j]),twosighigh[j]*(ggFSigma[j])-Exp[j]*(ggFSigma[j]));
    }
    
    multigraph->Add(grae,"3");
    
    grae = new TGraphAsymmErrors(22);
    grae->SetName("sig1");
    grae->SetTitle("Expected CLs #pm 1 #sigma");
    
    ci = TColor::GetColor("#00ff00");
    grae->SetFillColor(ci);
    
    for ( Int_t j=0; j<Npoints_lm; j++ ) {
      grae->SetPoint(j,mass[j],Exp[j]*(ggFSigma[j]));
      grae->SetPointError(j,0,0,Exp[j]*(ggFSigma[j])-onesiglow[j]*(ggFSigma[j]),onesighigh[j]*(ggFSigma[j])-Exp[j]*(ggFSigma[j]));
    }
    
    multigraph->Add(grae,"3");
    
    TGraph *graph = new TGraph(Npoints_lm);
    graph->SetName("");
    graph->SetTitle("Expected CLs - Median");
    graph->SetFillColor(1);
    graph->SetLineColor(1);
    graph->SetLineStyle(2);
    graph->SetLineWidth(2);
    
    for (Int_t j=0; j<Npoints_lm; j++) {
      graph->SetPoint(j,mass[j],Exp[j]*(ggFSigma[j]));
    }
    
    multigraph->Add(graph,"L");
    
    
    TGraph *graphN = new TGraph(Npoints_lm);
    graphN->SetName("SM_crosssections");
    graphN->SetTitle("Total SM Cross-Sections");
    graphN->SetFillColor(1);
    graphN->SetLineColor(2);
    graphN->SetLineStyle(2);
    graphN->SetLineWidth(2);
    
    for ( Int_t j=0; j<Npoints_lm; j++) {
      graphN->SetPoint(j,mass[j],(ggFSigma[j]));
    }
    // multigraph->Add(graphN,"L");
    
    
    multigraph->Draw("");
    
    
    gre = new TGraphErrors(Npoints_lm);
    gre->SetName("CLs_observed");
    gre->SetTitle("Observed CLs");
    gre->SetFillColor(1);
    gre->SetLineColor(1);
    gre->SetLineWidth(2);
    gre->SetMarkerStyle(20);
        
    for ( Int_t j=0; j<Npoints_lm; j++) {
      gre->SetPoint(j,mass[j],obs[j]*(ggFSigma[j]));
      gre->SetPointError(j,0,0);
    }
    
 
    
    gre->Draw("l");
    
    
  }
}