// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for ChiSquare fitting points with errors to a user defined function.

  This ROOT macro is a simple example of how one can fit points in ROOT and evaluate
  the fit using the ChiSquare method and test.

  References:
    Barlow: Chapter 6
    Cowan: Chapter 7
    Bevington: Chapter 9

  Author: Troels C. Petersen (NBI)
  Email:  petersen@nbi.dk
  Date:   12th of September 2011
*/
// ----------------------------------------------------------------------------------- //


// ----------------------------------------------------------------------------------- //
double sqr(double a) {
// ----------------------------------------------------------------------------------- //
  return a*a;
}


// ----------------------------------------------------------------------------------- //
double func_pol1(double *x, double *par) {
// ----------------------------------------------------------------------------------- //
  return par[0] + par[1]*x[0];
}

// ----------------------------------------------------------------------------------- //
double func_pol2(double *x, double *par) {
// ----------------------------------------------------------------------------------- //
  return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}


// ----------------------------------------------------------------------------------- //
void FitPoints() {
// ----------------------------------------------------------------------------------- //
  gROOT->Reset();

  // Setting of general plotting style:
  gStyle->SetCanvasColor(0);
  gStyle->SetFillColor(0);
  // Setting what to be shown in statistics box:
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);


  // Random numbers from the Mersenne-Twister:
  TRandom3 r;
  r.SetSeed(0);  // Initializing the random numbers!


  // ------------------------------------------------------------------ //
  // Define data:
  // ------------------------------------------------------------------ //

  // In this exercise, there are two choices of data:
  //  - Fixed data points (written into the program below)
  //  - Randomly generated data points (generated below).
  // To choose between the two, simply move the "comment-out" markers.

  // Fixed data points (could also be read from file):
  const int Npoints = 9;
  double x[Npoints]  = {1.05, 2.36, 2.94, 4.45, 5.21, 6.19, 6.50, 8.09, 9.42};
  double ex[Npoints] = {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00};
  double y[Npoints]  = {1.73, 3.13, 4.13, 4.93, 5.23, 5.13, 4.63, 4.03, 3.03};
  double ey[Npoints] = {0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13};

  /*
  // Generate the data on the spot using random numbers:
  const int Npoints = 9;
  double x[Npoints];
  double ex[Npoints];
  double y[Npoints];
  double ey[Npoints];

  for (int i=0; i < Npoints; i++) {
    x[i]  = double(i+1);
    ex[i] = 0.0;
    y[i]  = 0.2 * sqr(x[i]) - 1.6 * x[i] + 5.0 + r.Gaus(0.0,1.0);
    ey[i] = 1.0;
  }
  */


  // ------------------------------------------------------------------ //
  // Fit the data and plot the result:
  // ------------------------------------------------------------------ //

  // Make a white canvas:
  c1 = new TCanvas("c1","",650,20,600,450);
  c1->SetFillColor(0);

  // Define a histogram (range, axis titles, etc.), where the graph is to be plottet:
  TH1F *h1 = c1->DrawFrame(0.0,0.0,10.0,10.0);
  h1->SetXTitle("Something integer [in units of your imagination]");
  h1->SetTitleOffset(1.0,"y");
  h1->SetYTitle("Something continuous [perhaps]");

  // Make graphs (using a known number of points in four arrays):
  TGraphErrors* graph = new TGraphErrors(Npoints, x, y, ex, ey);

  // Fit graphs:
  int Npar = 2;
  TF1 *fit = new TF1("fit", "pol1", 0.5, 9.5);
  // TF1 *fit = new TF1("fit", "[0] + [1]*x", 0.5, 9.5);     // Another way of doing the fit!
  // TF1 *fit = new TF1("fit", func_pol1, 0.5, 9.5, Npar);   // Yet, another way...
  fit->SetParameters(0.0, 0.0, 0.0);
  fit->SetLineColor(4);

  graph->SetMarkerColor(1);
  graph->SetMarkerStyle(20);
  graph->SetMarkerSize(1.0);
  graph->Fit("fit","r");
  graph->Draw("P same");

  // Add legend:
  TLegend* legend = new TLegend(0.16,0.76,0.40,0.86);
  legend->AddEntry(graph," Data points ","p");
  legend->AddEntry(fit,  " Fit to data ","l");
  legend->SetFillColor(0);
  legend->Draw();

  // Save to file:
  c1->Update();
  // c1->SaveAs("GraphFitted.png");
}

//---------------------------------------------------------------------------------- 
/*

First of all note, that now we are not fitting a histogram (TH1), but a graph with
errors on the points (TGraphErrors). With graphs it is easier to make non-equidistant
points and they are not the result of filling something into a histogram!


Questions:
----------
 1) Start by using the numbers given, and fit these with a line. What is the Chi2
    and the number of degrees of freedom (Ndof)? What is the probability of this?

 2) Try to fit with a parabola instead. What is now the Chi2 and Ndof? And what is
    the probability of getting a fit that has a HIGHER Chi2, that is a worse fit?
    Is this fit then reasonable?
    How much do you need to inflate the errors to get a 50% chance? Before changing
    the errors, make an estimate yourself!

 3) Now switch to the randomly generated numbers (move the comment-out blocks).
    How good a fit do you get here? Is this to be expected? Try running the fit a
    few times (with different seeds for the random numbers), and see if you get
    what you expect. In fact, what do you expect?

 4) If you were to make a loop over the fit, and plot the probability of the Chi2
    for each fit, which distribution of probabilities would you expect.


Advanced questions:
-------------------
 1) Returning to the fixed numbers, what is the probability that the curve goes
    through origo? How to evaluate this?


*/
//----------------------------------------------------------------------------------