]>
Commit | Line | Data |
---|---|---|
4c039060 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
f531a546 | 16 | // $Id$ |
4c039060 | 17 | |
959fbac5 | 18 | /////////////////////////////////////////////////////////////////////////// |
19 | // Class AliRandom | |
20 | // Generate universal random numbers on all common machines. | |
21 | // Available distributions : Uniform, Gaussian, Poisson and | |
22 | // User defined function | |
23 | // | |
24 | // Features : | |
25 | // ---------- | |
26 | // 1) Period = 2**144 | |
27 | // 2) Same sequence of 24-bit real numbers on all common machines | |
28 | // | |
29 | // Reference : | |
30 | // ----------- | |
31 | // G.Marsaglia and A.Zaman, FSU-SCRI-87-50, Florida State University, 1987. | |
32 | // | |
33 | // Coding example : | |
34 | // ---------------- | |
35 | // | |
36 | // Float_t rndm; // Variable to hold a single random number | |
37 | // const Int_t n=1000; | |
38 | // Float_t rvec[n]; // Vector to hold n random numbers | |
39 | // | |
40 | // AliRandom r; // Create a Random object with default sequence | |
41 | // | |
42 | // rndm=r.Uniform(); // Provide a uniform random number in <0,1> | |
43 | // Float_t a=3.; | |
44 | // Float_t b=5.; | |
45 | // rndm=r.Uniform(a,b); // Provide a uniform random number in <a,b> | |
46 | // r.Uniform(rvec,n); // Provide n uniform randoms in <0,1> in rvec | |
47 | // r.Uniform(rvec,n,a,b); // Provide n uniform randoms in <a,b> in rvec | |
48 | // | |
49 | // rndm=r.Gauss(); // Provide a Gaussian random number with | |
50 | // // mean=0 and sigma=1 | |
51 | // Float_t mean=25.; | |
52 | // Float_t sigma=5.; | |
53 | // rndm=r.Gauss(mean,sigma); // Provide a Gaussian random number | |
54 | // // with specified mean and sigma | |
55 | // r.Gauss(rvec,n); // n Gaussian randoms mean=0 sigma=1 | |
56 | // r.Gauss(rvec,n,mean,sigma); // n Gaussian randoms with specified | |
57 | // // mean and sigma | |
58 | // | |
59 | // rndm=r.Poisson(mean); // Provide a Poisson random number with | |
60 | // // specified mean | |
61 | // r.Poisson(rvec,nmean); // n Poisson randoms with specified mean | |
62 | // | |
63 | // Int_t seed=1837724 | |
64 | // AliRandom p(seed); // Create a Random object with specified seed. | |
65 | // // The sequence is started from scratch. | |
66 | // Int_t cnt1=25; | |
67 | // Int_t cnt2=8; | |
68 | // AliRandom q(seed,cnt1,cnt2); // Create a Random object with specified seed | |
69 | // // The sequence is started at the location | |
70 | // // denoted by the counters cnt1 and cnt2. | |
71 | // | |
84bb7c66 | 72 | // q.Data(); // Print the current seed, cnt1 and cnt2 values. |
959fbac5 | 73 | // q.GetSeed(); // Provide the current seed value. |
74 | // q.GetCnt1(); // Provide the current cnt1 value. | |
75 | // q.GetCnt2(); // Provide the current cnt2 value. | |
76 | // | |
77 | // Float_t udist(Float_t x) // A user defined distribution | |
78 | // { | |
79 | // return x*x-4.*x; | |
80 | // } | |
81 | // | |
82 | // Int_t nbins=100; | |
83 | // q.SetUser(a,b,nbins,udist); // Initialise generator for udist distribution | |
84 | // q.User(); // Provide a random number according to the udist distribution | |
85 | // q.User(rvec,n); // Provide n randoms according to the udist distribution | |
86 | // | |
87 | // Float_t* x=new Float_t[nbins]; | |
88 | // Float_t* y=new Float_t[nbins]; | |
89 | // | |
90 | // ... code to fill x[] and y[] .. | |
91 | // | |
92 | // AliRandom s; | |
93 | // s.SetUser(x,y,nbins); // Initialise generator for (x[i],y[i]) distribution | |
94 | // s.User(); // Provide a random number according to the user distribution | |
95 | // s.User(rvec,n); // Provide n randoms according to the user distribution | |
96 | // | |
97 | // Notes : | |
98 | // ------- | |
99 | // 1) Allowed seed values : 0 <= seed <= 921350143 | |
100 | // Default seed = 53310452 | |
101 | // 2) To ensure a unique sequence for each run, one can automatically | |
102 | // construct a seed value by e.g. using the date and time. | |
103 | // 3) Using the rvec facility saves a lot of CPU time for large n values. | |
104 | // | |
105 | //--- Author: Nick van Eijndhoven 11-oct-1997 UU-SAP Utrecht | |
f531a546 | 106 | //- Modified: NvE $Date$ UU-SAP Utrecht |
959fbac5 | 107 | /////////////////////////////////////////////////////////////////////////// |
108 | ||
d88f97cc | 109 | #include "AliRandom.h" |
c72198f1 | 110 | #include "Riostream.h" |
d88f97cc | 111 | |
112 | ClassImp(AliRandom) // Class implementation to enable ROOT I/O | |
113 | ||
114 | AliRandom::AliRandom() | |
115 | { | |
959fbac5 | 116 | // Creation of an AliRandom object and default initialisation. |
d88f97cc | 117 | // |
118 | // A seed is used to create the initial u[97] table. | |
119 | // This seed is converted into four startup parameters i j k and l | |
120 | // (see member function "unpack"). | |
121 | // | |
122 | // Suggested test values : i=12 j=34 k=56 l=78 (see article) | |
123 | // which corresponds to : seed = 53310452 | |
124 | ||
125 | Int_t seed=53310452; // Default seed | |
126 | Start(seed,0,0); // Start the sequence for this seed from scratch | |
127 | } | |
128 | /////////////////////////////////////////////////////////////////////////// | |
129 | AliRandom::AliRandom(Int_t seed) | |
130 | { | |
131 | // Creation of an AliRandom object and user defined initialisation | |
132 | ||
133 | Start(seed,0,0); // Start the sequence for this seed from scratch | |
134 | } | |
135 | /////////////////////////////////////////////////////////////////////////// | |
136 | AliRandom::AliRandom(Int_t seed,Int_t cnt1,Int_t cnt2) | |
137 | { | |
138 | // Creation of an AliRandom object and user defined initialisation | |
139 | // | |
140 | // seed is the seed to create the initial u[97] table. | |
141 | // The range of the seed is : 0 <= seed <= 921350143 | |
142 | // | |
143 | // cnt1 and cnt2 are the parameters for the counting system | |
144 | // to enable a start of the sequence at a certain point. | |
145 | // The current values of seed, cnt1 and cnt2 can be obtained | |
146 | // via the member functions "GetSeed", "GetCnt1" and "GetCnt2" resp. | |
147 | // To start from scratch one should select : cnt1=0 and cnt2=0 | |
148 | ||
149 | Start(seed,cnt1,cnt2); // Start the sequence from a user defined point | |
150 | } | |
151 | /////////////////////////////////////////////////////////////////////////// | |
152 | AliRandom::~AliRandom() | |
153 | { | |
154 | // Destructor to delete memory allocated for the area function arrays | |
155 | if (fXa) delete [] fXa; | |
156 | fXa=0; | |
157 | if (fYa) delete [] fYa; | |
158 | fYa=0; | |
159 | if (fIbins) delete [] fIbins; | |
160 | fIbins=0; | |
161 | } | |
162 | /////////////////////////////////////////////////////////////////////////// | |
163 | void AliRandom::Start(Int_t seed,Int_t cnt1,Int_t cnt2) | |
164 | { | |
165 | // Start a certain sequence from scratch or from a user defined point | |
166 | // | |
167 | // The algorithm to start from scratch is based on the routine RSTART | |
168 | // as described in the report by G.Marsaglia and A.Zaman | |
169 | // (FSU-SCRI-87-50 Florida State University 1987). | |
170 | // | |
171 | // seed is the seed to create the initial u[97] table. | |
172 | // This seed is converted into four startup parameters i j k and l | |
173 | // (see the member function "unpack"). | |
174 | // | |
175 | // The range of the seed is : 0 <= seed <= 921350143 | |
176 | // | |
177 | // Suggested test values : i=12 j=34 k=56 l=78 (see article) | |
178 | // which corresponds to : seed = 53310452 | |
179 | // | |
180 | // cnt1 and cnt2 are the parameters for the counting system | |
181 | // to enable a start of the sequence at a certain point. | |
182 | // The current values of seed, cnt1 and cnt2 can be obtained | |
183 | // via the member functions "GetSeed", "GetCnt1" and "GetCnt2" resp. | |
184 | // To start from scratch one should select : cnt1=0 and cnt2=0 | |
185 | ||
186 | // Reset the area function | |
187 | fNa=0; | |
188 | fXa=0; | |
189 | fYa=0; | |
190 | fIbins=0; | |
191 | ||
192 | // Clipping parameter to prevent overflow of the counting system | |
193 | fClip=1000000; | |
194 | ||
195 | // Set the lags for the Fibonacci sequence of the first part | |
196 | // The sequence is set to F(97,33,*) (see article) | |
197 | fI=97; | |
198 | fJ=33; | |
199 | ||
200 | // Unpack the seed value and determine i, j, k and l | |
201 | fSeed=seed; | |
202 | Int_t i,j,k,l; | |
203 | Unpack(seed,i,j,k,l); | |
204 | ||
205 | // Reset counters | |
206 | fCnt1=0; | |
207 | fCnt2=0; | |
208 | ||
209 | // Fill the starting table for the first part of the combination | |
210 | Float_t s,t; | |
211 | Int_t m; | |
212 | for (Int_t ii=0; ii<97; ii++) | |
213 | { | |
214 | s=0.; | |
215 | t=0.5; | |
216 | ||
217 | for (Int_t jj=1; jj<=24; jj++) | |
218 | { | |
219 | m=(((i*j)%179)*k)%179; | |
220 | i=j; | |
221 | j=k; | |
222 | k=m; | |
223 | l=((53*l)+1)%169; | |
224 | if ((l*m)%64 >= 32) s+=t; | |
225 | t=0.5*t; | |
226 | } | |
227 | fU[ii]=s; | |
228 | } | |
229 | ||
230 | // Initialise the second part of the combination | |
231 | fC=362436./16777216.; | |
232 | fCd=7654321./16777216.; | |
233 | fCm=16777213./16777216.; | |
234 | ||
235 | // Generate random numbers upto the user selected starting point | |
236 | // on basis of the counting system | |
237 | if (cnt1 > 0) Uniform(cnt1); | |
238 | if (cnt2 > 0) | |
239 | { | |
240 | for (Int_t n=1; n<=cnt2; n++) | |
241 | { | |
242 | Uniform(fClip); | |
243 | } | |
244 | } | |
245 | } | |
246 | /////////////////////////////////////////////////////////////////////////// | |
247 | void AliRandom::Unpack(Int_t seed,Int_t& i,Int_t& j,Int_t& k,Int_t& l) | |
248 | { | |
249 | // Unpack the seed into the four startup parameters i,j,k and l | |
250 | // | |
251 | // The range of the seed is : 0 <= seed <= 921350143 | |
252 | // | |
253 | // From the article : | |
254 | // The i,j and k values may be chosen in the interval : 1 <= n <= 178 | |
255 | // The l value may be in the interval : 0 <= l <= 168 | |
256 | // | |
257 | // Taking into account the period of the 3-lagged Fibonacci sequence | |
258 | // The following "bad" combinations have to be ruled out : | |
259 | // | |
260 | // i j k l period | |
261 | // 1 1 1 X 1 | |
262 | // 178 1 1 X 4 | |
263 | // 1 178 1 X 2 | |
264 | // 1 1 178 X 4 | |
265 | // 178 178 1 X 4 | |
266 | // 178 1 178 X 2 | |
267 | // 1 178 178 X 4 | |
268 | // 178 178 178 X 1 | |
269 | // | |
270 | // To rule these "bad" combinations out all together, we choose | |
271 | // the following allowed ranges : | |
272 | // The i,j and k values may be chosen in the interval : 2 <= n <= 177 | |
273 | // The l value may be in the interval : 0 <= l <= 168 | |
274 | // | |
275 | // and use the formula : | |
276 | // seed = (i-2)*176*176*169 + (j-2)*176*169 + (k-2)*169 + l | |
277 | ||
278 | if ((seed >= 0) && (seed <= 921350143)) // Check seed value | |
279 | { | |
280 | Int_t idum=seed; | |
281 | Int_t imin2=idum/(176*176*169); | |
282 | idum=idum%(176*176*169); | |
283 | Int_t jmin2=idum/(176*169); | |
284 | idum=idum%(176*169); | |
285 | Int_t kmin2=idum/169; | |
286 | ||
287 | i=imin2+2; | |
288 | j=jmin2+2; | |
289 | k=kmin2+2; | |
290 | l=seed%169; | |
291 | } | |
292 | else | |
293 | { | |
294 | cout << " *AliRandom::unpack()* Unallowed seed value encountered." | |
295 | << " seed = " << seed << endl; | |
296 | cout << " Seed will be set to default value." << endl; | |
297 | ||
298 | seed=53310452; // Default seed | |
299 | Start(seed,0,0); // Start the sequence for this seed from scratch | |
300 | } | |
301 | } | |
302 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 303 | Int_t AliRandom::GetSeed() const |
d88f97cc | 304 | { |
305 | // Provide the current seed value | |
306 | return fSeed; | |
307 | } | |
308 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 309 | Int_t AliRandom::GetCnt1() const |
d88f97cc | 310 | { |
311 | // Provide the current value of the counter cnt1 | |
312 | return fCnt1; | |
313 | } | |
314 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 315 | Int_t AliRandom::GetCnt2() const |
d88f97cc | 316 | { |
317 | // Provide the current value of the counter cnt2 | |
318 | return fCnt2; | |
319 | } | |
320 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 321 | void AliRandom::Data() const |
d88f97cc | 322 | { |
323 | // Print the current seed, cnt1 and cnt2 values | |
324 | cout << " *Random* seed = " << fSeed | |
325 | << " cnt1 = " << fCnt1 << " cnt2 = " << fCnt2 << endl; | |
326 | } | |
327 | /////////////////////////////////////////////////////////////////////////// | |
328 | Float_t AliRandom::Uniform() | |
329 | { | |
330 | // Generate uniform random numbers in the interval <0,1> | |
331 | // | |
332 | // The algorithm is based on lagged Fibonacci sequences (first part) | |
333 | // combined with a congruential method (second part) | |
334 | // as described in the report by G.Marsaglia and A.Zaman | |
335 | // (FSU-SCRI-87-50 Florida State University 1987). | |
336 | // | |
337 | // Features : | |
338 | // 1) Period = 2**144 | |
339 | // 2) Same sequence of 24-bit real numbers on all common machines | |
340 | ||
341 | // First part of the combination : F(97,33,*) (see article for explanation) | |
342 | Float_t unirnu=fU[fI-1]-fU[fJ-1]; | |
343 | if (unirnu < 0) unirnu+=1.; | |
344 | fU[fI-1]=unirnu; | |
345 | fI-=1; | |
346 | if (fI == 0) fI=97; | |
347 | fJ-=1; | |
348 | if (fJ == 0) fJ=97; | |
349 | ||
350 | // Second part of the combination (see article for explanation) | |
351 | fC-=fCd; | |
352 | if (fC < 0.) fC+=fCm; | |
353 | unirnu-=fC; | |
354 | if (unirnu < 0.) unirnu+=1.; | |
355 | ||
356 | // Update the counting system to enable sequence continuation | |
357 | // at an arbitrary starting position. | |
358 | // Two counters have been introduced to avoid overflow | |
359 | // fCnt1 : Counter which goes up to fClip | |
360 | // and is reset when fClip is reached | |
361 | // fCnt2 : Counts the number of times fClip has been reached | |
362 | fCnt1+=1; | |
363 | if (fCnt1 >= fClip) | |
364 | { | |
365 | fCnt1=0; | |
366 | fCnt2+=1; | |
367 | } | |
368 | ||
369 | if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range | |
370 | ||
371 | return unirnu; | |
372 | } | |
373 | /////////////////////////////////////////////////////////////////////////// | |
374 | Float_t AliRandom::Uniform(Float_t a,Float_t b) | |
375 | { | |
376 | // Generate uniform random numbers in the interval <a,b> | |
377 | Float_t rmin=a; | |
378 | if (a > b) rmin=b; | |
379 | ||
380 | Float_t rndm=Uniform(); | |
381 | rndm=rmin+fabs(rndm*(a-b)); | |
382 | ||
383 | return rndm; | |
384 | } | |
385 | /////////////////////////////////////////////////////////////////////////// | |
386 | void AliRandom::Uniform(Float_t* vec,Int_t n,Float_t a,Float_t b) | |
387 | { | |
388 | // Generate a vector of uniform random numbers in the interval <a,b> | |
389 | // This saves lots of (member)function calls in case many random | |
390 | // numbers are needed at once. | |
391 | // | |
392 | // n = The number of random numbers to be generated | |
393 | // | |
394 | // The algorithm is based on lagged Fibonacci sequences (first part) | |
395 | // combined with a congruential method (second part) | |
396 | // as described in the report by G.Marsaglia and A.Zaman | |
397 | // (FSU-SCRI-87-50 Florida State University 1987). | |
398 | // | |
399 | // Features : | |
400 | // 1) Period = 2**144 | |
401 | // 2) Same sequence of 24-bit real numbers on all common machines | |
402 | ||
403 | // Determine the minimum of a and b | |
404 | Float_t rmin=a; | |
405 | if (a > b) rmin=b; | |
406 | ||
407 | // First generate random numbers within <0,1> | |
408 | if (n > 0) // Check n value | |
409 | { | |
410 | for (Int_t jvec=0; jvec<n; jvec++) | |
411 | { | |
412 | // First part of the combination : F(97,33,*) | |
413 | Float_t unirnu=fU[fI-1]-fU[fJ-1]; | |
414 | if (unirnu < 0) unirnu+=1.; | |
415 | fU[fI-1]=unirnu; | |
416 | fI-=1; | |
417 | if (fI == 0) fI=97; | |
418 | fJ-=1; | |
419 | if (fJ == 0) fJ=97; | |
420 | ||
421 | // Second part of the combination | |
422 | fC-=fCd; | |
423 | if (fC < 0.) fC+=fCm; | |
424 | unirnu-=fC; | |
425 | if (unirnu < 0.) unirnu+=1.; | |
426 | ||
427 | // Update the counting system to enable sequence continuation | |
428 | // at an arbitrary starting position. | |
429 | // Two counters have been introduced to avoid overflow | |
430 | // fCnt1 : Counter which goes up to fClip | |
431 | // and is reset when fClip is reached | |
432 | // fCnt2 : Counts the number of times fClip has been reached | |
433 | fCnt1+=1; | |
434 | if (fCnt1 >= fClip) | |
435 | { | |
436 | fCnt1=0; | |
437 | fCnt2+=1; | |
438 | } | |
439 | ||
440 | if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range | |
441 | ||
442 | // Fill the vector within the selected range | |
443 | vec[jvec]=rmin+fabs(unirnu*(a-b)); | |
444 | } | |
445 | } | |
446 | else | |
447 | { | |
448 | cout << " *AliRandom::Uniform* Invalid value n = " << n << endl; | |
449 | } | |
450 | } | |
451 | /////////////////////////////////////////////////////////////////////////// | |
452 | void AliRandom::Uniform(Float_t* vec,Int_t n) | |
453 | { | |
454 | // Generate a vector of uniform random numbers in the interval <0,1> | |
455 | // This saves lots of (member)function calls in case many random | |
456 | // numbers are needed at once. | |
457 | // | |
458 | // n = The number of random numbers to be generated | |
459 | ||
460 | Uniform(vec,n,0.,1.); | |
461 | } | |
462 | /////////////////////////////////////////////////////////////////////////// | |
463 | void AliRandom::Uniform(Int_t n) | |
464 | { | |
465 | // Generate n uniform random numbers in in one go. | |
466 | // This saves lots of (member)function calls in case one needs to skip | |
467 | // to a certain point in a sequence. | |
468 | // | |
469 | // n = The number of random numbers to be generated | |
470 | // | |
471 | // Note : No check is made here to exclude 0 from the range. | |
472 | // It's only the number of generated randoms that counts | |
473 | // | |
474 | // The algorithm is based on lagged Fibonacci sequences (first part) | |
475 | // combined with a congruential method (second part) | |
476 | // as described in the report by G.Marsaglia and A.Zaman | |
477 | // (FSU-SCRI-87-50 Florida State University 1987). | |
478 | // | |
479 | // Features : | |
480 | // 1) Period = 2**144 | |
481 | // 2) Same sequence of 24-bit real numbers on all common machines | |
482 | ||
483 | if (n > 0) // Check n value | |
484 | { | |
485 | for (Int_t jvec=0; jvec<n; jvec++) | |
486 | { | |
487 | // First part of the combination : F(97,33,*) | |
488 | Float_t unirnu=fU[fI-1]-fU[fJ-1]; | |
489 | if (unirnu < 0) unirnu+=1.; | |
490 | fU[fI-1]=unirnu; | |
491 | fI-=1; | |
492 | if (fI == 0) fI=97; | |
493 | fJ-=1; | |
494 | if (fJ == 0) fJ=97; | |
495 | ||
496 | // Second part of the combination | |
497 | fC-=fCd; | |
498 | if (fC < 0.) fC+=fCm; | |
499 | unirnu-=fC; | |
500 | if (unirnu < 0.) unirnu+=1.; | |
501 | ||
502 | // Update the counting system to enable sequence continuation | |
503 | // at an arbitrary starting position. | |
504 | // Two counters have been introduced to avoid overflow | |
505 | // fCnt1 : Counter which goes up to fClip | |
506 | // and is reset when fClip is reached | |
507 | // fCnt2 : Counts the number of times fClip has been reached | |
508 | fCnt1+=1; | |
509 | if (fCnt1 >= fClip) | |
510 | { | |
511 | fCnt1=0; | |
512 | fCnt2+=1; | |
513 | } | |
514 | } | |
515 | } | |
516 | else | |
517 | { | |
518 | cout << " *AliRandom::Uniform* Invalid value n = " << n << endl; | |
519 | } | |
520 | } | |
521 | /////////////////////////////////////////////////////////////////////////// | |
522 | Float_t AliRandom::Gauss(Float_t mean,Float_t sigma) | |
523 | { | |
524 | // Generate gaussian distributed random numbers with certain mean and sigma | |
525 | // | |
526 | // Method : | |
527 | // P(x) = The gaussian distribution function | |
528 | // --> ln(P) provides an expression for (x-xmean)**2 from which | |
529 | // the following expression for x can be obtained | |
530 | // | |
531 | // x = xmean +/- sigma * sqrt(-2*ln(q)) | |
532 | // | |
533 | // in which q is an expression in terms of pi, sigma and p and lies within | |
534 | // the interval <0,1>. | |
535 | // | |
536 | // The gaussian distributed x values are obtained as follows : | |
537 | // | |
538 | // 1) Two uniform random numbers q1 and q2 in <0,1> are generated. | |
539 | // 2) q1 is in fact a uniform generated value for P which is substituted | |
540 | // directly in the formula above. | |
541 | // 3) The value of q2 determines whether we use the + or - sign. | |
542 | ||
543 | // Generate the two uniform random numbers q1 and q2 in <0,1> | |
544 | Float_t q1,q2; | |
545 | q1=Uniform(); | |
546 | q2=Uniform(); | |
547 | ||
548 | // Use q1 and q2 to get the gaussian distributed random number | |
549 | Float_t pi=acos(-1.); | |
550 | Float_t unirng=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1)); | |
551 | ||
552 | return unirng; | |
553 | } | |
554 | /////////////////////////////////////////////////////////////////////////// | |
555 | Float_t AliRandom::Gauss() | |
556 | { | |
557 | // Generate gaussian distributed random numbers with mean=0 and sigma=1 | |
558 | ||
559 | return Gauss(0.,1.); | |
560 | } | |
561 | /////////////////////////////////////////////////////////////////////////// | |
562 | void AliRandom::Gauss(Float_t* vec,Int_t n,Float_t mean,Float_t sigma) | |
563 | { | |
564 | // Generate a vector of gaussian random numbers with certain mean and sigma | |
565 | // This saves lots of (member)function calls in case many random | |
566 | // numbers are needed at once. | |
567 | // | |
568 | // n = The number of random numbers to be generated | |
569 | ||
570 | if (n > 0) // Check n value | |
571 | { | |
572 | // The vector to hold the q1 and q2 values. | |
573 | // Two subsequent q[] values are used for q1 and q2 | |
574 | // in order to obtain identical random numbers in the vector | |
575 | // as when generating n single ones. | |
576 | Int_t m=2*n; | |
577 | Float_t* q=new Float_t[m]; | |
578 | Uniform(q,m); | |
579 | ||
580 | // Fill the vector with gaussian random numbers | |
581 | Float_t pi=acos(-1.); | |
582 | Float_t q1,q2; | |
583 | for (Int_t jvec=0; jvec<n; jvec++) | |
584 | { | |
585 | q1=q[jvec*2]; // use two subsequent q[] values | |
586 | q2=q[(jvec*2)+1]; | |
587 | vec[jvec]=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1)); | |
588 | } | |
589 | delete [] q; | |
590 | } | |
591 | else | |
592 | { | |
593 | cout << " *AliRandom::Gauss* Invalid value n = " << n << endl; | |
594 | } | |
595 | } | |
596 | /////////////////////////////////////////////////////////////////////////// | |
597 | void AliRandom::Gauss(Float_t* vec,Int_t n) | |
598 | { | |
599 | // Generate a vector of gaussian random numbers with mean=0 and sigma=1 | |
600 | // This saves lots of (member)function calls in case many random | |
601 | // numbers are needed at once. | |
602 | // | |
603 | // n = The number of random numbers to be generated | |
604 | ||
605 | Gauss(vec,n,0.,1.); | |
606 | } | |
607 | /////////////////////////////////////////////////////////////////////////// | |
608 | Float_t AliRandom::Poisson(Float_t mean) | |
609 | { | |
610 | // Generate Poisson distributed random numbers with certain mean | |
611 | // | |
612 | // Method : | |
613 | // | |
614 | // P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function | |
615 | // | |
616 | // with : n = 0,1,2,3,... and mean > 0 | |
617 | // | |
618 | // To generate the distribution, the "sum trick" is used as mentioned | |
619 | // in "Formulae and Methods in Experimental data Evaluation Vol. 1" | |
620 | ||
621 | Float_t unirnp=0.; // Initialise the random number value | |
622 | ||
623 | if (mean <= 0.) // Check the mean value | |
624 | { | |
625 | cout << " *AliRandom::Poisson* Invalid value mean = " << mean | |
626 | << " Should be positive." << endl; | |
627 | } | |
628 | else | |
629 | { | |
630 | if (mean > 80.) // Use gaussian dist. for high mean values to save time | |
631 | { | |
632 | Float_t grndm=Gauss(); | |
633 | Float_t rpois=mean+grndm*sqrt(mean); | |
634 | Int_t npois=int(rpois); | |
635 | if ((rpois-float(npois)) > 0.5) npois++; | |
636 | unirnp=float(npois); | |
637 | } | |
638 | else // Construct a Poisson random number from a uniform one | |
639 | { | |
640 | Int_t npois=-1; | |
641 | Float_t expxm=exp(-mean); | |
642 | Float_t poitst=1.; | |
643 | while (poitst > expxm) | |
644 | { | |
645 | Float_t rndm=Uniform(); | |
646 | npois++; | |
647 | poitst=poitst*rndm; | |
648 | } | |
649 | unirnp=float(npois); | |
650 | } // end of check for using Gauss method | |
651 | } // end of mean validity checkn | |
652 | return unirnp; | |
653 | } | |
654 | /////////////////////////////////////////////////////////////////////////// | |
655 | void AliRandom::Poisson(Float_t* vec,Int_t n,Float_t mean) | |
656 | { | |
657 | // Generate a vector of Poisson dist. random numbers with certain mean | |
658 | // This saves lots of (member)function calls in case many random | |
659 | // numbers are needed at once. | |
660 | // | |
661 | // n = The number of random numbers to be generated | |
662 | // | |
663 | // Method : | |
664 | // | |
665 | // P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function | |
666 | // | |
667 | // with : n = 0,1,2,3,... and mean > 0 | |
668 | // | |
669 | // To generate the distribution, the "sum trick" is used as mentioned | |
670 | // in "Formulae and Methods in Experimental data Evaluation Vol. 1" | |
671 | ||
672 | if (n <= 0) // Check n value | |
673 | { | |
674 | cout << " *AliRandom::Poisson* Invalid value n = " << n << endl; | |
675 | } | |
676 | else | |
677 | { | |
678 | if (mean <= 0.) // Check the mean value | |
679 | { | |
680 | cout << " *AliRandom::Poisson* Invalid value mean = " << mean | |
681 | << " Should be positive." << endl; | |
682 | } | |
683 | else | |
684 | { | |
685 | if (mean > 80.) // Use gaussian dist. for high mean values to save time | |
686 | { | |
687 | Float_t* grndm=new Float_t[n]; | |
688 | Gauss(grndm,n); | |
689 | Int_t npois; | |
690 | Float_t rpois; | |
691 | for (Int_t jvec=0; jvec<n; jvec++) | |
692 | { | |
693 | rpois=mean+grndm[jvec]*sqrt(mean); | |
694 | npois=int(rpois); | |
695 | if ((rpois-float(npois)) > 0.5) npois++; | |
696 | vec[jvec]=float(npois); | |
697 | } | |
698 | delete [] grndm; | |
699 | } | |
700 | else // Construct Poisson random numbers from a uniform ones | |
701 | { | |
702 | Float_t expxm=exp(-mean); | |
703 | Int_t npois; | |
704 | Float_t poitst; | |
705 | for (Int_t jvec=0; jvec<n; jvec++) | |
706 | { | |
707 | npois=-1; | |
708 | poitst=1.; | |
709 | while (poitst > expxm) | |
710 | { | |
711 | Float_t rndm=Uniform(); | |
712 | npois++; | |
713 | poitst=poitst*rndm; | |
714 | } | |
715 | vec[jvec]=float(npois); | |
716 | } | |
717 | } // end of check for using Gauss method | |
718 | } // end of mean validity check | |
719 | } // end of n validity check | |
720 | } | |
721 | /////////////////////////////////////////////////////////////////////////// | |
722 | void AliRandom::SetUser(Float_t a,Float_t b,Int_t n,Float_t (*f)(Float_t)) | |
723 | { | |
724 | // Determine the area under f(x) as a function of x | |
725 | // This is called the "area function" and serves as a basis to | |
726 | // provide random numbers in [a,b] according to the user defined | |
727 | // distribution f(x). | |
728 | // The area function is normalised such that the most extreme | |
729 | // value is 1 or -1. | |
730 | ||
731 | fNa=n+1; // The number of bins for the area function | |
732 | fXa=new Float_t[fNa]; // The binned x values of the area function | |
733 | fYa=new Float_t[fNa]; // The corresponding summed f(x) values | |
734 | fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates | |
735 | ||
736 | Float_t xmin=a; | |
737 | if (a > b) xmin=b; | |
738 | Float_t step=fabs(a-b)/float(n); | |
739 | ||
740 | Float_t x; | |
741 | Float_t extreme=0; | |
742 | for (Int_t i=0; i<fNa; i++) // Fill bins of the area function | |
743 | { | |
744 | x=xmin+float(i)*step; | |
745 | fXa[i]=x; | |
746 | fYa[i]=f(x); | |
747 | if (i > 0) fYa[i]+=fYa[i-1]; | |
748 | if (fabs(fYa[i]) > extreme) extreme=fabs(fYa[i]); | |
749 | } | |
750 | fYamin=fYa[0]/extreme; | |
751 | fYamax=fYa[0]/extreme; | |
752 | for (Int_t j=0; j<fNa; j++) // Normalise the area function | |
753 | { | |
754 | fYa[j]=fYa[j]/extreme; | |
755 | if (fYa[j] < fYamin) fYamin=fYa[j]; | |
756 | if (fYa[j] > fYamax) fYamax=fYa[j]; | |
757 | } | |
758 | } | |
759 | /////////////////////////////////////////////////////////////////////////// | |
760 | void AliRandom::SetUser(Float_t* x,Float_t* y,Int_t n) | |
761 | { | |
762 | // Determine the area under y[i] as a function of x[i] | |
763 | // This is called the "area function" and serves as a basis to | |
764 | // provide random numbers in x according to the user provided | |
765 | // distribution (x[i],y[i]). | |
766 | // The area function is normalised such that the most extreme | |
767 | // value is 1 or -1. | |
768 | ||
769 | fNa=n; // The number of bins for the area function | |
770 | fXa=new Float_t[fNa]; // The binned x values of the area function | |
771 | fYa=new Float_t[fNa]; // The corresponding summed y values | |
772 | fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates | |
773 | ||
774 | // Order input data with increasing x | |
775 | fXa[0]=x[0]; | |
776 | fYa[0]=y[0]; | |
777 | for (Int_t i=1; i<fNa; i++) // Loop over x[] | |
778 | { | |
779 | for (Int_t j=0; j<i; j++) // Loop over xa[] | |
780 | { | |
781 | if (x[i] < fXa[j]) | |
782 | { | |
783 | for (Int_t k=i; k>=j; k--) // Create empty position | |
784 | { | |
785 | fXa[k+1]=fXa[k]; | |
786 | fYa[k+1]=fYa[k]; | |
787 | } | |
788 | fXa[j]=x[i]; // Put x[] value in empty position | |
789 | fYa[j]=y[i]; // Put y[] value in empty position | |
790 | break; // Go for next x[] value | |
791 | } | |
792 | if (j == i-1) // This x[] value is the largest so far | |
793 | { | |
794 | fXa[i]=x[i]; // Put x[] value at the end of x[] | |
795 | fYa[i]=y[i]; // Put y[] value at the end of y[] | |
796 | } | |
797 | } // End loop over fXa[] | |
798 | } // End loop over x[] | |
799 | ||
800 | Float_t extreme=0; | |
801 | for (Int_t l=0; l<fNa; l++) // Fill area function | |
802 | { | |
803 | if (l > 0) fYa[l]+=fYa[l-1]; | |
804 | if (fabs(fYa[l]) > extreme) extreme=fabs(fYa[l]); | |
805 | } | |
806 | fYamin=fYa[0]/extreme; | |
807 | fYamax=fYa[0]/extreme; | |
808 | for (Int_t m=0; m<fNa; m++) // Normalise the area function | |
809 | { | |
810 | fYa[m]=fYa[m]/extreme; | |
811 | if (fYa[m] < fYamin) fYamin=fYa[m]; | |
812 | if (fYa[m] > fYamax) fYamax=fYa[m]; | |
813 | } | |
814 | } | |
815 | /////////////////////////////////////////////////////////////////////////// | |
816 | Float_t AliRandom::User() | |
817 | { | |
818 | // Provide a random number according to the user defined distribution | |
819 | // | |
820 | // Method : | |
821 | // -------- | |
822 | // Select by a uniform random number a certain area fraction (from fYa[]) | |
823 | // of the area function. | |
824 | // The required random number is given by the corresponding x value (fXa[]) | |
825 | // of the area function. | |
826 | // In case of more than one x value candidate, select randomly one of them. | |
827 | ||
828 | Float_t unirnf=0; | |
829 | ||
830 | Float_t ra=Uniform(fYamin,fYamax); // Random value of the area function | |
831 | Float_t dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered | |
832 | Int_t ncand=0; | |
833 | Float_t dist; | |
834 | for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra | |
835 | { | |
836 | dist=fabs(ra-fYa[i]); | |
837 | if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate | |
838 | { | |
839 | ncand++; | |
840 | if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate | |
841 | dmin=dist; | |
842 | fIbins[ncand-1]=i+1; | |
843 | } | |
844 | } | |
845 | ||
846 | Int_t jbin=0; // The bin number to hold the required x value | |
847 | if (ncand == 1) jbin=fIbins[0]; | |
848 | ||
849 | if (ncand > 1) // Multiple x value candidates --> pick one randomly | |
850 | { | |
851 | Float_t cand=Uniform(1.,float(ncand)); | |
852 | Int_t jcand=int(cand); | |
853 | if ((cand-float(jcand)) > 0.5) jcand++; | |
854 | jbin=fIbins[jcand-1]; | |
855 | } | |
856 | ||
857 | if (jbin > 0) // Pick randomly a value in this x-bin | |
858 | { | |
859 | Float_t xlow=fXa[jbin-1]; | |
860 | if (jbin > 1) xlow=fXa[jbin-2]; | |
861 | Float_t xup=fXa[jbin-1]; | |
862 | if (jbin < fNa-1) xup=fXa[jbin]; | |
863 | unirnf=Uniform(xlow,xup); | |
864 | } | |
865 | ||
866 | if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl; | |
867 | ||
868 | return unirnf; | |
869 | } | |
870 | /////////////////////////////////////////////////////////////////////////// | |
871 | void AliRandom::User(Float_t* vec,Int_t n) | |
872 | { | |
873 | // Generate a vector of random numbers according to a user defined dist. | |
874 | // This saves lots of (member)function calls in case many random | |
875 | // numbers are needed at once. | |
876 | // | |
877 | // n = The number of random numbers to be generated | |
878 | // | |
879 | // Method : | |
880 | // -------- | |
881 | // Select by a uniform random number a certain area fraction (from fYa[]) | |
882 | // of the area function. | |
883 | // The required random number is given by the corresponding x value (fXa[]) | |
884 | // of the area function. | |
885 | // In case of more than one x value candidate, select randomly one of them. | |
886 | ||
887 | Float_t unirnf,ra,dmin,dist; | |
888 | Int_t ncand,jbin; | |
889 | for (Int_t jvec=0; jvec<n; jvec++) | |
890 | { | |
891 | unirnf=0; | |
892 | ra=Uniform(fYamin,fYamax); // Random value of the area function | |
893 | dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered | |
894 | ncand=0; | |
895 | for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra | |
896 | { | |
897 | dist=fabs(ra-fYa[i]); | |
898 | if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate | |
899 | { | |
900 | ncand++; | |
901 | if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate | |
902 | dmin=dist; | |
903 | fIbins[ncand-1]=i+1; | |
904 | } | |
905 | } | |
906 | ||
907 | jbin=0; // The bin number to hold the required x value | |
908 | if (ncand == 1) jbin=fIbins[0]; | |
909 | ||
910 | if (ncand > 1) // Multiple x value candidates --> pick one randomly | |
911 | { | |
912 | Float_t cand=Uniform(1.,float(ncand)); | |
913 | Int_t jcand=int(cand); | |
914 | if ((cand-float(jcand)) > 0.5) jcand++; | |
915 | jbin=fIbins[jcand-1]; | |
916 | } | |
917 | ||
918 | if (jbin > 0) // Pick randomly a value in this x-bin | |
919 | { | |
920 | Float_t xlow=fXa[jbin-1]; | |
921 | if (jbin > 1) xlow=fXa[jbin-2]; | |
922 | Float_t xup=fXa[jbin-1]; | |
923 | if (jbin < fNa-1) xup=fXa[jbin]; | |
924 | unirnf=Uniform(xlow,xup); | |
925 | } | |
926 | ||
927 | if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl; | |
928 | ||
929 | vec[jvec]=unirnf; | |
930 | ||
931 | } | |
932 | } | |
933 | /////////////////////////////////////////////////////////////////////////// |