#include "Riostream.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TPostScript.h"
#include "TGraphAsymmErrors.h"
#include "TLine.h"
#include "TFile.h"
#include "TLatex.h"

void plot_ww_xsec() {

  // Plot WW cross-section
  // ---------------------
  // root[] gSystem->CompileMacro("plot_ww_xsec.C");
  // root[] plot_ww_xsec()
  //
  // See options below.
  // This macro read in data file "ww_xsec.dat"
  //
  // Macro originally developed by Jianrong Deng, August 2004
	
  gROOT->Reset();
  gROOT->SetStyle("Plain");

  gStyle->SetLineWidth(3);

  ifstream in;
  // we assume a file x_sect.dat in the current directory
  in.open("ww_xsec.dat");
  // number of points to plot
  Int_t npoints = 40;

  // Zoom in on a certain energy range ? [Select logScale for UNZOOMED plot]
  bool zoomPlot = false;
  bool LogScale = true;
  bool includeRun2 = true;

  // legend labels
  const int nCurve = 20;
  // char grLgd[nCurve][20] = {"0", "ppbar CTEQ6L", "ppbar CTEQ6M", "pp CTEQ6L", "pp CTEQ6M", "RunI", "RunII CDF", "RunII D0"};
  char grLgd[nCurve][50] = {"0", "ppbar CTEQ6L (LO)", "ppbar CTEQ6M (NLO)", "pp CTEQ6L (LO)", "pp CTEQ6M (NLO)", "Run I CDF", "Run II CDF DILEPTON", "Run II CDF LEPTON+TRACK   ", "Run II D0"};
  int numCurve = 0;
  Int_t points;
  Int_t n = 40;
  Float_t sqrt_s[n], sigma[n], dsigma[n],dsigmaNeg[n], dsigmaX[n],dsigmaXNeg[n];
  TFile *f = new TFile("x_sect.root","RECREATE");
  TPostScript *ps;
  if ( zoomPlot ) ps = new TPostScript("x_sect_1T2.eps",111);
  else            ps = new TPostScript("x_sect_log.eps",111);
  TCanvas *c1 = new TCanvas("c1","Graph Draw Options",700,700);

  TLegend* leg = NULL;
  if (zoomPlot) {
    leg = new TLegend(0.12,0.55,0.52,0.88);
  } else {
    leg = new TLegend(0.11,0.63,0.45,0.89);
  }

  TGraph *gr[20];
  TGraphAsymmErrors *grData[20];

  // read data from ascii file: x_sect.dat
  int nCurvePlot = nCurve;
  if (zoomPlot) nCurvePlot = 2; // don't want proton curves for zoomed plot
  for (Int_t j=1; j!=nCurvePlot + 1; ++j) 
    {
      in >> points;
      if (in)
	{
	  printf("points=%d\n",points);
	  ++numCurve;
	  for (Int_t i=0; i!=points + 1; ++i) {
	    int tmp;
	    in >> tmp >> sqrt_s[i] >> sigma[i] >> dsigma[i] ;
	    sqrt_s[i] = sqrt_s[i]/1000.0;
	    sigma[i] = sigma[i]/1000.0;
	    dsigma[i] = dsigma[i]/1000.0;
	    if (i < n) printf("point= %d, sqrt_s = %f, sigma= %f, dsigma=%f\n",tmp, sqrt_s[i],sigma[i],dsigma[i]);
	  }   
	  points = points + 1;
	  gr[j] = new TGraph (npoints, sqrt_s, sigma);
	  // gr[j]->GetXaxis()->SetTitle("sqrt(s) (TeV)");
	  gr[j]->GetXaxis()->SetTitle("Center of Mass Energy     (TeV)");
	  // gr[j]->GetXaxis()->SetTitleOffset(1.2);
	  gr[j]->GetXaxis()->CenterTitle();
	  gr[j]->GetYaxis()->SetTitle("WW Cross-Section      (pb)");
	  gr[j]->GetYaxis()->SetTitleOffset(1.4);
	  gr[j]->GetYaxis()->CenterTitle();
	  // gr[j]->SetTitle("WW cross-section (pb)");
	  gr[j]->SetTitle("");
	  printf(" found %d points\n",points);
	}
    }
  
  
  // RunI and RunII data points:
  //Run1: 10.2 +6.3-5.1 (stat) +- 1.6 (syst)  ~
  //10.2 +6.5 -5.34
  //(not sure about the luminosity error in Run1)
  //sqrt(s)=1.8
  //                                                                                                                             
  //Run2 (CDF) blessed: 14.3 +5.6-4.9 (stat) +- 1.6 (syst) +- 0.9 (lum)
  // ~14.3 +5.89 -5.23   sqrt(s)=1.96  
  //Run2 (D0)  blessed: 13.8 +4.3 -3.8 (stat) +1.0 -0.8 (syst) +-0.9 (lum)
  //~13.8 +4.5 -3.98   sqrt(s)=1.96
   
  // RunI
  sqrt_s[0] = 1.8;
  sigma[0] = 10.2;
  dsigma[0] = 6.5;
  dsigmaNeg[0] = 5.34;
  dsigmaX[0] = 0;
  dsigmaXNeg[0] = 0;
  grData[1] = new TGraphAsymmErrors (1, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);

  // RunII CDF (DILEPTON, ERRORS ALREADY ADDED IN QUADRATURE) :
  //sqrt_s[0] = 1.96; separated the CDF and D0 points
  sqrt_s[0] = 2.02;
  if ( zoomPlot ) sqrt_s[0] = 1.95;
  sigma[0] = 14.6;
  dsigma[0] = 6.139;
  dsigmaNeg[0] = 5.98;
  dsigmaX[0] = 0;
  dsigmaXNeg[0] = 0;
  grData[2] = new TGraphAsymmErrors (1, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);

  // RunII CDF (LEPTON+TRACK,ERRORS ALREADY ADDED IN QUADRATURE) :
  //sqrt_s[0] = 1.96;
  sqrt_s[0] = 1.93;
  if ( zoomPlot ) sqrt_s[0] = 1.96;
  sigma[0] = 24.2;
  dsigma[0] = 8.769;
  dsigmaNeg[0] = 9.074;
  dsigmaX[0] = 0;
  dsigmaXNeg[0] = 0;
  grData[3] = new TGraphAsymmErrors (1, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);
  // TGraphAsymmErrors *grData = new TGraphAsymmErrors (3, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);

  // RunII D0
  //sqrt_s[0] = 1.96;
  sqrt_s[0] = 2.13;
  if ( zoomPlot ) sqrt_s[0] = 1.97;
  sigma[0] = 13.8;
  dsigma[0] = 4.55;
  dsigmaNeg[0] = 4.01;
  dsigmaX[0] = 0;
  dsigmaXNeg[0] = 0;
  grData[4] = new TGraphAsymmErrors (1, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);
  //TGraphAsymmErrors *grData = new TGraphAsymmErrors (3, sqrt_s, sigma,dsigmaX,dsigmaXNeg,dsigma,dsigmaNeg);

  c1->Update();
  ps->NewPage();
  c1->cd();
  // set logscale
  // set axis range
  if ( zoomPlot ) 
    {
      gr[numCurve]->SetMinimum(0.3);
      gr[numCurve]->SetMaximum(34);
      if (includeRun2) {
	gr[numCurve]->GetXaxis()->SetRangeUser(1.7,2.05);
      } else {
	gr[numCurve]->GetXaxis()->SetRangeUser(1.7,2.00);
      }
      gr[numCurve]->SetLineStyle(numCurve);
      gr[numCurve]->SetLineWidth(3);
    }
  else
    {
      gr[numCurve]->SetMinimum(0.3);
      gr[numCurve]->SetMaximum(500);
      gr[numCurve]->GetXaxis()->SetRangeUser(1,100);
      gr[numCurve]->SetLineStyle(numCurve);
      gr[numCurve]->SetLineWidth(3);
      if (LogScale) {
	gPad->SetLogx();
	gPad->SetLogy();
      }
    }

  // Ticks on all 4 axes (without labels) :
  gPad->SetTicks(1,1);
  // Ticks on all 4 axes (with labels) :
  // gPad->SetTicks(2,2);
   
  // plot the LO and NLO curves in the same color
  gr[numCurve]->SetLineColor((numCurve+1)/2);
  gr[numCurve]->Draw("AL"); 
  gr[numCurve]->Write();
  leg->AddEntry(gr[numCurve],grLgd[numCurve],"L");
  for (Int_t j=numCurve-1; j!=0; --j) 
    {
      leg->AddEntry(gr[j],grLgd[j],"L");
      gr[j]->SetLineStyle(j);
      gr[j]->SetLineColor((j+1)/2);
      gr[j]->SetLineWidth(3);
      gr[j]->Draw("L"); 
      gr[j]->Write();
    }
  // plot run data points 
  double markerSize = 1.2;
  if (zoomPlot) markerSize = 1.5;
  int legendId = numCurve;
  if (zoomPlot) legendId += 2;   // skip proton curve legends
  grData[1] ->SetMarkerStyle(20);
  grData[1] ->SetMarkerSize(markerSize);
  grData[1] ->SetLineWidth(3);
  leg->AddEntry(grData[1],grLgd[legendId+1],"P");
  grData[1] ->Draw("P");
  if (includeRun2) {
    grData[2] ->SetMarkerStyle(24);
    grData[2] ->SetMarkerSize(markerSize);
    grData[2] ->SetLineWidth(3);
    leg->AddEntry(grData[2],grLgd[legendId+2],"P");
    grData[2] ->Draw("P");
    grData[3] ->SetMarkerStyle(26);
    grData[3] ->SetMarkerSize(markerSize);
    grData[3] ->SetLineWidth(3);
    leg->AddEntry(grData[3],grLgd[legendId+3],"P");
    grData[3] ->Draw("P");
    grData[4] ->SetMarkerStyle(25);
    grData[4] ->SetMarkerSize(markerSize);
    grData[4] ->SetLineWidth(3);
    leg->AddEntry(grData[4],grLgd[legendId+4],"P");
    grData[4] ->Draw("P");
  }
  leg->Draw();


  // c1->SetGridx();
  // c1->SetGridy();

  // Add LHC line :
  TLine* l = new TLine(14.0,0.3,14.0,500.0);
  l->SetLineWidth(1);
  l->SetLineStyle(1);
  l->SetLineColor(2);

  if(!zoomPlot) {l->Draw("same");}
  
  TLatex* tl = new TLatex(13.0,10.0,"LHC");
  tl->SetTextAngle(90.0);
  tl->SetTextSize(0.03);
  tl->SetTextColor(2);

  if(!zoomPlot) {tl->Draw();}


  c1->Update();
  c1->Write();
  ps->Write();
  ps->Close();
  in.close();
  f->Write();
  if(!zoomPlot) {c1->SaveAs("wwxs_tev_lhc.eps");}
  if(zoomPlot) {c1->SaveAs("wwxs_tev.eps");}
  // gROOT->ProcessLine(".q");
}