// C. Delaere, 02/10/2012

#include "TRandom.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TStopwatch.h"
#include <iostream>

void genRandom() {

  // various initialization values
  float x_max = 15.;
  int nbins = 150;
  int Nevents = 1000000;
  
  // create histograms
  TH1F* histoRejectionExponential = new TH1F("histoRejectionExponential", "histoRejectionExponential", nbins, 0, x_max);
  TH1F* histoInversionExponential = new TH1F("histoInversionExponential", "histoInversionExponential", nbins, 0, x_max);
  TH1F* histoRootExponential      = new TH1F("histoRootExponential",      "histoRootExponential",      nbins, 0, x_max);

  // timing service
  TStopwatch *st=new TStopwatch();

  // loop over N events (if it is too slow use smaller number of events)
  // *** Rejection method ***
  double random_x_value = 0;
  st->Start();
  for(int i=0; i<Nevents; i++){
    do {
      random_x_value = gRandom->Uniform(0, x_max);
    } while(gRandom->Uniform(0,1) > TMath::Exp(-random_x_value));
    histoRejectionExponential->Fill(random_x_value);
  }
  st->Stop();
  std::cout << "Time for rejection implementation: " << st->CpuTime() << std::endl;
  // *** Inversion method ***
  st->Start();
  for(int i=0; i<Nevents; i++){
    histoInversionExponential->Fill( - TMath::Log( gRandom->Uniform(0,1) ) );
  }
  st->Stop();
  std::cout << "Time for inversion implementation: " << st->CpuTime() << std::endl;
  // *** Simply from ROOT ***
  st->Start();
  for(int i=0; i<Nevents; i++){
    histoRootExponential->Fill( gRandom->Exp(1) );
  }
  st->Stop();
  std::cout << "Time for ROOT implementation: " << st->CpuTime() << std::endl;

  // draw the histograms
  TCanvas *cv1 = new TCanvas("Exponential");
  cv1->Divide(3,1);
  cv1->cd(1);
  histoRejectionExponential->Draw();
  gPad->SetLogy();
  cv1->cd(2);
  histoInversionExponential->Draw();
  gPad->SetLogy();
  cv1->cd(3);
  histoRootExponential->Draw();
  gPad->SetLogy();

  // compare with each other
  std::cout << "compatibility between the Exponential distributions: " << std::endl;
  std::cout << "ROOT vs Inversion: " << histoRootExponential->KolmogorovTest(histoInversionExponential) << std::endl;
  std::cout << "ROOT vs Rejection: " << histoRootExponential->KolmogorovTest(histoRejectionExponential) << std::endl;
  std::cout << "Rejection vs Inversion: " << histoRejectionExponential->KolmogorovTest(histoInversionExponential) << std::endl;
  TCanvas *cv2 = new TCanvas("Differences");
  cv2->Divide(2,1);
  cv2->cd(1);
  TH1F* histoExponentialDiff1 = (TH1F*)histoRootExponential->Clone("histoExponentialDiff1");
  histoExponentialDiff1->Add(histoInversionExponential,-1 * histoRootExponential->GetEntries()/histoInversionExponential->GetEntries());
  histoExponentialDiff1->Divide(histoRootExponential);
  histoExponentialDiff1->SetTitle("Difference between ROOT and Inversion");
  histoExponentialDiff1->Draw();
  cv2->cd(2);
  TH1F* histoExponentialDiff2 = (TH1F*)histoRootExponential->Clone("histoExponentialDiff2");
  histoExponentialDiff2->Add(histoRejectionExponential,-1 * histoRootExponential->GetEntries()/histoRejectionExponential->GetEntries());
  histoExponentialDiff2->Divide(histoRootExponential);
  histoExponentialDiff2->SetTitle("Difference between ROOT and Rejection");
  histoExponentialDiff2->Draw();
}