// C. Delaere, 06/09/2010
// from a code byAlexander.Schmidt@cern.ch 03/20/2010

#include "TRandom.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TMath.h"

void genRandom() {
 
  // select a range on the x axis
  double x_min_max = 2.5;

  // create histogram
  TH1F * histoBruteForceGaussian = new TH1F("histoBruteForceGaussian",
  "histoBruteForceGaussian", 100, -10, 10);
  TH1F * histoRootGaussian = new TH1F("histoRootGaussian",
  "histoRootGaussian", 100, -10, 10);
  TH1F * histoInversionExponential = new TH1F("histoInversionExponential",
  "histoInversionExponential", 150, 0,15);
  TH1F * histoRootExponential = new TH1F("histoRootExponential",
  "histoRootExponential", 150, 0,15);

  // loop over N events (if it is too slow use smaller number of events)
  int Nevents = 10000000;
  for(int i=0; i<Nevents; i++){
    // get a random x value between -x_min_max and +x_min_max
    double random_x_value = gRandom->Uniform(-x_min_max, x_min_max);
    // decide if the x value should be filled into the histogram or NOT !
    // this needs to be done by checking the y value of the
    // normalized gauss function
    // draw another random number and fill the x-value into histogram
    if( gRandom->Uniform(0,1) < TMath::Gaus(random_x_value, 0 , 1, kTRUE)){
      histoBruteForceGaussian->Fill(random_x_value);
    }
    histoRootGaussian->Fill(gRandom->Gaus(0,1));
    // the second part of the exercise: exponential distribution
    // simply get it from root
    histoRootExponential->Fill( gRandom->Exp(1) );
    // or do it yourself with inversion
    histoInversionExponential->Fill( - TMath::Log( gRandom->Uniform(0,1) ) );
  }

  // draw the histograms, compare with each other
  TCanvas *cv1 = new TCanvas("Gaussian");
  histoBruteForceGaussian->Draw();
  histoRootGaussian->Draw("same");
  gPad->SetLogy();
  histoRootGaussian->GetYaxis()->SetRangeUser(1e-1,1e6);
  std::cout << "compatibility between the two Gaussian distributions: " << histoRootGaussian->KolmogorovTest(histoBruteForceGaussian) << std::endl;

  TCanvas *cv2 = new TCanvas("Exponential");
  histoInversionExponential->Draw();
  histoRootExponential->Draw("same");
  gPad->SetLogy();
  std::cout << "compatibility between the two Exponential distributions: " << histoRootExponential->KolmogorovTest(histoInversionExponential) << std::endl;

  TCanvas *cv2 = new TCanvas("Differences");
  cv2->Divide(2,1);
  cv2->cd(1);
  TH1F* histoGaussianDiff = histoRootGaussian->Clone("histoGaussianDiff");
  histoGaussianDiff->Add(histoBruteForceGaussian,-1 * histoRootGaussian->GetEntries()/histoBruteForceGaussian->GetEntries());
  histoGaussianDiff->Divide(histoRootGaussian);
  histoGaussianDiff->Draw();
  cv2->cd(2);
  TH1F* histoExponentialDiff = histoRootExponential->Clone("histoExponentialDiff");
  histoExponentialDiff->Add(histoInversionExponential,-1 * histoRootExponential->GetEntries()/histoInversionExponential->GetEntries());
  histoExponentialDiff->Divide(histoRootExponential);
  histoExponentialDiff->Draw();

}