// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for reading text (ASCII) files and fitting data written.

  Author: Troels C. Petersen (NBI/CERN)
  Email:  Troels.Petersen@cern.ch
  Date:   8th of September 2012
*/
// ----------------------------------------------------------------------------------- //

#include <TFile.h>
#include <fstream>
#include <iostream>

using namespace std;


// ----------------------------------------------------------------------------------- //
// Define a function if a variable(s) "x" and parameter(s) "par" (standard way):
double func_OscExp(double *x, double *par) {
// ----------------------------------------------------------------------------------- //
  return par[0] + par[1] * exp(par[2]*x[0]) * cos(par[3] + par[4]*x[0]);
}


// ----------------------------------------------------------------------------------- //
void HarmonicOscillator_Fys2Lab_FitData() {
// ----------------------------------------------------------------------------------- //

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

  int NeventsMax = 3617;           // OK, I read this before writing the macro!
  int Nevents = 0;
  double t, x;

  // Defining the histogram to contain the events. It could equally well have been a TGraph,
  // and if the time measurements were not equidistant, this would have been the way to go!
  // (In a TGraph, one specifies x_i, error_x_i, y_i, and error_y_i for N points.)
  double t_min = 0.0;
  double t_max = 0.01 * double(NeventsMax-1);
  TH1F* Hist_Osc = new TH1F("Hist_Osc", "Hist_Osc", NeventsMax, t_min-0.005, t_max+0.005);

  // File to read and fill into histogram:
  ifstream infile("data_HarmonicOscillator_V1.txt");
  if (infile.is_open()) {
    while ((!infile.eof()) && (Nevents < NeventsMax)) {
      infile >> t >> x;
      if (infile.eof()) break;
      Hist_Osc->SetBinContent(Nevents+1, x);       // Value measured at time t
      Hist_Osc->SetBinError(Nevents+1, 0.005);     // Error ESTIMATED for this value
                                                   // (probably it is closer to 0.0035)!

      // Write out the first 10 measurements (always a good idea!):
      if (Nevents < 10) printf("  %2d. events read: %6.3f  %6.3f \n", Nevents, t, x);
      Nevents++;
    }
  }
  infile.close();
  printf("  Number of entries in total: %6d \n", Nevents);


//---------------------------------------------------------------------------------- 
// Fit and Plot:
//---------------------------------------------------------------------------------- 

  TCanvas* c0 = new TCanvas("c0", "", 50, 20, 1200, 500);

  // Fitting histogram with oscillating exponential:
  // The first name should match that in front, while the second "func_OscExp"
  // refers to the function defined at the top. The function is defined in the
  // range [0.005; 36.005], and it has 5 paramters (this last information is
  // required, when the function is defined elsewhere - ROOT wants to know!).
  TF1 *fitOscExp = new TF1("fitOscExp", func_OscExp, 0.005, 36.005, 5);

  // For a complex function of many parameters like this, we are forced to give
  // it "reasonable" starting parameters (which is not always easy!).
  fitOscExp->SetParameters(0.037, 0.20, -0.1, -1.0, 9.0);

  // To fix a parameter (make function simpler, i.e. lower dimensionality), do:
  // fitOscExp->FixParameter(3, 0.0);
  fitOscExp->SetLineWidth(2);     // This is just linewidth
  fitOscExp->SetNpx(1000);        // This is how finely the line should be plottet

  // Setting range, line width and color and axis titles of histograms:
  Hist_Osc->GetXaxis()->SetRangeUser(-0.005, 36.165);
  Hist_Osc->GetXaxis()->SetTitle("Time (seconds)");
  Hist_Osc->GetYaxis()->SetTitle("Oscillation (arbitrary scale)");
  Hist_Osc->SetLineWidth(2);
  Hist_Osc->SetLineColor(4);

  // Fitting and drawing histograms with errors:
  Hist_Osc->Fit("fitOscExp", "R");
  Hist_Osc->Draw("e");

  c0->Update();
  // c0->SaveAs("AppliedStatistics_HarmonicOscillator_OscExpFit.png");
}



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


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