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