Logo ROOT   6.10/09
Reference Guide
testRandom.cxx
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 <cstdlib>
10 #include <typeinfo>
11 #include <iomanip>
12 #include <fstream>
13 
14 #include <limits>
15 
16 #ifdef HAVE_CLHEP
17 #include "CLHEP/Random/RandFlat.h"
18 #include "CLHEP/Random/RanluxEngine.h"
19 #include "CLHEP/Random/Ranlux64Engine.h"
20 #endif
21 
22 
23 
24 #ifndef PI
25 #define PI 3.14159265358979323846264338328 /* pi */
26 #endif
27 
28 //#define NEVT 10000
29 #ifndef NEVT
30 #define NEVT 10000000
31 #endif
32 
33 using namespace ROOT::Math;
34 
35 #ifdef HAVE_CLHEP
36 using namespace CLHEP;
37 #endif
38 
39 
40 // wrapper for stdrand
41 
42 class RandomStd {
43 public:
44  RandomStd() {
45  fScale = 1./(double(RAND_MAX) + 1.);
46  }
47 
48  inline void RndmArray(int n, double * x) {
49  for ( double * itr = x; itr != x+n; ++itr)
50  *itr= fScale*double(std::rand());
51  }
52  inline double Uniform() {
53  return fScale*double(std::rand());
54  }
55 
56  std::string Type() const { return "std::rand"; }
57  unsigned int EngineSize() const { return 0; }
58 
59  void SetSeed(int seed) { std::srand(seed); }
60 
61 private:
62  double fScale;
63 
64 };
65 
66 
67 
68 // wrapper for clhep
69 
70 template <class Engine>
71 class RandomCLHEP {
72 public:
73  RandomCLHEP(Engine & e) :
74  fRand(e)
75  {
76  }
77 
78  inline void RndmArray(int n, double * x) {
79  fRand.flatArray(n,x);
80  }
81  inline double Uniform() {
82  return fRand.flat();
83  }
84 
85  std::string Type() const { return std::string("CLHEP ") + Engine::engineName(); }
86  unsigned int EngineSize() const { return 0; }
87 
88  void SetSeed(int seed) { fRand.setSeed(seed); }
89 
90 private:
91  Engine & fRand;
92 
93 };
94 
95 
96 
97 
98 
99 template <class R>
100 void printName( const R & r) {
101  std::cout << "\nRandom :\t " << r.Type() << " \t size of state = " << r.EngineSize() << std::endl;
102 }
103 
104 // specializations for TRandom's
105 void printName( const TRandom & r) {
106  std::cout << "\nRandom :\t " << r.ClassName() << std::endl;
107 }
108 // specializations for TRandom's
109 void printName( const TRandom1 & r) {
110  std::cout << "\nRandom :\t " << r.ClassName() << std::endl;
111 }
112 // specializations for TRandom's
113 void printName( const TRandom2 & r) {
114  std::cout << "\nRandom :\t " << r.ClassName() << std::endl;
115 }
116 // specializations for TRandom's
117 void printName( const TRandom3 & r) {
118  std::cout << "\nRandom :\t " << r.ClassName() << std::endl;
119 }
120 
121 template <class R>
122 void generate( R & r, bool array=true) {
123 
124  TStopwatch w;
125 
126  int n = NEVT;
127  // estimate PI
128  double n1=0;
129  double x,y;
130  w.Start();
131  // use default seeds
132  // r.SetSeed(0);
133  //r.SetSeed(int(pow(2,28)) );
134  if (array) {
135  double ax[1000];
136  double ay[1000];
137  for (int i = 0; i < n; i+=1000 ) {
138  r.RndmArray(1000,ax);
139  r.RndmArray(1000,ay);
140  for (int j = 0; j < 1000; ++j)
141  if ( ( ax[j]*ax[j] + ay[j]*ay[j] ) <= 1.0 ) n1++;
142  }
143  }
144  else {
145  for (int i = 0; i < n; ++i) {
146  x=r.Uniform();
147  y=r.Uniform();
148  if ( ( x*x + y*y ) <= 1.0 ) n1++;
149  }
150  }
151  w.Stop();
152 
153  printName(r);
154  std::cout << "\tTime = " << w.RealTime()*1.0E9/NEVT << " "
155  << w.CpuTime()*1.0E9/NEVT
156  << " (ns/call)" << std::endl;
157  double piEstimate = 4.0 * double(n1)/double(n);
158  double delta = piEstimate-PI;
159  double sigma = std::sqrt( PI * (4 - PI)/double(n) );
160  std::cout << "\t\tDeltaPI = " << delta/sigma << " (sigma) " << std::endl;
161 }
162 
163 int main() {
164 
165  std::cout << "***************************************************\n";
166  std::cout << " TEST RANDOM NEVT = " << NEVT << std::endl;
167  std::cout << "***************************************************\n\n";
168 
169 
170 
184  RandomStd r10;
185 
186  TRandom tr0;
187  TRandom1 tr1;
188  TRandom1 tr1a(0,0);
189  TRandom1 tr1b(0,1);
190  TRandom1 tr1c(0,2);
191  TRandom1 tr1d(0,3);
192  TRandom1 tr1e(0,4);
193  TRandom2 tr2;
194  TRandom3 tr3;
195 
196 
197  generate(tr0);
198  generate(tr1);
199  generate(tr1a);
200  generate(tr1b);
201  generate(tr1c);
202  generate(tr1d);
203  generate(tr1e);
204  generate(tr2);
205  generate(tr3);
206 
207  generate(r10);
208 
209  generate(r1);
210  generate(r2);
211  generate(r3);
212  generate(r3s1);
213  generate(r3s2);
214  generate(r3d1);
215  generate(r3d2);
216  generate(r4);
217  generate(r5);
218  generate(r6);
219  generate(r7);
220  generate(r8);
221  generate(r9);
222 
223 
224 #ifdef HAVE_CLHEP
225  RanluxEngine e1(1,3);
226  RanluxEngine e2(1,4);
227  Ranlux64Engine e3(1,0);
228  Ranlux64Engine e4(1,1);
229  Ranlux64Engine e5(1,2);
230 
231  RandomCLHEP<RanluxEngine> crlx3(e1);
232  RandomCLHEP<RanluxEngine> crlx4(e2);
233 
234  RandomCLHEP<Ranlux64Engine> crlx64a(e3);
235  RandomCLHEP<Ranlux64Engine> crlx64b(e4);
236  RandomCLHEP<Ranlux64Engine> crlx64c(e5);
237 
238  generate(crlx3);
239  generate(crlx4);
240 
241  generate(crlx64a);
242  generate(crlx64b);
243  generate(crlx64c);
244 
245 #endif
246 
247  // generate 1000 number with GSL MT and check with TRandom3
248  int n = 1000;
249  std::vector<double> v1(n);
250  std::vector<double> v2(n);
251 
252  Random<GSLRngMT> gslRndm(4357);
253  TRandom3 rootRndm(4357);
254 
255 
256  gslRndm.RndmArray(n,&v1[0]);
257  rootRndm.RndmArray(n,&v2[0]);
258 
259  int nfail=0;
260  for (int i = 0; i < n; ++i) {
261  double d = std::fabs(v1[i] - v2[i] );
262  if (d > std::numeric_limits<double>::epsilon()*v1[i] ) nfail++;
263  }
264  if (nfail > 0) {
265  std::cout << "ERROR: Test failing comparing TRandom3 with GSL MT" << std::endl;
266  return -1;
267  }
268  // save the generated number
269  std::ofstream file("testRandom.out");
270  std::ostream & out = file;
271  int j = 0;
272  int prec = std::cout.precision(9);
273  while ( j < n) {
274  for (int l = 0; l < 8; ++l) {
275  out << std::setw(12) << v1[j+l] << ",";
276 // int nws = int(-log10(v1[j+l]));
277 // for (int k = nws; k >= 0; --k)
278 // out << " ";
279  }
280  out << std::endl;
281  j+= 8;
282  }
283  std::cout.precision(prec);
284 
285  return 0;
286 
287 }
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 generate(R &r, bool array=true)
Definition: testRandom.cxx:122
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:27
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 CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
Float_t delta
double sqrt(double)
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:135
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
const Double_t sigma
#define NEVT
Definition: testRandom.cxx:30
void printName(const R &r)
Definition: testRandom.cxx:100
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TRandom2 r(17)
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
TLine * l
Definition: textangle.C:4
int main()
Definition: testRandom.cxx:163
REAL epsilon
Definition: triangle.c:617
Type
enumeration specifying the integration types.
virtual void RndmArray(Int_t n, Float_t *array)
Return an array of n random numbers uniformly distributed in ]0,1].
Definition: TRandom3.cxx:138
#define PI
Definition: testRandom.cxx:25
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Definition: file.py:1
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
Documentation for the Random class.
Definition: Random.h:39
Stopwatch class.
Definition: TStopwatch.h:28