Logo ROOT   6.12/07
Reference Guide
RandomFunctions.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: L. Moneta 8/2015
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2015 , ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // file for random class
12 //
13 //
14 // Created by: Lorenzo Moneta : Tue 4 Aug 2015
15 //
16 //
17 #include "Math/RandomFunctions.h"
18 
19 #include "Math/DistFunc.h"
20 
21 #include "TMath.h"
22 
23 namespace ROOT {
24 namespace Math {
25 
26 
28 {
29  if (prob < 0 || prob > 1) return 0;
30  Int_t n = 0;
31  for (Int_t i=0;i<ntot;i++) {
32  if (Rndm() > prob) continue;
33  n++;
34  }
35  return n;
36 }
37 
38  ////////////////////////////////////////////////////////////////////////////////
39 /// Return a number distributed following a BreitWigner function with mean and gamma.
40 
42 {
43  Double_t rval, displ;
44  rval = 2*Rndm() - 1;
45  displ = 0.5*gamma*TMath::Tan(rval*TMath::PiOver2());
46 
47  return (mean+displ);
48 }
49 
50 ////////////////////////////////////////////////////////////////////////////////
51 /// Generates random vectors, uniformly distributed over a circle of given radius.
52 /// Input : r = circle radius
53 /// Output: x,y a random 2-d vector of length r
54 
56 {
57  Double_t phi = Uniform(0,TMath::TwoPi());
58  x = r*TMath::Cos(phi);
59  y = r*TMath::Sin(phi);
60 }
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// Returns an exponential deviate.
64 ///
65 /// exp( -t/tau )
66 
68 {
69  Double_t x = Rndm(); // uniform on ] 0, 1 ]
70  Double_t t = -tau * TMath::Log( x ); // convert to exponential distribution
71  return t;
72 }
73 
74 
75 
76 
78  double y = Rndm();
79  double z = Rndm();
80  double x = z * 6.28318530717958623;
81  double radius = std::sqrt(-2*std::log(y));
82  double g = radius * std::sin(x);
83  return mean + g * sigma;
84 }
85 
86 
87  // double GausImpl(TRandomEngine * r, double mean, double sigma) {
88  // double y = r->Rndm();
89  // double z = r->Rndm();
90  // double x = z * 6.28318530717958623;
91  // double radius = std::sqrt(-2*std::log(y));
92  // double g = radius * std::sin(x);
93  // return mean + g * sigma;
94  // }
95 
97 {
98  const Double_t kC1 = 1.448242853;
99  const Double_t kC2 = 3.307147487;
100  const Double_t kC3 = 1.46754004;
101  const Double_t kD1 = 1.036467755;
102  const Double_t kD2 = 5.295844968;
103  const Double_t kD3 = 3.631288474;
104  const Double_t kHm = 0.483941449;
105  const Double_t kZm = 0.107981933;
106  const Double_t kHp = 4.132731354;
107  const Double_t kZp = 18.52161694;
108  const Double_t kPhln = 0.4515827053;
109  const Double_t kHm1 = 0.516058551;
110  const Double_t kHp1 = 3.132731354;
111  const Double_t kHzm = 0.375959516;
112  const Double_t kHzmp = 0.591923442;
113  /*zhm 0.967882898*/
114 
115  const Double_t kAs = 0.8853395638;
116  const Double_t kBs = 0.2452635696;
117  const Double_t kCs = 0.2770276848;
118  const Double_t kB = 0.5029324303;
119  const Double_t kX0 = 0.4571828819;
120  const Double_t kYm = 0.187308492 ;
121  const Double_t kS = 0.7270572718 ;
122  const Double_t kT = 0.03895759111;
123 
124  Double_t result;
125  Double_t rn,x,y,z;
126 
127  do {
128  y = Rndm();
129 
130  if (y>kHm1) {
131  result = kHp*y-kHp1; break; }
132 
133  else if (y<kZm) {
134  rn = kZp*y-1;
135  result = (rn>0) ? (1+rn) : (-1+rn);
136  break;
137  }
138 
139  else if (y<kHm) {
140  rn = Rndm();
141  rn = rn-1+rn;
142  z = (rn>0) ? 2-rn : -2-rn;
143  if ((kC1-y)*(kC3+TMath::Abs(z))<kC2) {
144  result = z; break; }
145  else {
146  x = rn*rn;
147  if ((y+kD1)*(kD3+x)<kD2) {
148  result = rn; break; }
149  else if (kHzmp-y<exp(-(z*z+kPhln)/2)) {
150  result = z; break; }
151  else if (y+kHzm<exp(-(x+kPhln)/2)) {
152  result = rn; break; }
153  }
154  }
155 
156  while (1) {
157  x = Rndm();
158  y = kYm * Rndm();
159  z = kX0 - kS*x - y;
160  if (z>0)
161  rn = 2+y/x;
162  else {
163  x = 1-x;
164  y = kYm-y;
165  rn = -(2+y/x);
166  }
167  if ((y-kAs+x)*(kCs+x)+kBs<0) {
168  result = rn; break; }
169  else if (y<x+kT)
170  if (rn*rn<4*(kB-log(x))) {
171  result = rn; break; }
172  }
173  } while(0);
174 
175  return mean + sigma * result;
176 }
177 
178 ////////////////////////////////////////////////////////////////////////////////
179 /// Generate a random number following a Landau distribution
180 /// with location parameter mu and scale parameter sigma:
181 /// Landau( (x-mu)/sigma )
182 /// Note that mu is not the mpv(most probable value) of the Landa distribution
183 /// and sigma is not the standard deviation of the distribution which is not defined.
184 /// For mu =0 and sigma=1, the mpv = -0.22278
185 ///
186 /// The Landau random number generation is implemented using the
187 /// function landau_quantile(x,sigma), which provides
188 /// the inverse of the landau cumulative distribution.
189 /// landau_quantile has been converted from CERNLIB ranlan(G110).
190 
192 {
193  if (sigma <= 0) return 0;
194  Double_t x = Rndm();
195  Double_t res = mu + ROOT::Math::landau_quantile(x, sigma);
196  return res;
197 }
198 
199 ////////////////////////////////////////////////////////////////////////////////
200 /// Generates a random integer N according to a Poisson law.
201 /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
202 ///
203 /// Use a different procedure according to the mean value.
204 /// The algorithm is the same used by CLHEP.
205 /// For lower value (mean < 25) use the rejection method based on
206 /// the exponential.
207 /// For higher values use a rejection method comparing with a Lorentzian
208 /// distribution, as suggested by several authors.
209 /// This routine since is returning 32 bits integer will not work for values
210 /// larger than 2*10**9.
211 /// One should then use the Trandom::PoissonD for such large values.
212 
214 {
215  Int_t n;
216  if (mean <= 0) return 0;
217  if (mean < 25) {
218  Double_t expmean = TMath::Exp(-mean);
219  Double_t pir = 1;
220  n = -1;
221  while(1) {
222  n++;
223  pir *= Rndm();
224  if (pir <= expmean) break;
225  }
226  return n;
227  }
228  // for large value we use inversion method
229  else if (mean < 1E9) {
230  Double_t em, t, y;
231  Double_t sq, alxm, g;
232  Double_t pi = TMath::Pi();
233 
234  sq = TMath::Sqrt(2.0*mean);
235  alxm = TMath::Log(mean);
236  g = mean*alxm - TMath::LnGamma(mean + 1.0);
237 
238  do {
239  do {
240  y = TMath::Tan(pi*Rndm());
241  em = sq*y + mean;
242  } while( em < 0.0 );
243 
244  em = TMath::Floor(em);
245  t = 0.9*(1.0 + y*y)* TMath::Exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
246  } while( Rndm() > t );
247 
248  return static_cast<Int_t> (em);
249 
250  }
251  else {
252  // use Gaussian approximation vor very large values
253  n = Int_t(Gaus(0,1)*TMath::Sqrt(mean) + mean +0.5);
254  return n;
255  }
256 }
257 
258 ////////////////////////////////////////////////////////////////////////////////
259 /// Generates a random number according to a Poisson law.
260 /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
261 ///
262 /// This function is a variant of RandomFunctionsImpl<TRandomEngine>::Poisson returning a double
263 /// instead of an integer.
264 
266 {
267  Int_t n;
268  if (mean <= 0) return 0;
269  if (mean < 25) {
270  Double_t expmean = TMath::Exp(-mean);
271  Double_t pir = 1;
272  n = -1;
273  while(1) {
274  n++;
275  pir *= Rndm();
276  if (pir <= expmean) break;
277  }
278  return static_cast<Double_t>(n);
279  }
280  // for large value we use inversion method
281  else if (mean < 1E9) {
282  Double_t em, t, y;
283  Double_t sq, alxm, g;
284  Double_t pi = TMath::Pi();
285 
286  sq = TMath::Sqrt(2.0*mean);
287  alxm = TMath::Log(mean);
288  g = mean*alxm - TMath::LnGamma(mean + 1.0);
289 
290  do {
291  do {
292  y = TMath::Tan(pi*Rndm());
293  em = sq*y + mean;
294  } while( em < 0.0 );
295 
296  em = TMath::Floor(em);
297  t = 0.9*(1.0 + y*y)* TMath::Exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
298  } while( Rndm() > t );
299 
300  return em;
301 
302  } else {
303  // use Gaussian approximation vor very large values
304  return Gaus(0,1)*TMath::Sqrt(mean) + mean +0.5;
305  }
306 }
307 
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
311 
313 {
314  Double_t r, x, y, z;
315 
316  y = Rndm();
317  z = Rndm();
318  x = z * 6.28318530717958623;
319  r = TMath::Sqrt(-2*TMath::Log(y));
320  a = r * TMath::Sin(x);
321  b = r * TMath::Cos(x);
322 }
323 
324 ////////////////////////////////////////////////////////////////////////////////
325 /// Generates random vectors, uniformly distributed over the surface
326 /// of a sphere of given radius.
327 /// Input : r = sphere radius
328 /// Output: x,y,z a random 3-d vector of length r
329 /// Method: (based on algorithm suggested by Knuth and attributed to Robert E Knop)
330 /// which uses less random numbers than the CERNLIB RN23DIM algorithm
331 
333 {
334  Double_t a=0,b=0,r2=1;
335  while (r2 > 0.25) {
336  a = Rndm() - 0.5;
337  b = Rndm() - 0.5;
338  r2 = a*a + b*b;
339  }
340  z = r* ( -1. + 8.0 * r2 );
341 
342  Double_t scale = 8.0 * r * TMath::Sqrt(0.25 - r2);
343  x = a*scale;
344  y = b*scale;
345 }
346 
347 ////////////////////////////////////////////////////////////////////////////////
348 /// Returns a uniform deviate on the interval (0, x1).
349 
351 {
352  Double_t ans = Rndm();
353  return x1*ans;
354 }
355 
356 ////////////////////////////////////////////////////////////////////////////////
357 /// Returns a uniform deviate on the interval (x1, x2).
358 
360  return (b-a) * Rndm() + a;
361 }
362 
363 
364  } // namespace Math
365 } // namespace ROOT
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculate a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:441
Double_t Floor(Double_t x)
Definition: TMath.h:599
static constexpr double pi
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
Double_t Log(Double_t x)
Definition: TMath.h:648
constexpr Double_t TwoPi()
Definition: TMath.h:44
Double_t Landau(Double_t x, Double_t mean, Double_t sigma)
Double_t Gaus(Double_t x, Double_t mean, Double_t sigma)
Gauss.
int Int_t
Definition: RtypesCore.h:41
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
you should not use this method at all Int_t Int_t Double_t Double_t em
Definition: TRolke.cxx:630
const Double_t sigma
constexpr Double_t Pi()
Definition: TMath.h:40
double sin(double)
double gamma(double x)
ROOT::R::TRInterface & r
Definition: Object.C:4
auto * a
Definition: textangle.C:12
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2079
Double_t Cos(Double_t)
Definition: TMath.h:550
Double_t Exp(Double_t x)
Definition: TMath.h:621
static const double x1[5]
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
Definition of the generic impelmentation class for the RandomFunctions.
Double_t Sin(Double_t)
Definition: TMath.h:547
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:489
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
double exp(double)
const Int_t n
Definition: legend1.C:16
Double_t Tan(Double_t)
Definition: TMath.h:553
constexpr Double_t PiOver2()
Definition: TMath.h:48
double log(double)
static constexpr double g