ROOT   6.10/09 Reference Guide
piRandom.C
Go to the documentation of this file.
1 #include "Math/Random.h"
2 #include "Math/GSLRndmEngines.h"
3 #include "TStopwatch.h"
4 #include "TRandom1.h"
5 #include "TRandom2.h"
6 #include "TRandom3.h"
7 #include <iostream>
8 #include <cmath>
9 #include <typeinfo>
10 #include "TH1D.h"
11 #include "TStyle.h"
12 #include "TF1.h"
13 #include "TPaveLabel.h"
14
15 #include "TCanvas.h"
16
17 #ifndef PI
18 #define PI 3.14159265358979323846264338328 /* pi */
19 #endif
20
21 #define NLOOP 1000;
22 //#define NEVT 20000000;
23 #define NEVT 10000;
24
25 using namespace ROOT::Math;
26
27 template <class R>
28 void generate( R & r, TH1D * h) {
29
30  TStopwatch w;
31
32  r.SetSeed(0);
33  //r.SetSeed(int(std::pow(2.0,28)));
34  int m = NLOOP;
35  int n = NEVT;
36  for (int j = 0; j < m; ++j) {
37
38  //std::cout << r.GetSeed() << " ";
39
40  w.Start();
41 // if ( n < 40000000) iseed = std::rand();
42 // iseed = 0;
43  //TRandom3 r3(0);
44  //r.SetSeed( 0 ); // generate random seeds
45  //TRandom3 r3(0);
46  //r.SetSeed (static_cast<UInt_t> (4294967296.*r3.Rndm()) );
47
48  // estimate PI
49  double n1=0;
50  double rn[2000];
51  double x;
52  double y;
53  for (int ievt = 0; ievt < n; ievt+=1000 ) {
54  r.RndmArray(2000,rn);
55  for (int i=0; i < 1000; i++) {
56  x=rn[2*i];
57  y=rn[2*i+1];
58  if ( ( x*x + y*y ) <= 1.0 ) n1++;
59  }
60  }
61  double piEstimate = 4.0 * double(n1)/double(n);
62  double delta = piEstimate-PI;
63  h->Fill(delta);
64  }
65
66  w.Stop();
67  std::cout << std::endl;
68  std::cout << "Random: " << typeid(r).name()
69  << "\n\tTime = " << w.RealTime() << " " << w.CpuTime() << std::endl;
70  std::cout << "Time/call: " << w.CpuTime()/(2*n)*1.0E9 << std::endl;
71 }
72
73 int piRandom() {
74
75  TRandom r0;
76  TRandom1 r1;
77  TRandom2 r2;
78  TRandom3 r3;
79  //Random<GSLRngRand> r1;
80 // Random<GSLRngTaus> r2;
81 // Random<GSLRngRanLux> r3;
82
83  double n = NEVT;
84  int nloop = NLOOP;
85  double dy = 15/std::sqrt(n);
86
87  TH1D * h0 = new TH1D("h0","TRandom delta",100,-dy,dy);
88  TH1D * h1 = new TH1D("h1","TRandom1 delta",100,-dy,dy);
89  TH1D * h2 = new TH1D("h2","TRandom2 delta",100,-dy,dy);
90  TH1D * h3 = new TH1D("h3","TRandom3 delta",100,-dy,dy);
91
92  double sigma = std::sqrt( PI * (4 - PI)/n );
93  std::cout << "**********************************************************" << std::endl;
94  std::cout << " Generate " << int(n) << " for " << nloop << " times " << std::endl;
95  std::cout << "**********************************************************" << std::endl;
96  std::cout << "\tExpected Sigma = " << sigma << std::endl;
97
98 #define INTERACTIVE
99 #ifdef INTERACTIVE
100
101  double del, err;
102  TCanvas * c1 = new TCanvas("c1_piRandom","PI Residuals");
103  gStyle->SetOptFit(1111);
104  gStyle->SetOptLogy();
105  c1->Divide(2,2);
106  c1->cd(1);
107  generate(r0,h0);
108  h0->Fit("gaus");
109  h0->Draw();
110  TF1 * fg = (TF1*) h0->FindObject("gaus");
111  if (fg) {
112  del = (fg->GetParameter(2)-sigma);
113  err = fg->GetParError(2);
114  }
115  else { del = -999; err = 1; }
116
117  char text[20];
119  TPaveLabel * pl0 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
120  pl0->Draw();
121
122
123
124
125  c1->cd(2);
126  generate(r1,h1);
127  h1->Fit("gaus");
128  h1->Draw();
129  fg = (TF1*) h1->FindObject("gaus");
130  if (fg) {
131  del = (fg->GetParameter(2)-sigma);
132  err = fg->GetParError(2);
133  }
134  else { del = -999; err = 1; }
135
137  TPaveLabel * pl1 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
138  pl1->Draw();
139
140
141  c1->cd(3);
142  generate(r2,h2);
143  h2->Fit("gaus");
144  h2->Draw();
145  fg = (TF1*) h2->FindObject("gaus");
146  if (fg) {
147  del = (fg->GetParameter(2)-sigma);
148  err = fg->GetParError(2);
149  }
150  else { del = -999; err = 1; }
151
153  TPaveLabel * pl2 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
154  pl2->Draw();
155
156
157  c1->cd(4);
158  generate(r3,h3);
159  h3->Fit("gaus");
160  h3->Draw();
161  fg = (TF1*) h3->FindObject("gaus");
162  if (fg) {
163  del = (fg->GetParameter(2)-sigma);
164  err = fg->GetParError(2);
165  }
166  else { del = -999; err = 1; }
167
169  TPaveLabel * pl3 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
170  pl3->Draw();
171
172 #else
173  generate(r0,h0);
174  generate(r1,h1);
175  generate(r2,h2);
176  generate(r3,h3);
177 #endif
178
179  return 0;
180
181 }
182
183 int main() {
184
185  piRandom();
186  return 0;
187 }
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3126
Random number generator class based on M.
Definition: TRandom3.h:27
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
return c1
Definition: legend1.C:41
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:27
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
TH1 * h
Definition: legend2.C:5
Definition: TCanvas.cxx:679
#define PI
Definition: piRandom.C:18
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
void generate(R &r, TH1D *h)
Definition: piRandom.C:28
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1656
Float_t delta
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
The Ranlux Random number generator class.
Definition: TRandom1.h:27
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
unsigned int r3[N_CITIES]
Definition: simanTSP.cxx:323
#define NLOOP
Definition: piRandom.C:21
int piRandom()
Definition: piRandom.C:73
const Double_t sigma
TH1F * h1
Definition: legend1.C:5
A Pave (see TPave) with a text centered in the Pave.
Definition: TPaveLabel.h:20
TRandom2 r(17)
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1219
TMarker * m
Definition: textangle.C:8
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:594
virtual TObject * FindObject(const char *name) const
Search object named name in the list of functions.
Definition: TH1.cxx:3537
The Canvas class.
Definition: TCanvas.h:31
virtual void Draw(Option_t *option="")
Draw this pavelabel with its current attributes.
Definition: TPaveLabel.cxx:77
TText * text
Double_t y[n]
Definition: legend1.C:17
int main()
Definition: piRandom.C:183
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)
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:466
1-Dim function class
Definition: TF1.h:150
void SetOptLogy(Int_t logy=1)
Definition: TStyle.h:305
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3564
#define NEVT
Definition: piRandom.C:23
const Int_t n
Definition: legend1.C:16
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
Stopwatch class.
Definition: TStopwatch.h:28