// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for analysing measurements of the length of the lecture table in Auditorium
  A at NBI.

  There are four measurements of the table length:
    1: Estimated by eye.
    2: Estimated using human scale (Troels).
    3: Measurement with a 30cm ruler.
    4: Measurement with a 2m folding ruler.
  Each person was asked not only to state the measurement, but also their estimated
  uncertainty, though some persons did not include an uncertainty, and in some cases
  measurements are missing (indicated by the value -1). As both sides of the table
  were used for measurements, this was also indicated (0 for room side, 1 for board
  side).
  None of the persons could see others measurements in order to get the largest
  degree of independence.

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

// Include from C++:
#include <fstream>
#include <iostream>

// Include from ROOT:
#include <TROOT.h>
#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TFile.h>


//---------------------------------------------------------------------------------- 
// Simple example of defining your own function. Reading from an end, it returns a
// "double" (i.e. a number with (many) decimals), it's name is "sqr", and as input
// it takes a "double". And between the {}, you define what it should do!
double sqr(double a) {
  return a*a;
}


  // SolutionLevel2:
  // In the second iteration, I correct/exclude some measurements, and choose to look at
  // the errors, to see if they are realistic (and thus usable).
  //
  // The corrections/exclusions can be done either in the program or in the data file.
  // Both are OK, but in the latter case one should DEFINITELY KEEP THE ORIGINAL DATA FILE!!!


// ----------------------------------------------------------------------------------- //
// This is the main function. It does not return anything (void) and takes no input.
// It simply runs whatever is "inside".
void TableMeasurements_SolutionLevel2() {
// ----------------------------------------------------------------------------------- //

  // Set the showing of statistics and fitting results (0 means off, 1111 means all on):
  gROOT->Reset();
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit(1111);

  // Set the graphics:
  gStyle->SetStatBorderSize(1);
  gStyle->SetStatFontSize(0.055);
  gStyle->SetCanvasColor(0);
  gStyle->SetPalette(1);

  int NeventsMax = 100;  // Maximum number of events to be read.
  int Nevents = 0;       // Counter of events read so far.
  double x[9];           // Array of 9 numbers for storing the current measuremets read.

  // Define three sums to be increased with each event for calculation of mean (and error) and width:
  const int Ndata = 4;
  double sum0[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double sum1[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double sum2[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double mean[Ndata];
  double emean[Ndata];
  double width[Ndata];

  // SolutionLevel2: I define approximate means (from SolutionLevel1) to be used for the pull and accuracy:
  //                 They can be updated at any level of the analysis, as they are simply references.
  double approx_mean[Ndata] = {433.8, 377.9, 336.0, 336.0};

  // SolutionLevel2: OK - let us automate the histograms, and make new ones for the resuals/pulls,
  //                 which are defined as (measurement - average) / sigma. These are good for testing
  //                 if the errors are reasonable.
  char* name_data[Ndata] = {"ByEye", "ByHuman", "By30cmR", "By2mFR"};
                     // This defines an array of "char", i.e. several characters in a row, i.e. a name!

  // Histograms (and their specifications) for the measurements:
  int nbins_Meas[Ndata] = {50, 50, 100, 100};
  double min_Meas[Ndata] = {0.0, 0.0, 325.0, 325.0};
  double max_Meas[Ndata] = {1000.0, 1000.0, 350.0, 350.0};
  TH1F* Hist_Meas[Ndata];

  // Histograms for the pull distributions (to see if the errors seem realistic):
  int nbins_Pull[Ndata] = {100, 100, 100, 100};
  double min_Pull[Ndata] = {-10.0, -10.0, -10.0, -10.0};
  double max_Pull[Ndata] = { 10.0,  10.0,  10.0,  10.0};
  TH1F* Hist_Pull[Ndata];

  // Histograms for the accuracy (i.e. distance to mean) compared to the uncertainty (as reported), thus in 2D:
  int nbins_Ac[Ndata] = {150, 150, 100, 100};
  double min_Ac[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double max_Ac[Ndata] = {300.0, 300.0, 10.0, 10.0};
  int nbins_Er[Ndata] = {150, 150, 100, 100};
  double min_Er[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double max_Er[Ndata] = {300.0, 300.0, 10.0, 10.0};
  TH2F* Hist_AcVsEr[Ndata];

  char name[80];     // This is simply for the "output" of the sprintf command below, which "prints" into a string.
  for (int i=0; i < Ndata; i++) {
    sprintf(name, "Hist_Meas_%s", name_data[i]);        // When i=0, name="Hist_Meas_ByEye" and so forth.
    Hist_Meas[i] = new TH1F(name, name, nbins_Meas[i], min_Meas[i], max_Meas[i]);
    sprintf(name, "Hist_Pull_%s", name_data[i]);
    Hist_Pull[i] = new TH1F(name, name, nbins_Pull[i], min_Pull[i], max_Pull[i]);
    sprintf(name, "Hist_AcVsEr_%s", name_data[i]);
    Hist_AcVsEr[i] = new TH2F(name, name, nbins_Ac[i], min_Ac[i], max_Ac[i], nbins_Er[i], min_Er[i], max_Er[i]);
  }

  ifstream dataset("TableMeasurements_data.csv");      // Open the data file.


  // Start the loop over data:
  // -------------------------
  if (dataset.is_open()) {              // If opening the file succeeded, then...
    // While the reading is not at the End-of-File (eof),
    // and the number of events does not exceed the maximum, do...
    while ((!dataset.eof()) && (Nevents < NeventsMax)) {    

      // Read the dataset into the array "x" of length 9.
      dataset >> x[0] >> x[1] >> x[2] >> x[3] >> x[4] >> x[5] >> x[6] >> x[7] >> x[8];

      // SolutionLevel2: As there are 4 measurements and 4 errors, put these into vectors as well
      // so that any analysis step can be applied to all 4 data sets in a loop:
      double l[Ndata]  = {x[0], x[2], x[4], x[6]};
      double el[Ndata] = {x[1], x[3], x[5], x[7]};

      // Write out the first 5 events to the screen with some values.
      // Note how the slighly cryptic "%5.2f" refers to 5 spaces with 2 digits
      // for a float (i.e. double), thus printing out the value (here of x[3]) nicely.
      if (Nevents < 5)
	printf("  Read: %6d  %7.1f   %3.0f %2.0f  %5.2f \n", Nevents, x[0], x[1], x[2], x[3]);

      // For every 10 events (here the "%" is integer division, i.e. 12%5 = 2), print progress.
      if (Nevents%10 == 0)
	printf("  Read %6d events. \n", Nevents);

      // SolutionLevel2: Correct some of the measurements (i.e. l[2] and l[3])!!!
      if (fabs(l[2] - 306.0) < 6.0) l[2] +=  30.0;  // Correct 30cm up
      if (fabs(l[2] - 366.0) < 6.0) l[2] -=  30.0;  // Correct 30cm down
      if (fabs(l[3] - 136.0) < 6.0) l[3] += 200.0;  // Correct 2.0m up
      if (fabs(l[3] - 536.0) < 6.0) l[3] -= 200.0;  // Correct 2.0m down

      // If the value read is good, use this measurement:
      // SolutionLevel2: Loop over measurements and get the sums for each of them (and put them into histograms).
      //                 Note that some measurements are excluded, on the grounds of being obvious mistakes.
      //                 But they are kept in the histograms (for record), where we simply fit in a small range.
      for (int i=0; i < Ndata; i++) {
	if (l[i] > 0.0) {
	  if ((i<2) || ((i>1) && (fabs(l[i] - 336.0) < 6.0))) {
	    sum0[i] += 1.0;
	    sum1[i] += l[i];
	    sum2[i] += sqr(l[i]);
	  }
	  Hist_Meas[i]->Fill(l[i]);

	  // For the pull and Accuracy Vs. Error plots, the error is needed:
	  if (el[i] > 0.0) {
	    Hist_Pull[i]->Fill((l[i] - approx_mean[i]) / el[i]);
	    Hist_AcVsEr[i]->Fill(el[i], abs(l[i] - approx_mean[i]));
	  }
	}
      }

      Nevents++;                  // Increase the number of events read by one.
    }
  }

  // After reading the data, finish:
  dataset.close();                                          // Close the data file.
  printf("  Number of entries in total: %6d \n", Nevents);  // Write the total number of events read.

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

  // Calculate mean and width:
  // SolutionLevel2: Loop over sums to get means, widths, and errors on the means (include another decimal!):
  for (int i=0; i < Ndata; i++) {
    if (sum0[i] > 1.0) {
      mean[i] = sum1[i] / sum0[i];
      width[i] = sqrt(sum2[i] / sum0[i] - sqr(mean[i]));
      emean[i] = width[i]/sqrt(sum0[i]);
      printf("  Data %1d.   Mean = %7.2f +- %7.2f  Width = %7.2f \n", i, mean[i], emean[i], width[i]);
    } else {
      printf("  Error: Too little data to calculate mean and width! %1d Nevents = %3d \n", i, int(sum0[i]));
    }
  }
  // SolutionLevel2: Note how the errors on the measurement means is now getting small!



// ----------------------------------------------------------------------------------- //
  // Plot histogram on the screen:

  // Make canvas:
  // SolutionLevel2: Plot all the four measurement distributions:
  TCanvas* canvas = new TCanvas("canvas","",250,20,800,600);     // Make a new window
  canvas->Divide(2,2);

  TF1* func_gauss[Ndata];
  double min_fit[Ndata] = {   0.0,    0.0, 330.0, 330.0};
  double max_fit[Ndata] = {1000.0, 1000.0, 342.0, 342.0};
  for (int i=0; i < Ndata; i++) {
    canvas->cd(i+1);             // Note that the first canvas is labeled "1" (and not "0" as normally).
    sprintf(name, "func_gauss_%s", name_data[i]);
    func_gauss[i] = new TF1(name, "gaus", min_fit[i], max_fit[i]);  // Define a Gaussian function
    Hist_Meas[i]->Fit(name,"LR");                                        // Fit histogram with function
  }

  canvas->Update();
  // canvas->SaveAs("LengthOfTable.png");   // Alternative formats: jpg, eps, pdf, etc!


  // SolutionLevel2: Also plot the pull distributions and the accuracy vs. error:
  TCanvas* canvas2 = new TCanvas("canvas2","",280,40,800,600);     // Make a new window
  canvas2->Divide(2,2);

  TF1* func_pull[Ndata];
  for (int i=0; i < Ndata; i++) {
    canvas2->cd(i+1);
    sprintf(name, "func_pull_%s", name_data[i]);
    func_pull[i] = new TF1(name, "gaus", -3.0, 3.0);  // Define a Gaussian function
    Hist_Pull[i]->Fit(name,"LR");                     // Fit histogram with function
  }

  canvas2->Update();
  // canvas2->SaveAs("LengthOfTable_pull.png");   // Alternative formats: jpg, eps, pdf, etc!


  TCanvas* canvas3 = new TCanvas("canvas3","",310,60,800,600);     // Make a new window
  canvas3->Divide(2,2);

  TF1* func_acvser[Ndata];
  for (int i=0; i < Ndata; i++) {
    canvas3->cd(i+1);
    Hist_AcVsEr[i]->SetMarkerStyle(20);                    // Set the markers to be 
    Hist_AcVsEr[i]->SetMarkerSize(0.5);                    // Set the marker size.
    Hist_AcVsEr[i]->GetXaxis()->SetTitle("Error (cm)");    // Set the X-axis title.
    Hist_AcVsEr[i]->GetYaxis()->SetTitle("Accuracy = distance to mean (cm)");
    Hist_AcVsEr[i]->Draw();                                // Draw the scatter plot.

    TLine* unit = new TLine(0.0, 0.0, 200.0, 200.0);
    unit->Draw("same");
  }

  canvas3->Update();
  // canvas3->SaveAs("LengthOfTable_AcVsEr.png");   // Alternative formats: jpg, eps, pdf, etc!


  // Printing the correlation coefficients between accuracies and errors:
  printf("  Accuracy Vs. Error correlation:  30cm ruler:  %5.2f \n", Hist_AcVsEr[2]->GetCorrelationFactor());
  printf("  Accuracy Vs. Error correlation:  2m folding:  %5.2f \n", Hist_AcVsEr[3]->GetCorrelationFactor());


  double wmean0 =     1.0/sqr(emean[2]) +     1.0/sqr(emean[3]); 
  double wmean1 = mean[2]/sqr(emean[2]) + mean[3]/sqr(emean[3]); 

  printf("  Final result:     L = %7.3f +- %5.3f \n", wmean1/wmean0, sqrt(1.0/wmean0));

}


/*
SolutionLevel2: Discussion after analysis step
----------------------------------------------
  After second iteration, there are several conclusions:
  1: The two estimates have reasonable errors, with the exception of two measurements by eye, and
     one by human scale. A weighted average seems justified, if these measurements are excluded.

  2: The 30cm ruler measurements have passable errors, but they are generally overestimated, giving
     a small pull, and thus a large peak around 0. Again, two measurements do NOT lie within reasonable
     distance of the mean in terms of errors, and should be excluded.

  3: The 2m folding ruler measurements again have highly overestimated errors, and also some measurements
     are on the high side. Choosing to accept values within 6 sigma in pull (as was done for "by eye"),
     this would exclude 3 measurements.

   The issue with (somtimes highly) overestimated errors is a bit of a problem for a weighted
   average, as the uncertainty on the mean is calculated from the errors themselves. Thus, even
   if they weighted the measurements correctly according to their uncertainty, the error on this
   mean might be too large.
   So one possibility is NOT to do a weighted average for each data point. This gives:

     Data 0.   Nmeas = 43   Mean =  433.81 +- 14.58   Width = 95.60    
     Data 1.   Nmeas = 43   Mean =  377.88 +-  4.95   Width = 32.46    
     Data 2.   Nmeas = 35   Mean =  336.03 +-  0.19   Width =  1.13    
     Data 3.   Nmeas = 36   Mean =  335.97 +-  0.11   Width =  0.64    

   To get an overall estimate, I would conclude that the two estimates are both biased (if we
   accept that the table length is somewhere around 336cm, the two estimates er certainly not in
   agreement with this) and not at all comparing to the measurements in precision, so these will
   not be included in the average.
   The two measurement averages can be combined in a weighted average, which gives the final mean
   of (look in the code how this is calculated):

         Combined: Mean = 335.984 +- 0.093

   Thus, combining all the measurements, we find that we can measure the length of the table in
   Auditorium A with a precision of about 1mm.

   I continue at SolutionLevel3 for the weighted averages...
*/


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

Start by taking a close look at the data, both by inspecting the numbers, and then by
considering the histograms produced by running the macro.


Questions:
----------
First, consider the two estimates, first by eye and then by human scale (i.e. not the
measurements).

 1)  Calculate their mean and widths. Did the "help" of having a human
    scale decrease the width? Assuming that the table is 336cm long, did it also give a
    mean closer to this value?

 2) Do any of the measurements looks wrong/bad/suspecious? Do you find any reason to
    exclude any measurements?

 3) Were the uncertainties that people gave realistic? Try to plot the estimates with
    the mean subtracted and divided by the uncertainty. Such a plot is called a
    "residual" or "pull" plot, and should be a unit Gaussian if the measurements
    are unbiased and with correct errors. Are they?

Now consider the two measurements.

 4) Ask yourself again, if any of the measurements and/or errors looks wrong/bad?
    Would you correct or exclude any of the measurements and how would you justify this?
    If so, make sure that your program does this, or make a copy of the data file, and
    make the corrections here.

 5) Given the (possibly) updated dataset, fit the length measurements with a Gaussian
    distribution, both for the 30cm and the 2m ruler cases. Should the binning be changed
    for these fits? And is the Gaussian distribution justified?

 6) Given two succesful fits, the width of these Gaussians actually show the "correct"
    uncertainty in the measurements. How do they compare to the uncertainties estimated
    by the people making the measurements?
    How much was the estimated uncertainties lowered from first to second type of
    measurement? Was this justified?

 7) Is there any correlation between the errors on the measurements and the actual
    distance to the mean? That is, do the errors contain useful information, that
    should be used in the calculation of the most accurate estimate of the mean? 

 8) Try to calculate the mean and then the weighted mean including uncertainties.
    How much difference does this make?


Advanced questions:
-------------------
 1) The length of the table seems to be different when measured with a 30cm and a 2m
    ruler - quantify this statement! I.e. what is the difference, and what is the
    uncertainty on that difference? Is this significantly away from 0? How big an
    effect on the 30cm ruler would this correspond to?
    And, is there any sign of a difference in length between the two sides of the table?

 2) In question 3) you chose a binning - does this have any influence on the result?

 3) If you were asked for the best estimate of the length of the table, what would
    you do? (If posssible, read Bevington page 58!)


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