Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TRandom.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Author: Rene Brun, Lorenzo Moneta 15/12/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /**
13 
14 \class TRandom
15 
16 @ingroup Random
17 
18 This is the base class for the ROOT Random number generators.
19 This class defines the ROOT Random number interface and it should not be instantiated directly but used via its derived classes.
20 The generator provided in TRandom itself is a LCG (Linear Congruential Generator), the <a href="https://www.gnu.org/software/gsl/manual/html_node/Unix-random-number-generators.html">BSD `rand`
21 generator</a>, that it should not be used because its period is only 2**31, i.e. approximatly 2 billion events, that can be generated in just few seconds.
22 
23 To generate random numbers, one should use the derived class, which are :
24 - TRandom3: it is based on the "Mersenne Twister generator",
25 it is fast and a very long period of about \f$10^{6000}\f$. However it fails some of the most stringent tests of the
26 <a href="http://simul.iro.umontreal.ca/testu01/tu01.html">TestU01 suite</a>.
27 In addition this generator provide only numbers with 32 random bits, which might be not sufficient for some application based on double or extended precision.
28 This generator is however used in ROOT used to instantiate the global pointer to the ROOT generator, *gRandom*.
29 - ::TRandomMixMax: Generator based on the family of the MIXMAX matrix generators (see the <a href="https://mixmax.hepforge.org">MIXMAX HEPFORGE Web page</a> and the
30  the documentation of the class ROOT::Math::MixMaxEngine for more information), that are base on the Asanov dynamical C systems.
31 This generator has a state of N=240 64 bit integers, proof random properties, it provides 61 random bits and it has a very large period (\f$10^{4839}\f$).
32 Furthermore, it provides the capability to be seeded with the guarantee that, for each given different seed, a different sequence of random numbers will be generated.
33 The only drawback is that the seeding time is time consuming, of the order of 0.1 ms, while the time to generate a number is few ns (more than 10000 faster).
34 - ::TRandomMixMax17: Another MixMax generator, but with a smaller state, N=17, and this results in a smaller entropy than the generator with N=240. However, it has the same seeding capabilities, with a much faster seeding time (about 200 times less than TRandomMixMax240 and comparable to TRandom3).
35 - ::TRandomMixMax256 : A variant of the MIXMAX generators, based on a state of N=256, and described in the
36  <a href="http://arxiv.org/abs/1403.5355">2015 paper</a>. This implementation has been modified with respect to the paper, by skipping 2 internal interations,
37  to provide improved random properties.
38 - TRandomMT64 : Generator based on a the Mersenne-Twister generator with 64 bits,
39  using the implementation provided by the standard library ( <a href="http://www.cplusplus.com/reference/random/mt19937_64/">std::mt19937_64</a> )
40 - TRandom1 based on the RANLUX algorithm, has mathematically proven random proprieties
41  and a period of about \f$10{171}\f$. It is however much slower than the others and it has only 24 random bits. It can be constructed with different luxury levels.
42 - TRandomRanlux48 : Generator based on a the RanLux generator with 48 bits and highest luxury level
43  using the implementation provided by the standard library (<a href="http://www.cplusplus.com/reference/random/ranlux48/">std::ranlux48</a>). The drawback of this generator is its slow generation
44  time.
45 - TRandom2 is based on the Tausworthe generator of L'Ecuyer, and it has the advantage
46 of being fast and using only 3 words (of 32 bits) for the state. The period however is not impressively long, it is 10**26.
47 
48 Using the template TRandomGen class (template on the contained Engine type), it is possible to add any generator based on the standard C++ random library
49 (see the C++ <a href="http://www.cplusplus.com/reference/random/">random</a> documentation.) or different variants of the MIXMAX generator using the
50 ROOT::Math::MixMaxEngine. Some of the listed generator above (e.g. TRandomMixMax256 or TRandomMT64) are convenient typedef's of generator built using the
51 template TRandomGen class.
52 
53 Please note also that this class (TRandom) implements also a very simple generator (linear congruential) with period = \f$10^9\f$, known to have defects (the lower random bits are correlated) and it
54 is failing the majority of the random number generator tests. Therefore it should NOT be used in any statistical study.
55 
56 The following table shows some timings (in nanoseconds/call)
57 for the random numbers obtained using a macbookpro 2.6 GHz Intel Core i7 CPU:
58 
59 
60 - TRandom 3 ns/call (but this is a very BAD Generator, not to be used)
61 - TRandom2 5 ns/call
62 - TRandom3 5 ns/call
63 - ::TRandomMixMax 6 ns/call
64 - ::TRandomMixMax17 6 ns/call
65 - ::TRandomMT64 9 ns/call
66 - ::TRandomMixMax256 10 ns/call
67 - ::TRandom1 80 ns/call
68 - ::TRandomRanlux48 250 ns/call
69 
70 The following methods are provided to generate random numbers disctributed according to some basic distributions:
71 
72 - `::Exp(tau)`
73 - `::Integer(imax)`
74 - `::Gaus(mean,sigma)`
75 - `::Rndm()`
76 - `::Uniform(x1)`
77 - `::Landau(mpv,sigma)`
78 - `::Poisson(mean)`
79 - `::Binomial(ntot,prob)`
80 
81 Random numbers distributed according to 1-d, 2-d or 3-d distributions contained in TF1, TF2 or TF3 objects can also be generated.
82 For example, to get a random number distributed following abs(sin(x)/x)*sqrt(x)
83 you can do :
84 \code{.cpp}
85  TF1 *f1 = new TF1("f1","abs(sin(x)/x)*sqrt(x)",0,10);
86  double r = f1->GetRandom();
87 \endcode
88 or you can use the UNURAN package. You need in this case to initialize UNURAN
89 to the function you would like to generate.
90 \code{.cpp}
91  TUnuran u;
92  u.Init(TUnuranDistrCont(f1));
93  double r = u.Sample();
94 \endcode
95 
96 The techniques of using directly a TF1,2 or 3 function is powerful and
97 can be used to generate numbers in the defined range of the function.
98 Getting a number from a TF1,2,3 function is also quite fast.
99 UNURAN is a powerful and flexible tool which containes various methods for
100 generate random numbers for continuous distributions of one and multi-dimension.
101 It requires some set-up (initialization) phase and can be very fast when the distribution
102 parameters are not changed for every call.
103 
104 The following table shows some timings (in nanosecond/call)
105 for basic functions, TF1 functions and using UNURAN obtained running
106 the tutorial math/testrandom.C
107 Numbers have been obtained on an Intel Xeon Quad-core Harpertown (E5410) 2.33 GHz running
108 Linux SLC4 64 bit and compiled with gcc 3.4
109 
110 ~~~~
111 Distribution nanoseconds/call
112  TRandom TRandom1 TRandom2 TRandom3
113 Rndm.............. 5.000 105.000 7.000 10.000
114 RndmArray......... 4.000 104.000 6.000 9.000
115 Gaus.............. 36.000 180.000 40.000 48.000
116 Rannor............ 118.000 220.000 120.000 124.000
117 Landau............ 22.000 123.000 26.000 31.000
118 Exponential....... 93.000 198.000 98.000 104.000
119 Binomial(5,0.5)... 30.000 548.000 46.000 65.000
120 Binomial(15,0.5).. 75.000 1615.000 125.000 178.000
121 Poisson(3)........ 96.000 494.000 109.000 125.000
122 Poisson(10)....... 138.000 1236.000 165.000 203.000
123 Poisson(70)....... 818.000 1195.000 835.000 844.000
124 Poisson(100)...... 837.000 1218.000 849.000 864.000
125 GausTF1........... 83.000 180.000 87.000 88.000
126 LandauTF1......... 80.000 180.000 83.000 86.000
127 GausUNURAN........ 40.000 139.000 41.000 44.000
128 PoissonUNURAN(10). 85.000 271.000 92.000 102.000
129 PoissonUNURAN(100) 62.000 256.000 69.000 78.000
130 ~~~~
131 
132 Note that the time to generate a number from an arbitrary TF1 function
133 using TF1::GetRandom or using TUnuran is independent of the complexity of the function.
134 
135 TH1::FillRandom(TH1 *) or TH1::FillRandom(const char *tf1name)
136 can be used to fill an histogram (1-d, 2-d, 3-d from an existing histogram
137 or from an existing function.
138 
139 Note this interesting feature when working with objects.
140  You can use several TRandom objects, each with their "independent"
141  random sequence. For example, one can imagine
142 ~~~~
143  TRandom *eventGenerator = new TRandom();
144  TRandom *tracking = new TRandom();
145 ~~~~
146  `eventGenerator` can be used to generate the event kinematics.
147  tracking can be used to track the generated particles with random numbers
148  independent from eventGenerator.
149  This very interesting feature gives the possibility to work with simple
150  and very fast random number generators without worrying about
151  random number periodicity as it was the case with Fortran.
152  One can use TRandom::SetSeed to modify the seed of one generator.
153 
154 A TRandom object may be written to a Root file
155 
156 - as part of another object
157 - or with its own key (example: `gRandom->Write("Random")` ) ;
158 
159 */
160 
161 #include "TROOT.h"
162 #include "TMath.h"
163 #include "TRandom.h"
164 #include "TRandom3.h"
165 #include "TSystem.h"
166 #include "TDirectory.h"
167 #include "Math/QuantFuncMathCore.h"
168 #include "TUUID.h"
169 
170 ClassImp(TRandom);
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Default constructor. For seed see SetSeed().
174 
175 TRandom::TRandom(UInt_t seed): TNamed("Random","Default Random number generator")
176 {
177  SetSeed(seed);
178 }
179 
180 ////////////////////////////////////////////////////////////////////////////////
181 /// Default destructor. Can reset gRandom to 0 if gRandom points to this
182 /// generator.
183 
184 TRandom::~TRandom()
185 {
186  if (gRandom == this) gRandom = 0;
187 }
188 
189 ////////////////////////////////////////////////////////////////////////////////
190 /// Generates a random integer N according to the binomial law.
191 /// Coded from Los Alamos report LA-5061-MS.
192 ///
193 /// N is binomially distributed between 0 and ntot inclusive
194 /// with mean prob*ntot and prob is between 0 and 1.
195 ///
196 /// Note: This function should not be used when ntot is large (say >100).
197 /// The normal approximation is then recommended instead
198 /// (with mean =*ntot+0.5 and standard deviation sqrt(ntot*prob*(1-prob)).
199 
200 Int_t TRandom::Binomial(Int_t ntot, Double_t prob)
201 {
202  if (prob < 0 || prob > 1) return 0;
203  Int_t n = 0;
204  for (Int_t i=0;i<ntot;i++) {
205  if (Rndm() > prob) continue;
206  n++;
207  }
208  return n;
209 }
210 
211 ////////////////////////////////////////////////////////////////////////////////
212 /// Return a number distributed following a BreitWigner function with mean and gamma.
213 
214 Double_t TRandom::BreitWigner(Double_t mean, Double_t gamma)
215 {
216  Double_t rval, displ;
217  rval = 2*Rndm() - 1;
218  displ = 0.5*gamma*TMath::Tan(rval*TMath::PiOver2());
219 
220  return (mean+displ);
221 }
222 
223 ////////////////////////////////////////////////////////////////////////////////
224 /// Generates random vectors, uniformly distributed over a circle of given radius.
225 /// Input : r = circle radius
226 /// Output: x,y a random 2-d vector of length r
227 
228 void TRandom::Circle(Double_t &x, Double_t &y, Double_t r)
229 {
230  Double_t phi = Uniform(0,TMath::TwoPi());
231  x = r*TMath::Cos(phi);
232  y = r*TMath::Sin(phi);
233 }
234 
235 ////////////////////////////////////////////////////////////////////////////////
236 /// Returns an exponential deviate.
237 ///
238 /// exp( -t/tau )
239 
240 Double_t TRandom::Exp(Double_t tau)
241 {
242  Double_t x = Rndm(); // uniform on ] 0, 1 ]
243  Double_t t = -tau * TMath::Log( x ); // convert to exponential distribution
244  return t;
245 }
246 
247 ////////////////////////////////////////////////////////////////////////////////
248 /// Samples a random number from the standard Normal (Gaussian) Distribution
249 /// with the given mean and sigma.
250 /// Uses the Acceptance-complement ratio from W. Hoermann and G. Derflinger
251 /// This is one of the fastest existing method for generating normal random variables.
252 /// It is a factor 2/3 faster than the polar (Box-Muller) method used in the previous
253 /// version of TRandom::Gaus. The speed is comparable to the Ziggurat method (from Marsaglia)
254 /// implemented for example in GSL and available in the MathMore library.
255 ///
256 /// REFERENCE: - W. Hoermann and G. Derflinger (1990):
257 /// The ACR Method for generating normal random variables,
258 /// OR Spektrum 12 (1990), 181-185.
259 ///
260 /// Implementation taken from
261 /// UNURAN (c) 2000 W. Hoermann & J. Leydold, Institut f. Statistik, WU Wien
262 
263 Double_t TRandom::Gaus(Double_t mean, Double_t sigma)
264 {
265  const Double_t kC1 = 1.448242853;
266  const Double_t kC2 = 3.307147487;
267  const Double_t kC3 = 1.46754004;
268  const Double_t kD1 = 1.036467755;
269  const Double_t kD2 = 5.295844968;
270  const Double_t kD3 = 3.631288474;
271  const Double_t kHm = 0.483941449;
272  const Double_t kZm = 0.107981933;
273  const Double_t kHp = 4.132731354;
274  const Double_t kZp = 18.52161694;
275  const Double_t kPhln = 0.4515827053;
276  const Double_t kHm1 = 0.516058551;
277  const Double_t kHp1 = 3.132731354;
278  const Double_t kHzm = 0.375959516;
279  const Double_t kHzmp = 0.591923442;
280  /*zhm 0.967882898*/
281 
282  const Double_t kAs = 0.8853395638;
283  const Double_t kBs = 0.2452635696;
284  const Double_t kCs = 0.2770276848;
285  const Double_t kB = 0.5029324303;
286  const Double_t kX0 = 0.4571828819;
287  const Double_t kYm = 0.187308492 ;
288  const Double_t kS = 0.7270572718 ;
289  const Double_t kT = 0.03895759111;
290 
291  Double_t result;
292  Double_t rn,x,y,z;
293 
294  do {
295  y = Rndm();
296 
297  if (y>kHm1) {
298  result = kHp*y-kHp1; break; }
299 
300  else if (y<kZm) {
301  rn = kZp*y-1;
302  result = (rn>0) ? (1+rn) : (-1+rn);
303  break;
304  }
305 
306  else if (y<kHm) {
307  rn = Rndm();
308  rn = rn-1+rn;
309  z = (rn>0) ? 2-rn : -2-rn;
310  if ((kC1-y)*(kC3+TMath::Abs(z))<kC2) {
311  result = z; break; }
312  else {
313  x = rn*rn;
314  if ((y+kD1)*(kD3+x)<kD2) {
315  result = rn; break; }
316  else if (kHzmp-y<exp(-(z*z+kPhln)/2)) {
317  result = z; break; }
318  else if (y+kHzm<exp(-(x+kPhln)/2)) {
319  result = rn; break; }
320  }
321  }
322 
323  while (1) {
324  x = Rndm();
325  y = kYm * Rndm();
326  z = kX0 - kS*x - y;
327  if (z>0)
328  rn = 2+y/x;
329  else {
330  x = 1-x;
331  y = kYm-y;
332  rn = -(2+y/x);
333  }
334  if ((y-kAs+x)*(kCs+x)+kBs<0) {
335  result = rn; break; }
336  else if (y<x+kT)
337  if (rn*rn<4*(kB-log(x))) {
338  result = rn; break; }
339  }
340  } while(0);
341 
342  return mean + sigma * result;
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// Returns a random integer uniformly distributed on the interval [ 0, imax-1 ].
347 /// Note that the interfal contains the values of 0 and imax-1 but not imax.
348 
349 UInt_t TRandom::Integer(UInt_t imax)
350 {
351  UInt_t ui;
352  ui = (UInt_t)(imax*Rndm());
353  return ui;
354 }
355 
356 ////////////////////////////////////////////////////////////////////////////////
357 /// Generate a random number following a Landau distribution
358 /// with location parameter mu and scale parameter sigma:
359 /// Landau( (x-mu)/sigma )
360 /// Note that mu is not the mpv(most probable value) of the Landa distribution
361 /// and sigma is not the standard deviation of the distribution which is not defined.
362 /// For mu =0 and sigma=1, the mpv = -0.22278
363 ///
364 /// The Landau random number generation is implemented using the
365 /// function landau_quantile(x,sigma), which provides
366 /// the inverse of the landau cumulative distribution.
367 /// landau_quantile has been converted from CERNLIB ranlan(G110).
368 
369 Double_t TRandom::Landau(Double_t mu, Double_t sigma)
370 {
371  if (sigma <= 0) return 0;
372  Double_t x = Rndm();
373  Double_t res = mu + ROOT::Math::landau_quantile(x, sigma);
374  return res;
375 }
376 
377 ////////////////////////////////////////////////////////////////////////////////
378 /// Generates a random integer N according to a Poisson law.
379 /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
380 ///
381 /// Use a different procedure according to the mean value.
382 /// The algorithm is the same used by CLHEP.
383 /// For lower value (mean < 25) use the rejection method based on
384 /// the exponential.
385 /// For higher values use a rejection method comparing with a Lorentzian
386 /// distribution, as suggested by several authors.
387 /// This routine since is returning 32 bits integer will not work for values
388 /// larger than 2*10**9.
389 /// One should then use the Trandom::PoissonD for such large values.
390 
391 Int_t TRandom::Poisson(Double_t mean)
392 {
393  Int_t n;
394  if (mean <= 0) return 0;
395  if (mean < 25) {
396  Double_t expmean = TMath::Exp(-mean);
397  Double_t pir = 1;
398  n = -1;
399  while(1) {
400  n++;
401  pir *= Rndm();
402  if (pir <= expmean) break;
403  }
404  return n;
405  }
406  // for large value we use inversion method
407  else if (mean < 1E9) {
408  Double_t em, t, y;
409  Double_t sq, alxm, g;
410  Double_t pi = TMath::Pi();
411 
412  sq = TMath::Sqrt(2.0*mean);
413  alxm = TMath::Log(mean);
414  g = mean*alxm - TMath::LnGamma(mean + 1.0);
415 
416  do {
417  do {
418  y = TMath::Tan(pi*Rndm());
419  em = sq*y + mean;
420  } while( em < 0.0 );
421 
422  em = TMath::Floor(em);
423  t = 0.9*(1.0 + y*y)* TMath::Exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
424  } while( Rndm() > t );
425 
426  return static_cast<Int_t> (em);
427 
428  }
429  else {
430  // use Gaussian approximation vor very large values
431  n = Int_t(Gaus(0,1)*TMath::Sqrt(mean) + mean +0.5);
432  return n;
433  }
434 }
435 
436 ////////////////////////////////////////////////////////////////////////////////
437 /// Generates a random number according to a Poisson law.
438 /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
439 ///
440 /// This function is a variant of TRandom::Poisson returning a double
441 /// instead of an integer.
442 
443 Double_t TRandom::PoissonD(Double_t mean)
444 {
445  Int_t n;
446  if (mean <= 0) return 0;
447  if (mean < 25) {
448  Double_t expmean = TMath::Exp(-mean);
449  Double_t pir = 1;
450  n = -1;
451  while(1) {
452  n++;
453  pir *= Rndm();
454  if (pir <= expmean) break;
455  }
456  return static_cast<Double_t>(n);
457  }
458  // for large value we use inversion method
459  else if (mean < 1E9) {
460  Double_t em, t, y;
461  Double_t sq, alxm, g;
462  Double_t pi = TMath::Pi();
463 
464  sq = TMath::Sqrt(2.0*mean);
465  alxm = TMath::Log(mean);
466  g = mean*alxm - TMath::LnGamma(mean + 1.0);
467 
468  do {
469  do {
470  y = TMath::Tan(pi*Rndm());
471  em = sq*y + mean;
472  } while( em < 0.0 );
473 
474  em = TMath::Floor(em);
475  t = 0.9*(1.0 + y*y)* TMath::Exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
476  } while( Rndm() > t );
477 
478  return em;
479 
480  } else {
481  // use Gaussian approximation vor very large values
482  return Gaus(0,1)*TMath::Sqrt(mean) + mean +0.5;
483  }
484 }
485 
486 ////////////////////////////////////////////////////////////////////////////////
487 /// Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
488 
489 void TRandom::Rannor(Float_t &a, Float_t &b)
490 {
491  Double_t r, x, y, z;
492 
493  y = Rndm();
494  z = Rndm();
495  x = z * 6.28318530717958623;
496  r = TMath::Sqrt(-2*TMath::Log(y));
497  a = (Float_t)(r * TMath::Sin(x));
498  b = (Float_t)(r * TMath::Cos(x));
499 }
500 
501 ////////////////////////////////////////////////////////////////////////////////
502 /// Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
503 
504 void TRandom::Rannor(Double_t &a, Double_t &b)
505 {
506  Double_t r, x, y, z;
507 
508  y = Rndm();
509  z = Rndm();
510  x = z * 6.28318530717958623;
511  r = TMath::Sqrt(-2*TMath::Log(y));
512  a = r * TMath::Sin(x);
513  b = r * TMath::Cos(x);
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// Reads saved random generator status from filename.
518 
519 void TRandom::ReadRandom(const char *filename)
520 {
521  if (!gDirectory) return;
522  char *fntmp = gSystem->ExpandPathName(filename);
523  TDirectory *file = (TDirectory*)gROOT->ProcessLine(Form("TFile::Open(\"%s\");",fntmp));
524  delete [] fntmp;
525  if(file && file->GetFile()) {
526  gDirectory->ReadTObject(this,GetName());
527  delete file;
528  }
529 }
530 
531 ////////////////////////////////////////////////////////////////////////////////
532 /// Machine independent random number generator.
533 /// Based on the BSD Unix (Rand) Linear congrential generator.
534 /// Produces uniformly-distributed floating points between 0 and 1.
535 /// Identical sequence on all machines of >= 32 bits.
536 /// Periodicity = 2**31, generates a number in (0,1).
537 /// Note that this is a generator which is known to have defects
538 /// (the lower random bits are correlated) and therefore should NOT be
539 /// used in any statistical study).
540 
541 Double_t TRandom::Rndm( )
542 {
543 #ifdef OLD_TRANDOM_IMPL
544  const Double_t kCONS = 4.6566128730774E-10;
545  const Int_t kMASK24 = 2147483392;
546 
547  fSeed *= 69069;
548  UInt_t jy = (fSeed&kMASK24); // Set lower 8 bits to zero to assure exact float
549  if (jy) return kCONS*jy;
550  return Rndm();
551 #endif
552 
553  // kCONS = 1./2147483648 = 1./(RAND_MAX+1) and RAND_MAX= 0x7fffffffUL
554  const Double_t kCONS = 4.6566128730774E-10; // (1/pow(2,31)
555  fSeed = (1103515245 * fSeed + 12345) & 0x7fffffffUL;
556 
557  if (fSeed) return kCONS*fSeed;
558  return Rndm();
559 }
560 
561 ////////////////////////////////////////////////////////////////////////////////
562 /// Return an array of n random numbers uniformly distributed in ]0,1].
563 
564 void TRandom::RndmArray(Int_t n, Double_t *array)
565 {
566  const Double_t kCONS = 4.6566128730774E-10; // (1/pow(2,31))
567  Int_t i=0;
568  while (i<n) {
569  fSeed = (1103515245 * fSeed + 12345) & 0x7fffffffUL;
570  if (fSeed) {array[i] = kCONS*fSeed; i++;}
571  }
572 }
573 
574 ////////////////////////////////////////////////////////////////////////////////
575 /// Return an array of n random numbers uniformly distributed in ]0,1].
576 
577 void TRandom::RndmArray(Int_t n, Float_t *array)
578 {
579  const Double_t kCONS = 4.6566128730774E-10; // (1/pow(2,31))
580  Int_t i=0;
581  while (i<n) {
582  fSeed = (1103515245 * fSeed + 12345) & 0x7fffffffUL;
583  if (fSeed) {array[i] = Float_t(kCONS*fSeed); i++;}
584  }
585 }
586 
587 ////////////////////////////////////////////////////////////////////////////////
588 /// Set the random generator seed. Note that default value is zero, which is
589 /// different than the default value used when constructing the class.
590 /// If the seed is zero the seed is set to a random value
591 /// which in case of TRandom depends on the lowest 4 bytes of TUUID
592 /// The UUID will be identical if SetSeed(0) is called with time smaller than 100 ns
593 /// Instead if a different generator implementation is used (TRandom1, 2 or 3)
594 /// the seed is generated using a 128 bit UUID. This results in different seeds
595 /// and then random sequence for every SetSeed(0) call.
596 
597 void TRandom::SetSeed(ULong_t seed)
598 {
599  if( seed==0 ) {
600  TUUID u;
601  UChar_t uuid[16];
602  u.GetUUID(uuid);
603  fSeed = UInt_t(uuid[3])*16777216 + UInt_t(uuid[2])*65536 + UInt_t(uuid[1])*256 + UInt_t(uuid[0]);
604  } else {
605  fSeed = seed;
606  }
607 }
608 
609 ////////////////////////////////////////////////////////////////////////////////
610 /// Generates random vectors, uniformly distributed over the surface
611 /// of a sphere of given radius.
612 /// Input : r = sphere radius
613 /// Output: x,y,z a random 3-d vector of length r
614 /// Method: (based on algorithm suggested by Knuth and attributed to Robert E Knop)
615 /// which uses less random numbers than the CERNLIB RN23DIM algorithm
616 
617 void TRandom::Sphere(Double_t &x, Double_t &y, Double_t &z, Double_t r)
618 {
619  Double_t a=0,b=0,r2=1;
620  while (r2 > 0.25) {
621  a = Rndm() - 0.5;
622  b = Rndm() - 0.5;
623  r2 = a*a + b*b;
624  }
625  z = r* ( -1. + 8.0 * r2 );
626 
627  Double_t scale = 8.0 * r * TMath::Sqrt(0.25 - r2);
628  x = a*scale;
629  y = b*scale;
630 }
631 
632 ////////////////////////////////////////////////////////////////////////////////
633 /// Returns a uniform deviate on the interval (0, x1).
634 
635 Double_t TRandom::Uniform(Double_t x1)
636 {
637  Double_t ans = Rndm();
638  return x1*ans;
639 }
640 
641 ////////////////////////////////////////////////////////////////////////////////
642 /// Returns a uniform deviate on the interval (x1, x2).
643 
644 Double_t TRandom::Uniform(Double_t x1, Double_t x2)
645 {
646  Double_t ans= Rndm();
647  return x1 + (x2-x1)*ans;
648 }
649 
650 ////////////////////////////////////////////////////////////////////////////////
651 /// Writes random generator status to filename.
652 
653 void TRandom::WriteRandom(const char *filename) const
654 {
655  if (!gDirectory) return;
656  char *fntmp = gSystem->ExpandPathName(filename);
657  TDirectory *file = (TDirectory*)gROOT->ProcessLine(Form("TFile::Open(\"%s\",\"recreate\");",fntmp));
658  delete [] fntmp;
659  if(file && file->GetFile()) {
660  gDirectory->WriteTObject(this,GetName());
661  delete file;
662  }
663 }