#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
using namespace RooFit ;

double B_value(double x);
double B_dependance(double x);

void example()
{

   Double_t xAxis1[6] = {-4, -3.7, -3.3, -2.85, -2.4, -2}; 

   TH1D *hist_to_fit = new TH1D("hist_to_fit","Unfold Response",5, xAxis1);
   hist_to_fit->SetBinContent(1,503295.4);
   hist_to_fit->SetBinContent(2,539456.5);
   hist_to_fit->SetBinContent(3,549759.1);
   hist_to_fit->SetBinContent(4,499685);
   hist_to_fit->SetBinContent(5,478172.9);
   hist_to_fit->SetBinError(1,2208.979);
   hist_to_fit->SetBinError(2,1709.717);
   hist_to_fit->SetBinError(3,1434.279);
   hist_to_fit->SetBinError(4,1382.245);
   hist_to_fit->SetBinError(5,1460.313);
   hist_to_fit->SetMinimum(0);
   hist_to_fit->SetMaximum(773938.2);
   hist_to_fit->SetEntries(10);

   // Declare observable
   RooRealVar logXi("logXi","logXi",-4, -2) ;

   // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
   RooDataHist dh("dh","dh",logXi,Import(*hist_to_fit)) ;


   // Make plot of binned dataset
   RooPlot* frame = logXi.frame(Title("Imported TH1")) ;
   dh.plotOn(frame) ; 

   // Declare variables & set initial values
   RooRealVar C("C","C",0,pow(10,8)) ;
   RooRealVar D("D","D",-1.,1.) ;

   C.setVal(1.5*pow(10,5));
   D.setVal(-0.21);

   // Declare pdf
   RooGenericPdf genpdf("genpdf","genpdf","C * B_dependance(logXi)*pow(10, D*logXi)",RooArgSet(logXi,C, D)) ;

   // Fit pdf to data
   genpdf.fitTo(dh,SumW2Error(kTRUE));

   // // Make plot frame in observable, plot fitted pdf
   RooPlot* frame1 = logXi.frame(Title("genpdf")) ;
   genpdf.plotOn(frame) ;

   // Draw on canvas
  TCanvas* myCanvas = new TCanvas("myCanvas","myCanvas",800,400) ;
  myCanvas->Divide(2);
  myCanvas->cd(1);
  gPad->SetLeftMargin(0.15);
  frame->GetYaxis()->SetTitleOffset(1.6);
  frame->Draw();

  
  // Plot result of fit in TF1 & compare it to original plot

  myCanvas->cd(2);
  hist_to_fit->SetMarkerStyle(20);
  TF1* function = new TF1("function","[0] * B_dependance(x)*pow(10, [1]*x)",-4, -2.);
  function->SetParameter(0, C.getVal());
  function->SetParameter(1, D.getVal());
  function->Draw("");
  hist_to_fit->Draw("E1 same");

  std::cout<<"\n\n"<<C.getVal()<<"\t"<<D.getVal()<<"\n";
}



double B_dependance(double x){
   double t_high = -0.06;
   double t_low = -0.9;

   double B_dependance = (1. / B_value(x)) * (exp(t_high*B_value(x)) - exp(t_low*B_value(x)));
   return B_dependance;
}

double B_value(double x){
   // double B_measured = 7.6;
   double B_measured = 5.09;
   double alphaPrime = 0.25;
   double aver_log10Xi = -2.98;

   double B0 = B_measured + (2 * alphaPrime * log(pow(10, aver_log10Xi)));

   return B0 - (2 * 0.25 * log(pow(10, x)));
}