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