#include "TMath.h"
#include "TRandom3.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFitResult.h"

#include <iostream>
#include <fstream>
#include <cmath>
#include <cfloat>

void runtestfit(void){

  // Example commands to fit a line to data.
  // Once beyond Phys1001, your data points' errors vary from point to point
  // So how do you fit data then?
  // The example uses two artificial data sets designed ot have proper statistical properties
  // likely want to read data in from a file instead.
  
  // Force some plot style parameters for all plots
  // can override these as necessary on a per-plot basis
  gROOT->SetStyle("Plain");
  gStyle->SetLineWidth(3);
  gStyle->SetHistLineWidth(3);
  gStyle->SetLabelSize(0.05,"xyz");
  gStyle->SetTitleSize(0.05,"xyz");
  gStyle->SetOptStat(kFALSE);
  //gStyle->SetOptTitle(kFALSE);
  gStyle->SetOptTitle(kTRUE);
  gStyle->SetTitleBorderSize(-1);
  gStyle->SetPalette(1,0);  //(50, 0) for deep blue sea
  gROOT->ForceStyle();


  double xdata[6] = {1.0, 2.0, 3.0, 5.0, 7.0, 10.0};
  double  ydata[6] = {1.64, 3.22, 3.87, 6.85, 10.96, 14.37};
  double  sigma[6] = {0.2, 0.4, 0.5, 0.6, 0.7, 0.8};
  double zeros[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  // Alternate, less precise data gives larger uncertainties on slope
  //double ydata[6] = {1.43, 2.10, 3.20, 6.43, 12.33, 12.75};
  //double sigma[6] = {0.40, 0.80, 1.00, 1.20, 1.40, 1.60};


  TGraphErrors *tg = new TGraphErrors(6, xdata, ydata, zeros, sigma);
  tg->SetMarkerStyle(21);
  tg->GetXaxis()->SetTitle("time (milliseconds)");
  tg->GetYaxis()->SetTitle("distance (micrometers)");
  tg->SetTitle("illustration of weighted fit");
  // I'm surprised the title is left-justified and in a box.
  // Need to fix that.

  tg->Draw("AP");

  // now fit a first-order polynomial to the data with error bars.
  TF1 *fitfn = new TF1("fit1","pol1");  // a math function, to be fit.  1st order polynomial

  //TFitResultPtr myfit =  tg->Fit("pol1","S");   // a simpler one-line syntax
   
  TFitResultPtr myfit =  tg->Fit(fitfn,"S");
  // TFitResultPtr myfit =  tg->Fit(fit1,"WS"); // gives unweighted simple fit.

  // print the result.
  std::cout << " " << std::endl; 
  std::cout << "====================================" << std::endl;
  std::cout << "Approximate fit and one-sigma errors" << std::endl;
  std::cout << "from diagonal of covariance matrix" << std::endl;
  myfit->Print("V");

  // Root's default behavior is to treat them as error bars, not just weights
  // But warning, other statistical fit packages (R, scipy, Matlab)
  // treat them only as weights describing the relative importance of each point
  // Requiring you to scale the parameter errors yourself.
  
  // would be cool to draw parameters with higher and lower values.
  // there is still a funny thing here.  The two fit parameters are correlated
  // some uses of fits to data need to incorporate that information too.

}