CommitLineData
4c039060 1/**************************************************************************
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
112ClassImp(AliRandom) // Class implementation to enable ROOT I/O
113
114AliRandom::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///////////////////////////////////////////////////////////////////////////
129AliRandom::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///////////////////////////////////////////////////////////////////////////
136AliRandom::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///////////////////////////////////////////////////////////////////////////
152AliRandom::~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///////////////////////////////////////////////////////////////////////////
163void 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///////////////////////////////////////////////////////////////////////////
247void 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 303Int_t AliRandom::GetSeed() const
d88f97cc 304{
305// Provide the current seed value
306 return fSeed;
307}
308///////////////////////////////////////////////////////////////////////////
261c0caf 309Int_t AliRandom::GetCnt1() const
d88f97cc 310{
311// Provide the current value of the counter cnt1
312 return fCnt1;
313}
314///////////////////////////////////////////////////////////////////////////
261c0caf 315Int_t AliRandom::GetCnt2() const
d88f97cc 316{
317// Provide the current value of the counter cnt2
318 return fCnt2;
319}
320///////////////////////////////////////////////////////////////////////////
261c0caf 321void 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///////////////////////////////////////////////////////////////////////////
328Float_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///////////////////////////////////////////////////////////////////////////
374Float_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///////////////////////////////////////////////////////////////////////////
386void 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///////////////////////////////////////////////////////////////////////////
452void 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///////////////////////////////////////////////////////////////////////////
463void 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///////////////////////////////////////////////////////////////////////////
522Float_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///////////////////////////////////////////////////////////////////////////
555Float_t AliRandom::Gauss()
556{
557// Generate gaussian distributed random numbers with mean=0 and sigma=1
558
559 return Gauss(0.,1.);
560}
561///////////////////////////////////////////////////////////////////////////
562void 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///////////////////////////////////////////////////////////////////////////
597void 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///////////////////////////////////////////////////////////////////////////
608Float_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///////////////////////////////////////////////////////////////////////////
655void 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///////////////////////////////////////////////////////////////////////////
722void 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///////////////////////////////////////////////////////////////////////////
760void 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///////////////////////////////////////////////////////////////////////////
816Float_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///////////////////////////////////////////////////////////////////////////
871void 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///////////////////////////////////////////////////////////////////////////