// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for exercise on calibration.

  You're developping a new method for measuring the distance to stars from Earth based
  telescopes, and want to calibrate this method. You know that the method depends on
  several factors, such as:
   * Amount of signal light from the star
   * Amount of background light in the surrounding sky
   * Temperature of star
   * Transparency of sky

  In order to determine the influence of these factors, and how much you need to correct
  for each of them, you consider 10.000 stars with known distances. From these, you can
  find how well your own method works, make corrections to biases as needed, and finally
  find out how precise your calibrated method is. Happy calibration.


  Author: Troels C. Petersen (NBI)
  Email:  petersen@nbi.dk
  Date:   28th of October 2012
*/
// ----------------------------------------------------------------------------------- //

#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include "TMath.h"


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


// ----------------------------------------------------------------------------------- //
Calibration() {
// ----------------------------------------------------------------------------------- //
  gROOT->Reset();

  gStyle->SetOptStat(1111);
  gStyle->SetOptFit(1111);
  gStyle->SetCanvasColor(0);
  gStyle->SetPalette(1);

  bool SaveFigures = false;



  // ------------------------------------------------------------------ //
  // Make Histograms and vectors:
  // ------------------------------------------------------------------ //

  TH1F* Hist_Raw   = new TH1F("Hist_Raw",   "Hist_Raw",   80, -2.0, 2.0);
  TH1F* Hist_Calib = new TH1F("Hist_Calib", "Hist_Calib", 80, -2.0, 2.0);

  // Example histograms of "Luminosity Signal (lsig)" distribution
  TH1F* Hist_lsig = new TH1F("Hist_lsig", "Hist_lsig", 50, 0.0, 50.0);
  TH2F* Hist_lsig2D = new TH2F("Hist_lsig2D", "Hist_lsig2D", 50, 0.0, 50.0,  50, -1.0, 1.0);


  // ------------------------------------------------------------------ //
  // Loop over data and make plots:
  // ------------------------------------------------------------------ //

  int Nstars = 0;
  double dist_truth, dist_obs, lsig, lbkg, temp, tsky;

  ifstream datafile;
  datafile.open("data_calib.txt");
  if (datafile.is_open()) {
    while (!datafile.eof()) {
      datafile >> dist_truth >> dist_obs >> lsig >> lbkg >> temp >> tsky;
      if (Nstars < 10) printf("  %9.4f  %9.4f      %9.4f  %9.4f  %9.4f  %9.4f \n",
    	    dist_truth, dist_obs, lsig, lbkg, temp, tsky);

      // Note the initial relative resolution:
      double ddistrel = (dist_obs - dist_truth) / dist_truth;
      Hist_Raw->Fill(ddistrel);

      // Do calibration (expand this yourself!):
      double dist_obs_calib = dist_obs;    // Make corrections that gives a better resolution...

      // Note the final relative resolution:
      double ddistrel_calib = (dist_obs_calib - dist_truth) / dist_truth;
      Hist_Calib->Fill(ddistrel_calib);


      // Histograms for making calibration:
      Hist_lsig->Fill(lsig);
      Hist_lsig2D->Fill(lsig, ddistrel);

      Nstars++;
    }
  }
  datafile.close();


  printf("\n  The initial and final resolutions are:  %6.3f  and  %6.3f \n",
	 Hist_Raw->GetRMS(), Hist_Calib->GetRMS());


  // ------------------------------------------------------------------ //
  // Show the distribution of fitting results:
  // ------------------------------------------------------------------ //

  TCanvas* c0 = new TCanvas("c0","", 20, 50, 600, 400);
    
  Hist_Raw->SetLineWidth(2);
  Hist_Raw->SetLineColor(2);
  Hist_Raw->Draw();
  Hist_Calib->SetLineWidth(2);
  Hist_Calib->SetLineColor(4);
  Hist_Calib->Draw("same");

  c0->Update();
  if (SaveFigures) c0->SaveAs("UncalibratedCalibrated.png");


  // ------------------------------------------------------------------ //
  // Show the distribution of fitting results:
  // ------------------------------------------------------------------ //

  TCanvas* c1 = new TCanvas("c1","", 120, 100, 600, 400);
    
  Hist_lsig2D->GetXaxis()->SetRangeUser(0.0, 50.0);
  Hist_lsig2D->Draw();
  TProfile* ProfX_lsig = (TProfile*)Hist_lsig2D->ProfileX("ProfX_lsig");
  ProfX_lsig->SetLineWidth(2);
  ProfX_lsig->SetLineColor(2);
  ProfX_lsig->Draw("same");

  c1->Update();

}


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

As always look at the data and get a feel for each of the variables. A good idea might
be to plot them all to know what range to expect them in.

Next, consider the distribution of relative differences between the observed (by your
method!) and actual distance: (dist_obs - dist_truth) / dist_truth
The RMS is 0.27, i.e. a 27% precision, and it is neither centered at zero (as it
should not to be biassed), nor very Gaussian. This is what you want to improve upon.

Finally, there is the distribution of the bias and relative precision as a function of
the signal luminosity (lsig). As you can see, the response does not depend on lsig,
and so there seems to be no (varying) bias from this variable. Check the other three
variables, and if you find some bias, try to correct for it, such that you get the
most narrow Gaussian centered around 0.


Questions:
----------
 1) What corrections do you apply for each of the variables, and how much do each of
    them improve on the result?
 2) Try at the end to figure out, what the final resolution depends on. Can you give
    and estimate of the uncertainty for each single star?

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