// Solution to Computing and Statistics Problem Sheet 6
// Glen Cowan, RHUL, Physics, November 2009

#include <iostream>
#include <cmath>
#include <TH1D.h>
#include <TFile.h>
#include <TRandom3.h>

using namespace std;

int main(){

// Open output file (apparently needs to be done before booking)

  TFile* file = new TFile("simpleMC.root", "recreate");

// Book histograms

  TH1D* h_Uni = new TH1D("h_Uni", "uniform random numbers",  100,  0, 1.0);
  TH1D* h_Exp = new TH1D("h_Exp", "exponential random numbers",  100,  0, 5.0);
  TH1D* h_xa = new TH1D("h_xa", "r1 + r2 - 1", 100, -1., 1.);
  TH1D* h_xb = new TH1D("h_xb", "r1+r2+r3+r4-2", 100, -2., 2.);
  TH1D* h_xc = new TH1D("h_xc", "sum ri - 6", 100, -6., 6.);
  TH1D* h_trans = new TH1D("h_trans", "x", 100, 0., 1.);
  TH1D* h_accrej = new TH1D("h_accrej", "x", 100, 0., 1.);

// Create a TRandom3 object to generate random numbers

  int seed = 12345;
  TRandom3* ran = new TRandom3(seed);

// Generate some random numbers and fill histograms

  const int numValues = 100000;
  const double xi = 1.0;                // mean of exponential pdf

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

    double r = ran->Rndm();             // uniform in ]0,1]
    double x = - xi * log(r);

    double xa = ran->Rndm() + ran->Rndm() - 1.;

    double xb = 0.;
    for (int i=0; i<4; i++){
      xb += ran->Rndm(); 
    }
    xb -= 2.;
   
    double xc = 0.;
    for (int i=0; i<12; i++){
      xc += ran->Rndm(); 
    }
    xc -= 6.;

    h_Uni->Fill(r);
    h_Exp->Fill(x);
    h_xa->Fill(xa);
    h_xb->Fill(xb);    
    h_xc->Fill(xc);

  }

// Generate random values according to f(x) = 2x (0 < x < 1)
// with transformation method; here x = sqrt(r).

  for (int i=0; i<numValues; i++){
    double r = ran->Rndm();
    double x = sqrt(r);
    h_trans->Fill(x);
  }

// and again using acceptance rejection method

  const double xMin = 0.;
  const double xMax = 1.;
  const double fMax = 2.;
  int numAccepted = 0;
  while ( numAccepted < numValues ){
    double x = (xMax - xMin)*ran->Rndm() + xMin;
    double u = fMax * ran->Rndm();
    double f = 2.*x/(xMax*xMax);
    if ( u < f ){
      numAccepted++;
      h_accrej->Fill(x);
    }
  }

// Store all histograms in the output file and close up

  file->Write();
  file->Close();

  return 0;
}