27-may-2001 NvE New class AliEvent introduced and RALICELinkDef.h & co. updated accor...
[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.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
106 //- Modified: NvE $Date$ UU-SAP Utrecht
107 ///////////////////////////////////////////////////////////////////////////
108
109 #include "AliRandom.h"
110  
111 ClassImp(AliRandom) // Class implementation to enable ROOT I/O
112  
113 AliRandom::AliRandom()
114 {
115 // Creation of an AliRandom object and default initialisation.
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 ///////////////////////////////////////////////////////////////////////////
128 AliRandom::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 ///////////////////////////////////////////////////////////////////////////
135 AliRandom::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 ///////////////////////////////////////////////////////////////////////////
151 AliRandom::~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 ///////////////////////////////////////////////////////////////////////////
162 void 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 ///////////////////////////////////////////////////////////////////////////
246 void 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 ///////////////////////////////////////////////////////////////////////////
302 Int_t AliRandom::GetSeed()
303 {
304 // Provide the current seed value
305  return fSeed;
306 }
307 ///////////////////////////////////////////////////////////////////////////
308 Int_t AliRandom::GetCnt1()
309 {
310 // Provide the current value of the counter cnt1
311  return fCnt1;
312 }
313 ///////////////////////////////////////////////////////////////////////////
314 Int_t AliRandom::GetCnt2()
315 {
316 // Provide the current value of the counter cnt2
317  return fCnt2;
318 }
319 ///////////////////////////////////////////////////////////////////////////
320 void 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 ///////////////////////////////////////////////////////////////////////////
327 Float_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 ///////////////////////////////////////////////////////////////////////////
373 Float_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 ///////////////////////////////////////////////////////////////////////////
385 void 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 ///////////////////////////////////////////////////////////////////////////
451 void 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 ///////////////////////////////////////////////////////////////////////////
462 void 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 ///////////////////////////////////////////////////////////////////////////
521 Float_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 ///////////////////////////////////////////////////////////////////////////
554 Float_t AliRandom::Gauss()
555 {
556 // Generate gaussian distributed random numbers with mean=0 and sigma=1
557  
558  return Gauss(0.,1.);
559 }
560 ///////////////////////////////////////////////////////////////////////////
561 void 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 ///////////////////////////////////////////////////////////////////////////
596 void 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 ///////////////////////////////////////////////////////////////////////////
607 Float_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 ///////////////////////////////////////////////////////////////////////////
654 void 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 ///////////////////////////////////////////////////////////////////////////
721 void 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 ///////////////////////////////////////////////////////////////////////////
759 void 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 ///////////////////////////////////////////////////////////////////////////
815 Float_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 ///////////////////////////////////////////////////////////////////////////
870 void 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 ///////////////////////////////////////////////////////////////////////////