#include <algorithm>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>

#include "TString.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TPad.h"
#include "TAxis.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TProfile.h"
#include "TF1.h"
#include "TPaveStats.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TVectorD.h"
#include "TPad.h"
#include "TSpline.h"


void savePlot(TCanvas*, TString);
void ReproduceFitProblem() {

  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(   102210);
  gStyle->SetOptFit(1111);
  gStyle->SetOptTitle(0);
  gStyle->SetGridStyle(3);
  gROOT->UseCurrentStyle();


  std::vector<double> bestgainFS, bestgainHS, bestgerrFS, bestgerrHS;
  int nfields = 11;
  
  bestgainFS.push_back(1.92111); bestgerrFS.push_back(0.22271);
  bestgainFS.push_back(2.25945); bestgerrFS.push_back(0.117518);
  bestgainFS.push_back(4.96651); bestgerrFS.push_back(0.0969638);
  bestgainFS.push_back(13.8233); bestgerrFS.push_back(0.202495);
  bestgainFS.push_back(41.0857); bestgerrFS.push_back(0.52347);
  bestgainFS.push_back(128.55);  bestgerrFS.push_back(1.55124);
  bestgainFS.push_back(414.757); bestgerrFS.push_back(4.92111);
  bestgainFS.push_back(1376.4);  bestgerrFS.push_back(16.2247);
  bestgainFS.push_back(4463.23); bestgerrFS.push_back(74.5165);
  bestgainFS.push_back(15342.7); bestgerrFS.push_back(671.713);
  bestgainFS.push_back(56983.1); bestgerrFS.push_back(8321.13);

  bestgainHS.push_back(0.962373); bestgerrHS.push_back(0.140619);
  bestgainHS.push_back(1.05265);  bestgerrHS.push_back(0.869236);
  bestgainHS.push_back(1.62597);  bestgerrHS.push_back(0.167776);
  bestgainHS.push_back(2.53181);  bestgerrHS.push_back(0.087133);
  bestgainHS.push_back(4.11346);  bestgerrHS.push_back(0.0686325);
  bestgainHS.push_back(6.59672);  bestgerrHS.push_back(0.0873496);
  bestgainHS.push_back(11.2233);  bestgerrHS.push_back(0.131859);
  bestgainHS.push_back(18.9957);  bestgerrHS.push_back(0.21812);
  bestgainHS.push_back(33.6369);  bestgerrHS.push_back(0.378424);
  bestgainHS.push_back(58.4167);  bestgerrHS.push_back(0.652233);
  bestgainHS.push_back(101.658);  bestgerrHS.push_back(1.46026);
  
  TGraphErrors * gr_noerr = new TGraphErrors(nfields, &bestgainHS[0], &bestgainFS[0], 0, 0);
  TF1 * tf_noerr          = new TF1("tf_noerr", "pow(x,[0])", bestgainHS[2], bestgainHS[nfields-1]);
  gr_noerr->SetMarkerColor(kBlack); gr_noerr->SetMarkerStyle(20); gr_noerr->SetMarkerSize(1.5);
  tf_noerr->SetLineColor(kRed);     tf_noerr->SetLineStyle(1);    tf_noerr->SetLineWidth(2);
  tf_noerr->SetParameter(1,2.5);    tf_noerr->SetParName(0,"a");
  gr_noerr->Fit(tf_noerr, "BRM+");


  TGraphErrors * gr_err = new TGraphErrors(nfields, &bestgainHS[0], &bestgainFS[0], &bestgerrHS[0], &bestgerrFS[0]);
  TF1 * tf_err          = new TF1("tf_err", "pow(x,[0])", bestgainHS[2], bestgainHS[nfields-1]);
  gr_err->SetMarkerColor(kBlack); gr_err->SetMarkerStyle(20); gr_err->SetMarkerSize(1.5);
  tf_err->SetLineColor(kRed);     tf_err->SetLineStyle(2);    tf_err->SetLineWidth(2);
  tf_err->SetParameter(1,2.5);    tf_err->SetParName(0,"b");
  gr_err->Fit(tf_err, "BRM+");

  TCanvas *c0 = new TCanvas("c0","c0",0,0,600,600);
  c0->SetFillStyle(4000);    c0->SetFillColor(4000); c0->SetFrameFillColor(4000); c0->SetFrameFillStyle(4000);
  c0->SetFrameBorderMode(0); c0->SetTicks(1,1);      c0->SetGrid(1,1);            c0->SetLogx(1); c0->SetLogy(1);
  c0->cd();
  gPad->SetLeftMargin(0.125); gPad->SetBottomMargin(0.125);
  gr_noerr->Draw("AP");
  gr_err->Draw("PEsames");
  gr_noerr->GetXaxis()->SetTitle("x");  gr_noerr->GetXaxis()->SetTitleOffset(1.25);
  gr_noerr->GetYaxis()->SetTitle("y");  gr_noerr->GetYaxis()->SetTitleOffset(1.25);
  gr_noerr->GetYaxis()->SetRangeUser(0.1,2E5);
  gr_noerr->GetXaxis()->SetMoreLogLabels();
  gPad->Update(); // !!! FUNDAMENTAL for picking up statbox !!! 

  TLatex latex_fit;   latex_fit.SetNDC();   latex_fit.SetTextSize(0.025);    latex_fit.SetTextAlign(11);    latex_fit.SetTextColor(kRed);
  
  TPaveStats *st0A = (TPaveStats*)(gr_noerr->GetListOfFunctions()->FindObject("stats"));
  double newx1 = 0.15, newx2 = 0.40, newy1 = 0.76, newy2 = 0.86;
  st0A->SetX1NDC(newx1);    st0A->SetX2NDC(newx2);    st0A->SetY1NDC(newy1); st0A->SetY2NDC(newy2);
  st0A->SetTextColor(kRed); st0A->SetLineColor(kRed); st0A->SetLineStyle(1); st0A->SetLineWidth(2);
  std::stringstream st0Astream; st0Astream<<"y = x^{a} [w/o Errors - OK]"; std::string st0Astring = st0Astream.str();
  latex_fit.DrawLatex(0.15,0.87,st0Astring.c_str());
  st0A->Draw();

  TPaveStats *st0B = (TPaveStats*) gr_err->GetListOfFunctions()->FindObject("stats");
  double newx3 = 0.15, newx4 = 0.40, newy3 = 0.62, newy4 = 0.72;
  st0B->SetX1NDC(newx3);    st0B->SetX2NDC(newx4);    st0B->SetY1NDC(newy3); st0B->SetY2NDC(newy4);
  st0B->SetTextColor(kRed); st0B->SetLineColor(kRed); st0B->SetLineStyle(2); st0B->SetLineWidth(2);
  std::stringstream st0Bstream; st0Bstream<<"y = x^{b} [w/  Errors - PROBLEM]"; std::string st0Bstring = st0Bstream.str();
  latex_fit.DrawLatex(0.15,0.73,st0Bstring.c_str());
  st0B->Draw();
 
  c0->Update();
  savePlot(c0, "ReproduceFitProblem"); 
  // ---------------------------------------------------------------------------------------------------------


  
}

void savePlot(TCanvas * c, TString name) {
  c->SaveAs(name+".png");
  c->SaveAs(name+".eps");
  gSystem->Exec("epstopdf "+name+".eps");
  std::cout<<"Plot "<<name<<" saved"<<std::endl;
}