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 ///////////////////////////////////////////////////////////////////////////