// ----------------------------------------------------------------------------------- //
/*
  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 (by Ulam Stanislav and 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:   23rd of September 2012
*/
// ----------------------------------------------------------------------------------- //

#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("e");
  gStyle->SetOptFit(1111);

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

  // Make histograms:
  // Note how the choices 120
  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 (generate numbers according to PDF f(x) = 2x in [0,1]):
//----------------------------------------------------------------------------------

  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):
    double 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":
    double x_hitmiss;
    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 Transformation numbers:
  TF1 *fit_Trans = new TF1("fit_Trans","[0] + [1]*x", 0.0, 1.0);
  fit_Trans->SetLineColor(4);
  Hist_Trans->Fit("fit_Trans", "RL");
  Hist_Trans->Draw();
  gPad->Update();
  TPaveStats *st1 = (TPaveStats*)Hist_Trans->GetListOfFunctions()->FindObject("stats");
  st1->SetX1NDC(0.12); st1->SetX2NDC(0.40);
  st1->SetY1NDC(0.89); st1->SetY2NDC(0.70);
  st1->SetLineColor(4);

  // Fitting histogram with Hit/Miss numbers:
  TF1 *fit_HitMiss = new TF1("fit_HitMiss","[0] + [1]*x", 0.0, 1.0);
  fit_HitMiss->SetLineColor(4);
  Hist_HitMiss->Fit("fit_HitMiss", "RL");
  Hist_HitMiss->Draw();
  gPad->Update();
  TPaveStats *st2 = (TPaveStats*)Hist_HitMiss->GetListOfFunctions()->FindObject("stats");
  st2->SetX1NDC(0.12); st2->SetX2NDC(0.40);
  st2->SetY1NDC(0.69); st2->SetY2NDC(0.50);
  st2->SetLineColor(2);


  // 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_HitMiss, "Hit & Miss", "LE");
  legend->AddEntry(Hist_Trans, "Transformation", "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.
    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.

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