#include <TFile.h>
#include <TH1F.h>
#include <TCanvas.h>

#include <iostream>

#include <RooRealVar.h>
#include <RooArgSet.h>
#include <RooFitResult.h>
#include <RooPlot.h>
#include <RooDataHist.h>
#include <RooPolynomial.h>
#include <RooChi2Var.h>
#include <RooMinuit.h>

using std::cout;
using std::endl;


void run(){
  TFile fin("test_INPUT.root");
  TH1F *hRatio = dynamic_cast<TH1F*>(fin.Get("ratio"));
  
  TFile fout("test_OUTPUT.root","RECREATE");

  
  RooRealVar x("x","x",150,450.);
  RooDataHist ratio("Ratio", "Ratio",RooArgSet(x), hRatio);

  RooRealVar a("a","a",2.,-10.,10.); 
  RooPolynomial lin("lin", "linear", x, RooArgList(a));
  


  // ****  FIT Possibilities ****
  //RooChi2Var chi2("chi2","chi2",lin,ratio,RooFit::DataError(RooAbsData::SumW2)) ;
  //RooMinuit m(chi2) ;
  //m.migrad() ;
  //m.hesse() ;
  //RooFitResult* sigma1FitResult = m.save();

  // OR
  //RooFitResult *sigma1FitResult = lin.chi2FitTo(ratio, RooFit::Extended(kFALSE), RooFit::Save(kTRUE),RooFit::SumW2Error(kTRUE));

  // OR
  RooFitResult *sigma1FitResult = lin.fitTo(ratio, RooFit::Extended(kFALSE), RooFit::Save(kTRUE),RooFit::SumW2Error(kTRUE), RooFit::NumCPU(4), RooFit::Warnings(kTRUE), RooFit::PrintLevel(0), RooFit::Verbose(kTRUE));
  // ****


  // Plot&Save
  RooPlot *dataFrame;
  dataFrame = x.frame(RooFit::Title("x"));
  
  ratio.plotOn(dataFrame);
  lin.plotOn(dataFrame);
  lin.paramOn(dataFrame, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.1, 0.4,0.9));

  cout<<"Status of the fit: " << sigma1FitResult->status()<<endl;
  sigma1FitResult->Print();
    
  // Save the results in a root file
  TCanvas* c = new TCanvas("r","r", 800, 400);
  c->cd();
  dataFrame->GetYaxis()->SetTitleOffset(1.2); 
  dataFrame->GetYaxis()->SetLabelSize(0.04);
  dataFrame->Draw(); 
  
  fout.cd();
  c->Write();
  //sigma1FitResult->Write();


  // ===> Check what plain root gives:
  hRatio->Fit("pol1");
  hRatio->Write();


  fout.Close();fin.Close();
  c->Close();
}