3ccce9bdf60fb4be65ff928f1d5ddf1c9ba676b3
[u/mrichter/AliRoot.git] / RALICE / AliRandom.cxx
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$
18 Revision 1.2  1999/09/29 09:24:28  fca
19 Introduction of the Copyright and cvs Log
20
21 */
22
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
113 #include "AliRandom.h"
114  
115 ClassImp(AliRandom) // Class implementation to enable ROOT I/O
116  
117 AliRandom::AliRandom()
118 {
119 // Creation of an AliRandom object and default initialisation.
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 ///////////////////////////////////////////////////////////////////////////
132 AliRandom::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 ///////////////////////////////////////////////////////////////////////////
139 AliRandom::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 ///////////////////////////////////////////////////////////////////////////
155 AliRandom::~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 ///////////////////////////////////////////////////////////////////////////
166 void 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 ///////////////////////////////////////////////////////////////////////////
250 void 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 ///////////////////////////////////////////////////////////////////////////
306 Int_t AliRandom::GetSeed()
307 {
308 // Provide the current seed value
309  return fSeed;
310 }
311 ///////////////////////////////////////////////////////////////////////////
312 Int_t AliRandom::GetCnt1()
313 {
314 // Provide the current value of the counter cnt1
315  return fCnt1;
316 }
317 ///////////////////////////////////////////////////////////////////////////
318 Int_t AliRandom::GetCnt2()
319 {
320 // Provide the current value of the counter cnt2
321  return fCnt2;
322 }
323 ///////////////////////////////////////////////////////////////////////////
324 void 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 ///////////////////////////////////////////////////////////////////////////
331 Float_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 ///////////////////////////////////////////////////////////////////////////
377 Float_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 ///////////////////////////////////////////////////////////////////////////
389 void 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 ///////////////////////////////////////////////////////////////////////////
455 void 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 ///////////////////////////////////////////////////////////////////////////
466 void 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 ///////////////////////////////////////////////////////////////////////////
525 Float_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 ///////////////////////////////////////////////////////////////////////////
558 Float_t AliRandom::Gauss()
559 {
560 // Generate gaussian distributed random numbers with mean=0 and sigma=1
561  
562  return Gauss(0.,1.);
563 }
564 ///////////////////////////////////////////////////////////////////////////
565 void 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 ///////////////////////////////////////////////////////////////////////////
600 void 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 ///////////////////////////////////////////////////////////////////////////
611 Float_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 ///////////////////////////////////////////////////////////////////////////
658 void 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 ///////////////////////////////////////////////////////////////////////////
725 void 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 ///////////////////////////////////////////////////////////////////////////
763 void 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 ///////////////////////////////////////////////////////////////////////////
819 Float_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 ///////////////////////////////////////////////////////////////////////////
874 void 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 ///////////////////////////////////////////////////////////////////////////