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