// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for illustrating the two methods for generating random numbers according
  to a known PDF:
   - Transformation method (if function can be integrated and then inverted).
   - Hit & Miss (or Accept-Reject) method (invented by John Von Neumann).

  The function used for illustration is: f(x) = 2x, in the interval [0,1].

  For more information see:
    G. Cowan: Chapter 3
    P. R. Bevington: page 81-84
    
  Author: Troels C. Petersen (NBI/CERN)
  Email:  Troels.Petersen@cern.ch
  Date:   16th of September 2010
*/
// ----------------------------------------------------------------------------------- //

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



//---------------------------------------------------------------------------------- 
void MakingRandomNumbers()
//---------------------------------------------------------------------------------- 
{
  // 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);
  gStyle->SetStatX(0.48);
  gStyle->SetStatY(0.88);

  TRandom3 r;
 
 // Set parameters:
  int Nnumbers = 10000;       // Number of random numbers produced.

  double x_trans, x_hitmiss;

  // Make histograms:
  TH1F *Hist_Trans   = new TH1F("Hist_Trans",   "Hist_Trans",   120, -0.1, 1.1);
  TH1F *Hist_HitMiss = new TH1F("Hist_HitMiss", "Hist_HitMiss", 120, -0.1, 1.1);


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

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

    // Transformation method:
    // ----------------------
    // Integration gives the function F(x) = x^2, which inverted gives F-1(x) = sqrt(x):
    x_trans = sqrt(r.Uniform());
    Hist_Trans->Fill(x_trans);

    // Hit & Miss method:
    // ------------------
    // Generate two random numbers uniformly distributed in [0,1]x[0,2], until they
    // fulfill the "Hit requirement":
    do {
      x_hitmiss = 1.0 * r.Uniform();
      double y  = 2.0 * r.Uniform();
    } while (2.0*x_hitmiss < y);
    Hist_HitMiss->Fill(x_hitmiss);
  }


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

  TCanvas* c0 = new TCanvas("c0", "", 100, 20, 600, 400);

  // Setting line width and color and axis titles of histograms:
  Hist_Trans->GetXaxis()->SetTitle("Random number");
  Hist_Trans->GetYaxis()->SetTitle("Frequency");
  Hist_Trans->SetLineWidth(2);
  Hist_Trans->SetLineColor(2);
  Hist_HitMiss->SetLineWidth(2);
  Hist_HitMiss->SetLineColor(4);

  // Fitting histogram with HitMissson:
  TF1 *fit_Trans = new TF1("fit_Trans","[0] + [1]*x", 0.0, 1.0);
  fit_Trans->SetLineColor(4);
  Hist_Trans->Fit("fit_Trans", "R");

  TF1 *fit_HitMiss = new TF1("fit_HitMiss","[0] + [1]*x", 0.0, 1.0);
  fit_HitMiss->SetLineColor(4);
  Hist_HitMiss->Fit("fit_HitMiss", "R");


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

  // Legend:
  TLegend* legend = new TLegend(0.60, 0.18, 0.85, 0.33);
  legend->SetLineColor(1);
  legend->SetFillColor(0);
  legend->AddEntry(Hist_Trans, "Transformation", "LE");
  legend->AddEntry(Hist_HitMiss, "Hit & Miss", "LE");
  legend->Draw();

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


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

  Hist_Trans->Write();
  Hist_HitMiss->Write();

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



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


Questions:
----------
 1) Make sure that you understand how the two methods works! What is the
    efficiency of each of the methods? 

 2) Compare the two histograms - are they from the same distribution? Answer by
    first comparing the fit parameters, and see if they overlap. Afterwards,
    calculate the Chi2 between the two histograms (i.e. the difference divided
    by the uncertainty bin by bin), and compare this to the number of degrees
    of freedom.
    To compute the chance of observing this Chi2 given Ndof, use the function
    "Prob" in ROOT's TMath (http://root.cern.ch/root/html/TMath.html) as follows:
    double prob = TMath::Prob(Chi2,Ndof).

 3) Try to repeat the exercise with the function f(x) = x^3 in the range [1,C],
    (perhaps in a new file, so that you save this macro), where C is a number,
    which ensures that the PDF f(x) is normalized. Think about what changes.

 4) Again compare the two histograms, this time perhaps using the fact, that
    ROOT has a function for TH1, which gives the Chi2 to another TH1 (Chi2Test).


Advanced questions:
-------------------
 1) Try also to apply the run test (Barlow 8.3.2) also known as the Wald-Wolfowitz
    runs test. If this is succesful, then try to shift one distribution by a little
    and see if the runs test is more sensitive than the Chi2 test.

 2) ROOT technical question:
    With two histograms in one plot, how can you make it print two statistics boxes?


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