RALICE/AliRandom.cxx

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