RALICE/AliRandom.cxx

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