Logo ROOT   6.12/07
Reference Guide
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 
20 using namespace ROOT::Math;
21 
22 int 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 }
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
return c1
Definition: legend1.C:41
User class for MathMore random numbers template on the Engine type.
Definition: QuasiRandom.h:60
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
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:67
Double_t x[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
The Canvas class.
Definition: TCanvas.h:31
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:290
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:84
#define gPad
Definition: TVirtualPad.h:285
virtual Int_t GetNbinsX() const
Definition: TH1.h:291
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
static char * skip(char **buf, const char *delimiters)
Definition: civetweb.c:2039
const Int_t n
Definition: legend1.C:16
virtual Int_t GetNbinsY() const
Definition: TH1.h:292
Documentation for the Random class.
Definition: Random.h:39
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290
Stopwatch class.
Definition: TStopwatch.h:28