Logo ROOT  
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//
18
20
21#include "TMath.h"
22
23namespace ROOT {
24namespace 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();
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;
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;
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
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define g(i)
Definition: RSha256.hxx:105
static const double x1[5]
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
double sqrt(double)
double sin(double)
double exp(double)
double log(double)
Definition of the generic impelmentation class for the RandomFunctions.
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double gamma(double x)
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.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
static constexpr double pi
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2083
Double_t Exp(Double_t x)
Definition: TMath.h:717
Double_t Floor(Double_t x)
Definition: TMath.h:693
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:437
constexpr Double_t PiOver2()
Definition: TMath.h:52
Double_t Log(Double_t x)
Definition: TMath.h:750
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Double_t Cos(Double_t)
Definition: TMath.h:631
constexpr Double_t Pi()
Definition: TMath.h:38
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:486
Double_t Sin(Double_t)
Definition: TMath.h:627
Double_t Tan(Double_t)
Definition: TMath.h:635
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
constexpr Double_t TwoPi()
Definition: TMath.h:45
auto * a
Definition: textangle.C:12