#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

void BrazilPlot() {
  /*
    Script to make limits plots from text file with format:
    mH observed expected exp+2sigma exp+1sigma exp-1sigma exp-2sigma
    R.Goncalo 1/9/2011
  */  
  cout << "BrazilPlot: draw upper limit +- uncertainty from file" << endl;

  cout << "--- declaring internal variables" << endl;
  float mass[999], obs[999], exp[999], expU2[999], expU1[999];
  float expD1[999], expD2[999], zero[999];
  for (unsigned int m = 0; m<999; ++m) {
    mass[m]=0; 
    obs[m]=0; 
    exp[m]=0; 
    expU2[m]=0; 
    expU1[m]=0; 
    expD1[m]=0; 
    expD2[m]=0;
    zero[m]=0;
  }

  
  cout << "--- Input file : " ;
  TString file;
  cin >> file;

  string line;
  ifstream infile (file);

  int mass_points=0;

  if (infile.is_open()) {

    cout << "--- reading in variables from input file" << endl;
    unsigned int i=0;
    vector<float> variables;
    variables.resize(7, 0);

    while ( !infile.eof() ) {
      infile >> variables[i];
      i++;
      // reset i, increment mass points and move variables
      if ( i == 7 ) {
	for ( int j=0; j<7; ++j ) std::cout << variables[j] << " ";
	std::cout << std::endl;
	i=0;
	mass[mass_points]  = variables[0]; 
	obs[mass_points]   = variables[1]; 
	exp[mass_points]   = variables[2]; 
	expU2[mass_points] = fabs(variables[3] - exp[mass_points]);  
	expU1[mass_points] = fabs(variables[4] - exp[mass_points]);  
	expD1[mass_points] = fabs(exp[mass_points] - variables[5]); 
	expD2[mass_points] = fabs(exp[mass_points] - variables[6]); 
	mass_points++;
      }
    }
    // close file when done
    infile.close();
  } else {
    cout << "Unable to open file"; 
  }

  cout << "--- check inputs:" << endl;
  cout << "m(H) obs    exp   exp+2s exp+1s exp-1s exp-2s" << endl;
  for (int k=0; k < mass_points; ++k) {
    cout << mass[k] << " " << obs[k] << " " << exp[k] << " " 
	 << expU2[k] << " " << expU1[k] << " " << expD1[k] << " " << expD2[k] 
	 << endl;
  }

  cout << "--- opening output root file" << endl;
  TFile* outfile = new TFile("BrazilPlot.root", "recreate");

  cout << "--- creating graphs" << endl;
  TGraph* gObs = new TGraph(mass_points, mass, obs);
  TGraph* gExp = new TGraph(mass_points, mass, exp);
  TGraphAsymmErrors* gExp1s = new TGraphAsymmErrors(mass_points, mass, exp, zero, zero, expD1, expU1);
  TGraphAsymmErrors* gExp2s = new TGraphAsymmErrors(mass_points, mass, exp, zero, zero, expD2, expU2);

  cout << "--- making plot canvas" << endl;
  TCanvas* c1 = new TCanvas("c1","c1",800,600);
  c1->cd();
  //c1->SetLogy();
  gExp2s->GetYaxis()->SetTitle("95% C.L. limit on #sigma/#sigma_{SM}");
  gExp2s->GetXaxis()->SetTitle("m_{H} [GeV]");
  gExp2s->SetTitle("");
  gExp2s->SetFillColor(5);
  gExp2s->Draw("A3");
  gExp1s->SetFillColor(3);
  gExp->SetLineColor(1);
  gExp->SetLineStyle(2);
  gExp->SetLineWidth(4);
  gExp1s->Draw("3");
  gExp->SetLineColor(1); 
  gExp->Draw("LP");
  gObs->SetLineWidth(4);
  gObs->SetLineColor(2); 
  gObs->Draw("LP");
  c1->Print("limits.png");

  cout << "--- closing output file BrazilPlot.root" << endl;
  outfile->Close();
  
  cout << "--- done!" << endl;
}