ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
pirndm.C
Go to the documentation of this file.
1 //Test program for random number generators (spped and quality)
2 //The program get n random number pairs x and y i [0,1]
3 //It counts the ratio of pairs in the circle of diameter 1
4 //compared to the total number of pairs.
5 //This ratio must be Pi/4
6 //The test shows the graph of the difference of this ratio compared to PI.
7 //To test only the speed, call pirndm(1) (default)
8 //To test quality and speed call pirndm(50)
9 // root pirndm.C+
10 //or
11 // root "pirndm.C+(10)"
12 //
13 //Author: Rene Brun
14 
15 #include "TROOT.h"
16 #include "TStopwatch.h"
17 #include "TMath.h"
18 #include "TRandom1.h"
19 #include "TRandom2.h"
20 #include "TRandom3.h"
21 #include "TCanvas.h"
22 #include "TH2.h"
23 #include "TGraph.h"
24 #include "TSystem.h"
25 #include "TLegend.h"
26 #include "TPaveLabel.h"
27 #include "TSystem.h"
28 #ifdef USE_MATHMORE
29 #include "Math/Random.h"
30 #include "Math/GSLRndmEngines.h"
31 #endif
32 #include <vector>
33 #include <iostream>
34 
35 TLegend *legend = 0;
36 TCanvas *c1 = 0;
39 
40 std::vector<TH2D *> vh2;
41 
42 
43 //____________________________________________________________________
44 template<class Random>
45 void piRandom(const char *name, Random *r, Long64_t n, Int_t color) {
46 
47  TH2D * h2 = new TH2D("h2",name,300,0,0.000003,300,0.,1.);
48 
49  timer.Start();
50  TGraph *gr = new TGraph();
51  gr->SetMarkerStyle(20);
52  gr->SetMarkerSize(0.9);
53  gr->SetMarkerColor(color);
54  gr->SetLineColor(color);
55  gr->SetLineWidth(2);
56 
57  Int_t k = 0;
58  Double_t diffpi;
59  Long64_t npi = 0;
60  Double_t pi = TMath::Pi();
61  const Int_t NR = 20000;
62  const Int_t NR2 = NR/2;
63  Double_t rn[NR];
64  Long64_t i = 0;
65  //double r1,r2;
66  while (i<=n) {
67  i += NR2;
68  r->RndmArray(NR,rn);
69  for (Int_t j=0;j<NR;j+=2) {
70  if (rn[j]*rn[j]+rn[j+1]*rn[j+1] <= 1) npi++;
71  if (rn[j] < 0.001) h2->Fill(rn[j],rn[j+1]);
72  }
73 // r1 = r->Rndm();
74 // r2 = r->Rndm();
75  // if (r1*r1+r2*r2 <= 1) npi++;
76  if (i && i % (n/10) == 0) {
78  Double_t norm = 4./Double_t(i);
79  diffpi = norm*npi - pi;
80  gr->SetPoint(k,i,diffpi);
81  if (k ==0) gr->Draw("lp");
82  else {
83  c1->Modified();
84  c1->Update();
85  }
86  k++;
87  }
88  }
89  timer.Stop();
90  Double_t cpu = timer.CpuTime();
91  cputot += cpu;
92  Double_t nanos = 1.e9*cpu/Double_t(2*n);
93  legend->AddEntry(gr,Form("%-14s: %6.1f ns/call",name,nanos),"lp");
94  c1->Modified();
95  c1->Update();
96  printf("RANDOM = %s : RT=%7.3f s, Cpu=%7.3f s\n",name,timer.RealTime(),cpu);
97 
98 // TCanvas * c2 = new TCanvas();
99 // h2->Draw();
100 // c2->Update();
101 
102 // c1->SetSelected(c1);
103  vh2.push_back(h2);
104 
105 }
106 
107 //________________________________________________________________________
109  Int_t np = 40;
110  TGraph *g = new TGraph(2*np+2);
111  Double_t xmax = Double_t(n)/Double_t(np);
112  for (Int_t i=1;i<=np;i++) {
113  Double_t x = i*xmax;
114  //Double_t e = 1./TMath::Sqrt(x);
115  Double_t e = TMath::Sqrt( 2 * TMath::Pi() * (4 - TMath::Pi() )/x );
116  g->SetPoint(i,x,e);
117  g->SetPoint(2*np-i+1,x,-e);
118  }
119  Double_t x0 = 0.1*xmax;
120  Double_t e0 = 1./TMath::Sqrt(x0);
121  g->SetPoint(0,x0,e0);
122  g->SetPoint(2*np+1,0,-e0);
123  g->SetPoint(2*np+2,0,e0);
124  g->SetFillColor(1);
125  g->SetFillStyle(3002);
126  g->Draw("f");
127 }
128 
129 
130 //________________________________________________________________________
131 void pirndm(Long64_t n1=1, unsigned int seed = 0) {
132  Long64_t n = n1*20000000;
133  c1 = new TCanvas("c1");
134  c1->SetLeftMargin(0.12);
135  c1->SetFrameFillColor(41);
136  c1->SetFrameBorderSize(6);
137  c1->SetGrid();
138  Double_t dy = 10/TMath::Sqrt(n);
139  // Double_t dy = 1.5e-3;
140  //if (n1 < 4) dy *=2;
141  TH2F *frame = new TH2F("h","",100,0,1.1*n,100,-dy,dy);
142  frame->GetXaxis()->SetTitle("Number of TRandom calls");
143  frame->GetYaxis()->SetTitle("Difference with #pi");
144  frame->GetYaxis()->SetTitleOffset(1.3);
145  frame->GetYaxis()->SetDecimals();
146  frame->SetStats(0);
147  frame->Draw();
148  legend = new TLegend(0.6,0.7,0.88,0.88);
149  legend->Draw();
150 
151 
152  ErrorBand(n);
153  std::cout << "seed is " << seed << std::endl;
154 
155  piRandom("TRandom",new TRandom(seed),n,kYellow);
156  piRandom("TRandom2",new TRandom2(seed),n,kBlue);
157  piRandom("TRandom3",new TRandom3(seed),n,kRed);
158  piRandom("TRandom1",new TRandom1(seed),n,kGreen);
159 
160 #ifdef USE_MATHMORE
161 #define GSL2
162 #ifdef GSL1
163  piRandom("TRandom()",new TRandom(),n,kYellow);
164  piRandom("TRandom3(0)",new TRandom3(0),n,kBlack);
169 #endif
170 #ifdef GSL2
171  piRandom("TRandom(1)",new TRandom(1),n,kYellow);
172  piRandom("TRandom(2^30)",new TRandom(1073741824),n,kRed);
173  piRandom("TRandom(256)",new TRandom(256),n,kMagenta);
174  piRandom("TRandom(2)",new TRandom(2),n,kBlue);
176  // piRandom("RANMAR",new ROOT::Math::Random<ROOT::Math::GSLRngRanMar>(),n,kBlue);
177  //piRandom("MINSTD",new ROOT::Math::Random<ROOT::Math::GSLRngMinStd>(),n,kBlack);
178  piRandom("new TRandom2",new TRandom2(),n,kCyan);
179  //piRandom("TRandom3",new TRandom3(),n,kCyan);
180 // piRandom("CMRG",new ROOT::Math::Random<ROOT::Math::GSLRngCMRG>(),n,kRed);
181 // piRandom("MRG",new ROOT::Math::Random<ROOT::Math::GSLRngMRG>(),n,kMagenta);
182 #else
183  piRandom("new TRandom2",new TRandom2(std::rand()),n,kYellow);
184  piRandom("new TRandom2",new TRandom2(std::rand()),n,kRed);
185  piRandom("new TRandom2",new TRandom2(std::rand()),n,kMagenta);
186  piRandom("new TRandom2",new TRandom2(std::rand()),n,kBlack);
187  piRandom("new TRandom2",new TRandom2(std::rand()),n,kBlue);
188  piRandom("new TRandom2",new TRandom2(std::rand()),n,kRed);
189  piRandom("new TRandom2",new TRandom2(std::rand()),n,kGreen);
190 #endif
191 #endif
192 
193  // reftime calculated on MACOS Intel dualcore 2GHz
194  // time for TRandom + TRandom2 + TRandom3 + TRandom1 for n = 10**7 (n1=5000)
195  Double_t reftime = (4629.530 + 5358.100 + 5785.240 + 26012.17)/5000.;
196  const Double_t rootmarks = 900*Double_t(n1)*reftime/cputot;
197  TPaveLabel *pl = new TPaveLabel(0.2,0.92,0.8,0.98,Form("cpu time = %6.1fs - rootmarks = %6.1f",cputot,rootmarks),"brNDC");
198  pl->Draw();
199  printf("******************************************************************\n");
200  printf("* ROOTMARKS =%6.1f * Root%-8s %d/%d\n",rootmarks,gROOT->GetVersion(),
201  gROOT->GetVersionDate(),gROOT->GetVersionTime());
202  printf("******************************************************************\n");
203 
204  printf("Time at the end of job = %f seconds\n",cputot);
205  c1->Print("pirndm.root");
206  c1->Print("pirndm.gif");
207 
208 
209  // draw 2D histos
210  TCanvas * c2 = new TCanvas();
211  int nx = 0;
212  int ny = vh2.size();
213  for ( nx = 1; nx < ny; nx++) {
214  double r = double(vh2.size())/nx;
215  ny = int(r - 0.01) + 1;
216  }
217  nx--;
218  std::cout << nx << " " << ny << " " << vh2.size() << std::endl;
219 
220  c2->Divide(ny,nx);
221  for (unsigned int i = 0; i < vh2.size(); ++i) {
222  c2->cd (i+1);
223  vh2[i]->Draw();
224  }
225  c2->Update();
226 
227 }
const int nx
Definition: kalman.C:16
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title Offset is a correction factor with respect to the "s...
Definition: TAttAxis.cxx:245
virtual void SetLineWidth(Width_t lwidth)
Definition: TAttLine.h:57
virtual Bool_t ProcessEvents()
Process pending events (GUI, timers, sockets).
Definition: TSystem.cxx:420
Random number generator class based on M.
Definition: TRandom3.h:29
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
long long Long64_t
Definition: RtypesCore.h:69
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:35
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
const double pi
Definition: Rtypes.h:61
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:29
Double_t cputot
Definition: pirndm.C:38
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
Definition: Rtypes.h:60
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
#define gROOT
Definition: TROOT.h:344
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
int Int_t
Definition: RtypesCore.h:41
virtual void SetFillStyle(Style_t fstyle)
Definition: TAttFill.h:52
Definition: Rtypes.h:61
Definition: Rtypes.h:61
virtual void Print(const char *filename="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:4134
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:740
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
The Ranlux Random number generator class.
Definition: TRandom1.h:29
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:29
const int ny
Definition: kalman.C:17
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)
Definition: TPad.h:326
virtual void SetMarkerColor(Color_t mcolor=1)
Definition: TAttMarker.h:51
tuple np
Definition: multifit.py:30
TH2D * h2
Definition: fit2dHist.C:45
A Pave (see TPave) with a text centered in the Pave.
Definition: TPaveLabel.h:32
virtual void SetLineColor(Color_t lcolor)
Definition: TAttLine.h:54
void ErrorBand(Long64_t n)
Definition: pirndm.C:108
ROOT::R::TRInterface & r
Definition: Object.C:4
R__EXTERN TSystem * gSystem
Definition: TSystem.h:545
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2878
virtual void SetFillColor(Color_t fcolor)
Definition: TAttFill.h:50
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:256
char * Form(const char *fmt,...)
std::vector< TH2D * > vh2
Definition: pirndm.C:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Definition: TAttMarker.h:53
tuple pl
Definition: first.py:10
TAxis * GetYaxis()
Definition: TH1.h:320
float xmax
Definition: THbookFile.cxx:93
virtual void SetMarkerSize(Size_t msize=1)
Definition: TAttMarker.h:54
TGraphErrors * gr
Definition: legend1.C:25
void piRandom(const char *name, Random *r, Long64_t n, Int_t color)
Definition: pirndm.C:45
Double_t Pi()
Definition: TMath.h:44
The Canvas class.
Definition: TCanvas.h:48
return c2
Definition: legend2.C:14
virtual void Draw(Option_t *option="")
Draw this pavelabel with its current attributes.
Definition: TPaveLabel.cxx:77
void Random()
Definition: utils.cpp:154
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:280
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
void rootmarks()
Definition: rootmarks.C:1
Definition: Rtypes.h:61
void SetDecimals(Bool_t dot=kTRUE)
Sets the decimals flag By default, blank characters are stripped, and then the label is correctly ali...
Definition: TAxis.h:203
#define name(a, b)
Definition: linkTestLib0.cpp:5
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2127
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:1073
void SetFrameBorderSize(Width_t size=1)
Definition: TAttPad.h:88
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
Definition: Rtypes.h:61
void SetFrameFillColor(Color_t color=1)
Definition: TAttPad.h:83
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2179
void pirndm(Long64_t n1=1, unsigned int seed=0)
Definition: pirndm.C:131
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:287
virtual void SetTitle(const char *title="")
Change (i.e. set) the title of the TNamed.
Definition: TNamed.cxx:152
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
const Int_t n
Definition: legend1.C:16
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8320
void Modified(Bool_t flag=1)
Definition: TPad.h:407
TAxis * GetXaxis()
Definition: TH1.h:319
Documentation for the Random class.
Definition: Random.h:41
virtual void SetLeftMargin(Float_t leftmargin)
Set Pad left margin in fraction of the pad width.
Definition: TAttPad.cxx:107
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:297
Stopwatch class.
Definition: TStopwatch.h:30