RALICE/AliRandom.cxx

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