// ----------------------------------------------------------------------------------- //
/*
  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;
}


  // SolutionLevel3:
  // To finalize the result, I will calculate the weighted means, and compare them to
  // the unweighted. Though some measurements are discarted (either for being flawed, or
  // simply not have an error), the expectation is that the uncertainty on the mean
  // will be (even) smaller.


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

  // 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 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];

  // SolutionLevel3: Now also include weighted results (only two sums, but also a measurement counter):
  double wnmeas[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double wsum0[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double wsum1[Ndata] = {0.0, 0.0, 0.0, 0.0};
  double wmean[Ndata];
  double wemean[Ndata];

  double approx_mean[Ndata] = {433.8, 377.9, 336.0, 336.0};
  char* name_data[Ndata] = {"ByEye", "ByHuman", "By30cmR", "By2mFR"};

  // 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];

      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);

      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 (fabs(l[3] - 336.45) < 0.06 && el[3] > 0.0 && el[3] < 0.15)
	el[3] = 2.0*el[3]; // Inflate too small errors (in measurements 336.4 and 336.5).


      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]);
	  }

	  // SolutionLevel3: For the weighted mean, measurements are excluded based on their pull.
	  if (el[i] > 0.0 && fabs(l[i] - approx_mean[i]) / el[i] < 6.0) {
	    wnmeas[i] += 1.0;
	    wsum0[i]  += 1.0  / sqr(el[i]);
	    wsum1[i]  += l[i] / sqr(el[i]);
	  }
	  Hist_Meas[i]->Fill(l[i]);
	  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:
  // SolutionLevel3: Loop over sums to get means and widths (include another decimal!):
  for (int i=0; i < Ndata; i++) {
    if (sum0[i] > 1.0 && wsum0[i] > 0.0) {

      // Unweighted:
      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.   Nmeas = %2.0f   Mean = %7.2f +- %7.2f   Width = %7.2f    \n",
	     i, sum0[i], mean[i], emean[i], width[i]);

      // Weighted:
      wmean[i] = wsum1[i] / wsum0[i];
      wemean[i] = sqrt(1.0 / wsum0[i]);
      printf("            Nmeas = %2.0f   Mean = %7.2f +- %7.2f \n",
	     wnmeas[i], wmean[i], wemean[i]);
    } else {
      printf("  Error: Too little data to calculate mean and width! %1d. Nevents = %3d \n", i, int(sum0[i]));
    }
  }



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

  // Make canvas:
  // SolutionLevel3: 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!


  // SolutionLevel3: 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(wemean[2]) +      1.0/sqr(wemean[3]); 
  double wmean1 = wmean[2]/sqr(wemean[2]) + wmean[3]/sqr(wemean[3]); 

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

}


/*
SolutionLevel3: Discussion after analysis step
----------------------------------------------
  After including the weighted averages, we conclude that:
  1: The two estimates with reasonable errors improve their precision quite significantly, as the
     uncertainty is lowered by a factor of two (corresponding to 4 times as many unweighted
     measurements!!!). Here, a weighted average seems justified, if some measurements are excluded.

  2: For the 30cm measurements, the weighted mean does not improve the precision, while it does so
     for the 2m measurements (again almost a factor of two). However, the 2m weighted average is
     dominated by two measurements, which claim to have a precision of 0.1cm:
       336.4 +- 0.1   and   336.5 +- 0.1
     Is this really trustworthy? I somehow doubt it, especially since they are slightly in conflict
     with the mean from all the other measurements (we got 335.984 +- 0.093 unweighted).

  And so the analysis ends at a splitting point.
  Either, these two points are trusted, in which case the averages obtained are as follows:
    Data 0.   Nmeas = 43   Mean =  433.81 +-   14.58   Width =   95.60    
              Nmeas = 41   Mean =  416.68 +-    6.55 
    Data 1.   Nmeas = 43   Mean =  377.88 +-    4.95   Width =   32.46    
              Nmeas = 42   Mean =  379.01 +-    2.57 
    Data 2.   Nmeas = 35   Mean =  336.03 +-    0.19   Width =    1.13    
              Nmeas = 33   Mean =  336.16 +-    0.19 
    Data 3.   Nmeas = 36   Mean =  335.97 +-    0.11   Width =    0.64    
              Nmeas = 31   Mean =  336.38 +-    0.06 
  Note the shift in central value for the 2m weighted mean, as expected.

  Or, they are overestimates of how precisely this can be done, and they are excluded or perhaps
  included but with inflated uncertainties. The next highest precision is 0.2cm, which seems more
  reasonable, and this would give:
    Data 0.   Nmeas = 43   Mean =  433.81 +-   14.58   Width =   95.60    
              Nmeas = 41   Mean =  416.68 +-    6.55 
    Data 1.   Nmeas = 43   Mean =  377.88 +-    4.95   Width =   32.46    
              Nmeas = 42   Mean =  379.01 +-    2.57 
    Data 2.   Nmeas = 35   Mean =  336.03 +-    0.19   Width =    1.13    
              Nmeas = 33   Mean =  336.16 +-    0.19 
    Data 3.   Nmeas = 36   Mean =  335.97 +-    0.11   Width =    0.64    
              Nmeas = 31   Mean =  336.29 +-    0.10 
  Now the shift in central value is less, but notice that now the improvement by doing a weighted
  average is almost gone! The weighted average between the two is now: 336.259 +- 0.087 cm

  Having tried the weighted average and seen the results, I would conclude that it does not
  really improve the precision, as the uncertainties cannot be fully trusted.
  I would also report these difficulties in combining the data, so that the reader can know,
  what was tried and with which result.

*/


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

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!)


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