// ----------------------------------------------------------------------------------- //
/*
  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 small simple function:
double sqr(double a) {return a*a;}
// ----------------------------------------------------------------------------------- //



// ----------------------------------------------------------------------------------- //
void RollingBall_Fys2Lab_FitData() {
// ----------------------------------------------------------------------------------- //

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


  // Set constants and define variables and histograms:
  int NeventsMax = 56001;
  int Nevents = 0;
  double t, vol;

  int lowhigh = 0;            // To begin with, the voltage (vol) is low (i.e. no ball passing!)
  double limit_vol  = 2.0;    // Define SOME LEVEL (your choice!) at which "the ball passes".

  double t_min = 0.0;
  double t_max = 0.000025 * double(NeventsMax-1);
  TH1F* Hist_vol = new TH1F("Hist_vol", "Hist_vol;Voltage (arbitrary scale);Number of entries", 6400, 0.0, 4.0);
  TH1F* Hist_tvol = new TH1F("Hist_tvol", "Hist_tvol;time (s);Voltage (arbitrary scale)", NeventsMax, t_min-0.0000125, t_max+0.0000125);


//---------------------------------------------------------------------------------- 
  // File to read:
//---------------------------------------------------------------------------------- 

  ifstream infile("data_RollingBall_V1.txt");
  if (infile.is_open()) {
    while ((!infile.eof()) && (Nevents < NeventsMax)) {
      infile >> t >> vol;
      if (infile.eof()) break;
      Hist_vol->Fill(vol);                        // Voltage measured.
      Hist_tvol->SetBinContent(Nevents+1, vol);   // Voltage measured vs. time.
      Hist_tvol->SetBinError(Nevents+1, 0.002);   // Voltage uncertainty on that measurement.
                                                  // Note that bin 0 is the "underflow" bin!

      // Register when the voltage rises and falls:
      // ------------------------------------------
      // Going from low -> high:
      if ((vol > limit_vol) && (lowhigh == 0)) {
	lowhigh = 1;   // OK, now we are in a "high voltage" interval!
	printf("  Pass low->high:  %8.6f  %7.4f \n", t, vol);
      }

      // Going from low -> high:
      if ((vol < limit_vol) && (lowhigh == 1)) {
	lowhigh = 0;   // OK, now we are (back) in a "low voltage" interval!
	printf("  Pass high->low:  %8.6f  %7.4f \n", t, vol);
      }

      Hist_vol->Fill(vol);

      // Write out the first 10 measurements:
      if (Nevents < 10) printf("  %2d. events read: %9.6f  %12.8f \n", Nevents, t, vol);
      Nevents++;
    }
  }
  infile.close();
  printf("  Number of entries in total: %6d \n", Nevents);



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

  // Draw the voltage distribution (to see the precision and to define "low" and "high" voltage):
  TCanvas* c0 = new TCanvas("c0", "", 50, 20, 600, 400);

  Hist_vol->GetXaxis()->SetRangeUser(0.225, 0.245);   // Zoom in on the "low voltage" part.
  // Hist_vol->GetXaxis()->SetRangeUser(3.60,3.65);
  TF1 *fitGauss = new TF1("fitGauss", "gaus", 0.23, 0.24);   // Gaussian fit (NOTE: range).
  fitGauss->SetLineColor(4);
  fitGauss->SetLineWidth(2);
  Hist_vol->Fit("fitGauss", "RL");
  Hist_vol->Draw();

  c0->Update();


  //---------------------------------------------------------------------------------- 

  // Drawing time series:
  TCanvas* c1 = new TCanvas("c1", "", 100, 50, 1200, 500);
  Hist_tvol->Draw("e");
  c1->Update();
  // c1->SaveAs("Fig_RollingBall_Fys2Lab_TimeSeries.png");



  //---------------------------------------------------------------------------------- 
  // Your analysis here:
  //---------------------------------------------------------------------------------- 

  double g = 9.8;   // Get your own number!
  double stat_g = 0.1;
  double syst_g = 0.1;

  printf("  Result: g = %6.4f +- %6.4f (stat) +- %6.4f (syst) m/s^-2 \n", g, stat_g, syst_g);
}



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

Consideration:
--------------
Little intro:
   Get a feel for the precision on the vol-measurement.
   (They are in fact discrete with 1600 possible values per unit [0,1]).
   Measure the uncertainty by plotting the distribution of "low" (or "high") level
   measurements (i.e. before/between peaks), and look at the RMS. Is the guess good?
   (It turns out that the precision is better at low level, than at high!)

Main analysis:
   Consider the time series - where would you define the voltage at which the ball passes?
   Does your final result depend on this definition? That would be an obvious systematic
   uncertainty!

Suggestion:
   Since this program provides the crossing times for the ball, and also how they
   change with your definition of crossing, it might be a good idea to include the
   subsequent calculations/analysis in the program as well (as has already been
   indicated).



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