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