Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
quasirandom.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook -js
4/// Example of generating quasi-random numbers
5///
6/// \macro_image
7/// \macro_output
8/// \macro_code
9///
10/// \author Lorenzo Moneta
11
12#include "Math/QuasiRandom.h"
13#include "Math/Random.h"
14#include "TH2.h"
15#include "TCanvas.h"
16#include "TStopwatch.h"
17
18#include <iostream>
19
20using namespace ROOT::Math;
21
22int quasirandom(int n = 10000, int skip = 0) {
23
24 TH2D * h0 = new TH2D("h0","Pseudo-random Sequence",200,0,1,200,0,1);
25 TH2D * h1 = new TH2D("h1","Sobol Sequence",200,0,1,200,0,1);
26 TH2D * h2 = new TH2D("h2","Niederrer Sequence",200,0,1,200,0,1);
27
28 RandomMT r0;
29 // quasi random numbers need to be created giving the dimension of the sequence
30 // in this case we generate 2-d sequence
31
32 QuasiRandomSobol r1(2);
34
35 // generate n random points
36
37 double x[2];
38 TStopwatch w; w.Start();
39 for (int i = 0; i < n; ++i) {
40 r0.RndmArray(2,x);
41 h0->Fill(x[0],x[1]);
42 }
43 std::cout << "Time for gRandom ";
44 w.Print();
45
46 w.Start();
47 if( skip>0) r1.Skip(skip);
48 for (int i = 0; i < n; ++i) {
49 r1.Next(x);
50 h1->Fill(x[0],x[1]);
51 }
52 std::cout << "Time for Sobol ";
53 w.Print();
54
55 w.Start();
56 if( skip>0) r2.Skip(skip);
57 for (int i = 0; i < n; ++i) {
58 r2.Next(x);
59 h2->Fill(x[0],x[1]);
60 }
61 std::cout << "Time for Niederreiter ";
62 w.Print();
63
64 TCanvas * c1 = new TCanvas("c1","Random sequence",600,1200);
65 c1->Divide(1,3);
66 c1->cd(1);
67 h0->Draw("COLZ");
68 c1->cd(2);
69
70 // check uniformity
71 h1->Draw("COLZ");
72 c1->cd(3);
73 h2->Draw("COLZ");
74 gPad->Update();
75
76 // test number of empty bins
77
78 int nzerobins0 = 0;
79 int nzerobins1 = 0;
80 int nzerobins2 = 0;
81 for (int i = 1; i <= h1->GetNbinsX(); ++i) {
82 for (int j = 1; j <= h1->GetNbinsY(); ++j) {
83 if (h0->GetBinContent(i,j) == 0 ) nzerobins0++;
84 if (h1->GetBinContent(i,j) == 0 ) nzerobins1++;
85 if (h2->GetBinContent(i,j) == 0 ) nzerobins2++;
86 }
87 }
88
89 std::cout << "number of empty bins for pseudo-random = " << nzerobins0 << std::endl;
90 std::cout << "number of empty bins for " << r1.Name() << "\t= " << nzerobins1 << std::endl;
91 std::cout << "number of empty bins for " << r2.Name() << "\t= " << nzerobins2 << std::endl;
92
93 int iret = 0;
94 if (nzerobins1 >= nzerobins0 ) iret += 1;
95 if (nzerobins2 >= nzerobins0 ) iret += 2;
96 return iret;
97
98}
#define gPad
User class for MathMore random numbers template on the Engine type.
Definition QuasiRandom.h:60
Documentation for the Random class.
Definition Random.h:42
void RndmArray(int n, double *array)
Generate an array of random numbers between ]0,1] 0 is excluded and 1 is included Function to preserv...
Definition Random.h:70
The Canvas class.
Definition TCanvas.h:23
virtual Int_t GetNbinsY() const
Definition TH1.h:297
virtual Int_t GetNbinsX() const
Definition TH1.h:296
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4993
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:292
Int_t Fill(Double_t)
Invalid Fill method.
Definition TH2.cxx:294
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH2.h:88
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5