// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for illustrating error propagation using random Gaussian numbers.

  For more information on error propagation, see:
    R. J. Barlow: page 48-61 
    P. R. Bevington: page 36-48

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

#include <TFile.h>
#include <TH1F.h>
#include <TRandom3.h>


//---------------------------------------------------------------------------------- 
// Defining a new function (here squaring!), which as input takes a double (called a),
// and returns a new double (a*a). The function "sqr" can then be used througout
// the remains of the program/macro.
double sqr(double a) {
  return a*a;
}


//---------------------------------------------------------------------------------- 
void ErrorPropagation()
//---------------------------------------------------------------------------------- 
{
  // Setting of general plotting style:
  gStyle->SetCanvasColor(0);
  gStyle->SetFillColor(0);
  // Setting what to be shown in statistics box:
  gStyle->SetOptStat(1110);   // Print: Entries, Mean, and RMS
  gStyle->SetOptFit(1111);    // Print: Everything!


  // Set parameters:
  int Nexp = 10000;
  double pi = 3.14159265358979323846264338328;
  TRandom3 r;

  // Make histograms:
  TH1F *Hist_x1 = new TH1F( "Hist_x1", "Hist_x1", 100,  -5.0,  5.0);
  TH1F *Hist_x2 = new TH1F( "Hist_x2", "Hist_x2", 100,  -5.0,  5.0);
  TH1F *Hist_y  = new TH1F( "Hist_y",  "Hist_y",  200, -10.0, 10.0);


//----------------------------------------------------------------------------------
// Loop over process:
//----------------------------------------------------------------------------------

  for (int i=0; i < Nexp; i++) {

    // Get random Gaussian numbers:
    double x1 = r.Gaus( 3.0, 1.5);
    Hist_x1->Fill(x1);
    double x2 = r.Gaus(-0.8, 0.2);
    Hist_x2->Fill(x1);

    // Calculate y:
    double y = x1 + 5.0 * x2;
    Hist_y->Fill(y);

    if (i < 5) printf("  Gaussian: %5.2f %5.2f  ->  %5.2f \n", x1, x2, y);
  }


//---------------------------------------------------------------------------------- 
// Plot histograms on screen:
//---------------------------------------------------------------------------------- 

  // Making a new window (canvas):
  // --------------------------------------------------------------------
  TCanvas* c0 = new TCanvas("c0", "", 100, 20, 600, 400);
  // c0->SetLogy();

  // Setting line width and color of histograms:
  Hist_y->SetLineWidth(2);
  Hist_y->SetLineColor(4);

  // Fitting histogram with Gaussian:
  // TF1 *fitGauss = new TF1("fitGauss","gaus", -8.0, 8.0);
  // Hist_y->Fit("fitGauss", "RL");

  // Drawing histograms with errors in the same plot:
  Hist_y->Draw("e");

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



//---------------------------------------------------------------------------------- 
// Write histograms to file:
//---------------------------------------------------------------------------------- 
  /*
  TFile *file = new TFile("results.root","RECREATE","Histograms");

  Hist_y->Write();

  file->Write();
  file->Close();
  */
}



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

Before running the program, consider our variable of interest, y, as a function of
x1 and x2 in three cases:
  a) y = x1 + 5*x2
  b) y = x1 * sqr(x2)
  c) y = log(sqr(x1*tan(x2))+sqrt(abs(x1-x2)/(cos(x2)+1.0+x1)))

In each of the three cases, x1 and x2 has been measured to be (with Gaussian errors):
     x1 =  3.0 +- 1.5
     x2 = -0.8 +- 0.2


Questions:
----------
 1) Calculate (i.e. analytically!) what the central value and uncertainty on y will
    be for case a) and b) using the error propagation formula.

 2) Next, convince yourself that if you choose x1 and x2 according to the distributions
    above and combine them following equation a), then that should tell you what the
    central value and uncertainty of y will be. Run the program, and see if it does!

 3) Repeat this for case b).

 4) In example c), try to see if you can quickly differentiate this with respect to
    x1 and x2, and thus predict what uncertainty to expect for y using the error
    propagation formula. It is OK to give up :-)
    Next, try to estimate it using random numbers like above - do you trust this
    number?


Advanced questions:
-------------------
 1) Repeat question 1) and 2), but this time given that x1 and x2 has a correlation
    of 0.6.

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