//
// Part of this is inspired from code of Vittorio Boccone 05/04/07
// C. Delaere, 05/09/2010
//

//TODO: do the tests for all 3 generators

// (Pseudo)Random numbers in Root
#include "TRandom.h"
#include "TRandom3.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TTree.h"
#include "TCanvas.h"

// The power of c/c++... Randu in 1 line! (1977 primitive computers)
// n_i=65539*n_(i-1) mod 2**31
int Randu() { static int n=1; return n = ( n + (n<<1) + (n<<16) ) & 0x7fffffff; }

bool RandomNumber()
{
  const int Events = 1e4;

  // I initialize TRandom3 and TRandom
  TRandom3 *myRandom3 = new TRandom3();
  myRandom3->SetSeed(0);
  TRandom *myRandom = new TRandom();
  myRandom->SetSeed(0);

  // Initialize 3 canvas one for each (pseudo)randomizer
  TCanvas *c1 = new TCanvas("c1","RANDU");
  c1->Divide(2,2);
  TCanvas *c2 = new TCanvas("c2","TRandom");
  c2->Divide(2,2);
  TCanvas *c3 = new TCanvas("c3","TRandom3");
  c3->Divide(2,2);

  // Create 3 trees (usefull to plot 1D,2D and 3D)
  TTree *tree1 = new TTree("T1","RANDU");
  TTree *tree2 = new TTree("T2","TRandom");
  TTree *tree3 = new TTree("T3","TRadnom3");

  // Create Branches, easy to plot afterwards
  Float_t ax,ay,az;
  tree1->Branch("ax",&ax,"ax\F");
  tree1->Branch("ay",&ay,"ay\F");
  tree1->Branch("az",&az,"az\F");
  Float_t bx,by,bz;
  tree2->Branch("bx",&bx,"bx\F");
  tree2->Branch("by",&by,"by\F");
  tree2->Branch("bz",&bz,"bz\F");
  Float_t cx,cy,cz;
  tree3->Branch("cx",&cx,"cx\F");
  tree3->Branch("cy",&cy,"cy\F");
  tree3->Branch("cz",&cz,"cz\F");

  // Loop on the number of events
  for (int i=0;i<Events;i++){
    ax = ((Double_t ) Randu())/0x7fffffff;
    ay = ((Double_t ) Randu())/0x7fffffff;
    az = ((Double_t ) Randu())/0x7fffffff;
    tree1->Fill();
    bx = myRandom->Uniform();
    by = myRandom->Uniform();
    bz = myRandom->Uniform();
    tree2->Fill();
    cx = myRandom3->Uniform();
    cy = myRandom3->Uniform();
    cz = myRandom3->Uniform();
    tree3->Fill();
  
   }
  tree1->Scan();

  // Draw 1D, 2D and 3D
  c1->cd(1);
  tree1->Draw("ax");
  c1->cd(2);
  tree1->Draw("ax:ay");
  c1->cd(3);
  tree1->Draw("ax:ay:az");
  c1->cd(4);
  tree1->Draw("ax:ay:az");
  c2->cd(1);
  tree2->Draw("bx");
  c2->cd(2);
  tree2->Draw("bx:by");
  c2->cd(3);
  tree2->Draw("bx:by:bz");
  c3->cd(1);
  tree3->Draw("cx");
  c3->cd(2);
  tree3->Draw("cx:cy");
  c3->cd(3);
  tree3->Draw("cx:cy:cz");

  // we continue by validating the three generators using some of the 
  // basic methods introduced at the lecture. We don't expect any failure
  // for such simple tests... let see.

  // test on equal distribution
  TH1F* h = new TH1F("toed","toed",100,0,1);
  TH1F* toed_out = new TH1F("toed_result","toed_result",100,0,200);

  for(int pseudo = 0; pseudo<1000; ++pseudo) {
    // fill an histogram
    h->Reset();
    h->Clear();
    for(int i=0;i<10000; ++i) { h->Fill(((Double_t ) Randu())/0x7fffffff); }
    
    // then compute the test...
    float chi2test = 0;
    int Nok = h->GetEntries()/h->GetNbinsX();
    for(int i=1;i<=h->GetNbinsX();++i) {
      chi2test += pow(h->GetBinContent(i)-Nok,2)/Nok;
    }
    // ...and fill the summary histogram
    toed_out->Fill(chi2test);
  }
  // conclude the test by checking that we have a chi2 distribution
  TCanvas *c4 = new TCanvas("c4","test on equal distribution");
  toed_out->DrawNormalized();
  // compare to the chi2
  TH1F* chi2 = new TH1F("chi2","chi2",100,0,200);
  int k = 100;
  for(int i=0;i<200;++i) {
    Float_t x = chi2->GetBinCenter(i+1);
    chi2->SetBinContent(i+1,pow(x,k/2-1)*exp(-x/2)/(pow(2,k/2)*TMath::Gamma(k/2)));
  }
  chi2->DrawNormalized("same")->SetLineColor(4);
  std::cout << "TOED: The Kolmogorov test output for the comparison between the two is: " << chi2->KolmogorovTest(toed_out) << std::endl;
  
  // gap test: only the last on n in [a,b] -> P= p(1-p)^(n-1) with p=b-a
  // we will work with a=0.3, b=0.4, p=0.1, n=50 -> P = 0.00057264
  // we will therefore generate 100000 sequences of 50 numbers
  //TCanvas *c5 = new TCanvas("c5","gap test");
  int count = 0;
  for(int i=0;i<100000; ++i) {
    bool goodcase = false;
    for(int j=0;j<49;++j) {
      float x = (((Double_t ) Randu())/0x7fffffff);
      goodcase |= (x>0.3 && x<0.4);
    }
    float x = (((Double_t ) Randu())/0x7fffffff);
    if((!goodcase) && (x>0.3 && x<0.4)) count++;
  }
  std::cout << "Probability is " << count/100000. << ", to be compared to the theoretical 0.00057264168970223546" << std::endl;
  
  // random walk test
  // here we count the number of numbers among 1000 below 0.005 and compare that to a binomial
  TCanvas *c6 = new TCanvas("c6","random walk test");
  TH1F* rwt_out = new TH1F("rwt_result","rwt_result",25,0,25);
  for(int i=0;i<1000; ++i) {
    int count = 0;
    for(int j=0;j<1000;++j) {
      if((((Double_t ) Randu())/0x7fffffff)<0.005) count++;
    }
    rwt_out->Fill(count);
  }
  rwt_out->Draw();
  // compare to the binomial
  TH1F* binomial = new TH1F("binomial","binomial",25,0,25);
  for(int i=0;i<25;++i) {
    binomial->SetBinContent(i+1,1000*TMath::Binomial(1000, i)*TMath::Power(0.005, i)*TMath::Power(1-0.005, 1000-i));
  }
  binomial->Draw("same");
  binomial->SetLineColor(4);
  std::cout << "RWT: The Kolmogorov test output for the comparison between the two is: " << binomial->KolmogorovTest(rwt_out) << std::endl;

}