1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // Generate universal random numbers on all common machines.
21 // Available distributions : Uniform, Gaussian, Poisson and
22 // User defined function
27 // 2) Same sequence of 24-bit real numbers on all common machines
31 // G.Marsaglia and A.Zaman, FSU-SCRI-87-50, Florida State University, 1987.
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
40 // AliRandom r; // Create a Random object with default sequence
42 // rndm=r.Uniform(); // Provide a uniform random number in <0,1>
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
49 // rndm=r.Gauss(); // Provide a Gaussian random number with
50 // // mean=0 and sigma=1
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
59 // rndm=r.Poisson(mean); // Provide a Poisson random number with
61 // r.Poisson(rvec,nmean); // n Poisson randoms with specified mean
64 // AliRandom p(seed); // Create a Random object with specified seed.
65 // // The sequence is started from scratch.
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.
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.
77 // Float_t udist(Float_t x) // A user defined distribution
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
87 // Float_t* x=new Float_t[nbins];
88 // Float_t* y=new Float_t[nbins];
90 // ... code to fill x[] and y[] ..
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
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.
105 //--- Author: Nick van Eijndhoven 11-oct-1997 UU-SAP Utrecht
106 //- Modified: NvE $Date$ UU-SAP Utrecht
107 ///////////////////////////////////////////////////////////////////////////
109 #include "AliRandom.h"
110 #include "Riostream.h"
112 ClassImp(AliRandom) // Class implementation to enable ROOT I/O
114 AliRandom::AliRandom()
116 // Creation of an AliRandom object and default initialisation.
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").
122 // Suggested test values : i=12 j=34 k=56 l=78 (see article)
123 // which corresponds to : seed = 53310452
125 Int_t seed=53310452; // Default seed
126 Start(seed,0,0); // Start the sequence for this seed from scratch
128 ///////////////////////////////////////////////////////////////////////////
129 AliRandom::AliRandom(Int_t seed)
131 // Creation of an AliRandom object and user defined initialisation
133 Start(seed,0,0); // Start the sequence for this seed from scratch
135 ///////////////////////////////////////////////////////////////////////////
136 AliRandom::AliRandom(Int_t seed,Int_t cnt1,Int_t cnt2)
138 // Creation of an AliRandom object and user defined initialisation
140 // seed is the seed to create the initial u[97] table.
141 // The range of the seed is : 0 <= seed <= 921350143
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
149 Start(seed,cnt1,cnt2); // Start the sequence from a user defined point
151 ///////////////////////////////////////////////////////////////////////////
152 AliRandom::~AliRandom()
154 // Destructor to delete memory allocated for the area function arrays
155 if (fXa) delete [] fXa;
157 if (fYa) delete [] fYa;
159 if (fIbins) delete [] fIbins;
162 ///////////////////////////////////////////////////////////////////////////
163 void AliRandom::Start(Int_t seed,Int_t cnt1,Int_t cnt2)
165 // Start a certain sequence from scratch or from a user defined point
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).
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").
175 // The range of the seed is : 0 <= seed <= 921350143
177 // Suggested test values : i=12 j=34 k=56 l=78 (see article)
178 // which corresponds to : seed = 53310452
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
186 // Reset the area function
192 // Clipping parameter to prevent overflow of the counting system
195 // Set the lags for the Fibonacci sequence of the first part
196 // The sequence is set to F(97,33,*) (see article)
200 // Unpack the seed value and determine i, j, k and l
203 Unpack(seed,i,j,k,l);
209 // Fill the starting table for the first part of the combination
212 for (Int_t ii=0; ii<97; ii++)
217 for (Int_t jj=1; jj<=24; jj++)
219 m=(((i*j)%179)*k)%179;
224 if ((l*m)%64 >= 32) s+=t;
230 // Initialise the second part of the combination
231 fC=362436./16777216.;
232 fCd=7654321./16777216.;
233 fCm=16777213./16777216.;
235 // Generate random numbers upto the user selected starting point
236 // on basis of the counting system
237 if (cnt1 > 0) Uniform(cnt1);
240 for (Int_t n=1; n<=cnt2; n++)
246 ///////////////////////////////////////////////////////////////////////////
247 void AliRandom::Unpack(Int_t seed,Int_t& i,Int_t& j,Int_t& k,Int_t& l)
249 // Unpack the seed into the four startup parameters i,j,k and l
251 // The range of the seed is : 0 <= seed <= 921350143
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
257 // Taking into account the period of the 3-lagged Fibonacci sequence
258 // The following "bad" combinations have to be ruled out :
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
275 // and use the formula :
276 // seed = (i-2)*176*176*169 + (j-2)*176*169 + (k-2)*169 + l
278 if ((seed >= 0) && (seed <= 921350143)) // Check seed value
281 Int_t imin2=idum/(176*176*169);
282 idum=idum%(176*176*169);
283 Int_t jmin2=idum/(176*169);
285 Int_t kmin2=idum/169;
294 cout << " *AliRandom::unpack()* Unallowed seed value encountered."
295 << " seed = " << seed << endl;
296 cout << " Seed will be set to default value." << endl;
298 seed=53310452; // Default seed
299 Start(seed,0,0); // Start the sequence for this seed from scratch
302 ///////////////////////////////////////////////////////////////////////////
303 Int_t AliRandom::GetSeed() const
305 // Provide the current seed value
308 ///////////////////////////////////////////////////////////////////////////
309 Int_t AliRandom::GetCnt1() const
311 // Provide the current value of the counter cnt1
314 ///////////////////////////////////////////////////////////////////////////
315 Int_t AliRandom::GetCnt2() const
317 // Provide the current value of the counter cnt2
320 ///////////////////////////////////////////////////////////////////////////
321 void AliRandom::Data() const
323 // Print the current seed, cnt1 and cnt2 values
324 cout << " *Random* seed = " << fSeed
325 << " cnt1 = " << fCnt1 << " cnt2 = " << fCnt2 << endl;
327 ///////////////////////////////////////////////////////////////////////////
328 Float_t AliRandom::Uniform()
330 // Generate uniform random numbers in the interval <0,1>
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).
338 // 1) Period = 2**144
339 // 2) Same sequence of 24-bit real numbers on all common machines
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.;
350 // Second part of the combination (see article for explanation)
352 if (fC < 0.) fC+=fCm;
354 if (unirnu < 0.) unirnu+=1.;
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
369 if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
373 ///////////////////////////////////////////////////////////////////////////
374 Float_t AliRandom::Uniform(Float_t a,Float_t b)
376 // Generate uniform random numbers in the interval <a,b>
380 Float_t rndm=Uniform();
381 rndm=rmin+fabs(rndm*(a-b));
385 ///////////////////////////////////////////////////////////////////////////
386 void AliRandom::Uniform(Float_t* vec,Int_t n,Float_t a,Float_t b)
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.
392 // n = The number of random numbers to be generated
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).
400 // 1) Period = 2**144
401 // 2) Same sequence of 24-bit real numbers on all common machines
403 // Determine the minimum of a and b
407 // First generate random numbers within <0,1>
408 if (n > 0) // Check n value
410 for (Int_t jvec=0; jvec<n; jvec++)
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.;
421 // Second part of the combination
423 if (fC < 0.) fC+=fCm;
425 if (unirnu < 0.) unirnu+=1.;
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
440 if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
442 // Fill the vector within the selected range
443 vec[jvec]=rmin+fabs(unirnu*(a-b));
448 cout << " *AliRandom::Uniform* Invalid value n = " << n << endl;
451 ///////////////////////////////////////////////////////////////////////////
452 void AliRandom::Uniform(Float_t* vec,Int_t n)
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.
458 // n = The number of random numbers to be generated
460 Uniform(vec,n,0.,1.);
462 ///////////////////////////////////////////////////////////////////////////
463 void AliRandom::Uniform(Int_t n)
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.
469 // n = The number of random numbers to be generated
471 // Note : No check is made here to exclude 0 from the range.
472 // It's only the number of generated randoms that counts
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).
480 // 1) Period = 2**144
481 // 2) Same sequence of 24-bit real numbers on all common machines
483 if (n > 0) // Check n value
485 for (Int_t jvec=0; jvec<n; jvec++)
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.;
496 // Second part of the combination
498 if (fC < 0.) fC+=fCm;
500 if (unirnu < 0.) unirnu+=1.;
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
518 cout << " *AliRandom::Uniform* Invalid value n = " << n << endl;
521 ///////////////////////////////////////////////////////////////////////////
522 Float_t AliRandom::Gauss(Float_t mean,Float_t sigma)
524 // Generate gaussian distributed random numbers with certain mean and sigma
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
531 // x = xmean +/- sigma * sqrt(-2*ln(q))
533 // in which q is an expression in terms of pi, sigma and p and lies within
534 // the interval <0,1>.
536 // The gaussian distributed x values are obtained as follows :
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.
543 // Generate the two uniform random numbers q1 and q2 in <0,1>
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));
554 ///////////////////////////////////////////////////////////////////////////
555 Float_t AliRandom::Gauss()
557 // Generate gaussian distributed random numbers with mean=0 and sigma=1
561 ///////////////////////////////////////////////////////////////////////////
562 void AliRandom::Gauss(Float_t* vec,Int_t n,Float_t mean,Float_t sigma)
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.
568 // n = The number of random numbers to be generated
570 if (n > 0) // Check n value
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.
577 Float_t* q=new Float_t[m];
580 // Fill the vector with gaussian random numbers
581 Float_t pi=acos(-1.);
583 for (Int_t jvec=0; jvec<n; jvec++)
585 q1=q[jvec*2]; // use two subsequent q[] values
587 vec[jvec]=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1));
593 cout << " *AliRandom::Gauss* Invalid value n = " << n << endl;
596 ///////////////////////////////////////////////////////////////////////////
597 void AliRandom::Gauss(Float_t* vec,Int_t n)
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.
603 // n = The number of random numbers to be generated
607 ///////////////////////////////////////////////////////////////////////////
608 Float_t AliRandom::Poisson(Float_t mean)
610 // Generate Poisson distributed random numbers with certain mean
614 // P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function
616 // with : n = 0,1,2,3,... and mean > 0
618 // To generate the distribution, the "sum trick" is used as mentioned
619 // in "Formulae and Methods in Experimental data Evaluation Vol. 1"
621 Float_t unirnp=0.; // Initialise the random number value
623 if (mean <= 0.) // Check the mean value
625 cout << " *AliRandom::Poisson* Invalid value mean = " << mean
626 << " Should be positive." << endl;
630 if (mean > 80.) // Use gaussian dist. for high mean values to save time
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++;
638 else // Construct a Poisson random number from a uniform one
641 Float_t expxm=exp(-mean);
643 while (poitst > expxm)
645 Float_t rndm=Uniform();
650 } // end of check for using Gauss method
651 } // end of mean validity checkn
654 ///////////////////////////////////////////////////////////////////////////
655 void AliRandom::Poisson(Float_t* vec,Int_t n,Float_t mean)
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.
661 // n = The number of random numbers to be generated
665 // P(n) = exp(-mean)*mean**n/n! is the Poisson distribution function
667 // with : n = 0,1,2,3,... and mean > 0
669 // To generate the distribution, the "sum trick" is used as mentioned
670 // in "Formulae and Methods in Experimental data Evaluation Vol. 1"
672 if (n <= 0) // Check n value
674 cout << " *AliRandom::Poisson* Invalid value n = " << n << endl;
678 if (mean <= 0.) // Check the mean value
680 cout << " *AliRandom::Poisson* Invalid value mean = " << mean
681 << " Should be positive." << endl;
685 if (mean > 80.) // Use gaussian dist. for high mean values to save time
687 Float_t* grndm=new Float_t[n];
691 for (Int_t jvec=0; jvec<n; jvec++)
693 rpois=mean+grndm[jvec]*sqrt(mean);
695 if ((rpois-float(npois)) > 0.5) npois++;
696 vec[jvec]=float(npois);
700 else // Construct Poisson random numbers from a uniform ones
702 Float_t expxm=exp(-mean);
705 for (Int_t jvec=0; jvec<n; jvec++)
709 while (poitst > expxm)
711 Float_t rndm=Uniform();
715 vec[jvec]=float(npois);
717 } // end of check for using Gauss method
718 } // end of mean validity check
719 } // end of n validity check
721 ///////////////////////////////////////////////////////////////////////////
722 void AliRandom::SetUser(Float_t a,Float_t b,Int_t n,Float_t (*f)(Float_t))
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
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
738 Float_t step=fabs(a-b)/float(n);
742 for (Int_t i=0; i<fNa; i++) // Fill bins of the area function
744 x=xmin+float(i)*step;
747 if (i > 0) fYa[i]+=fYa[i-1];
748 if (fabs(fYa[i]) > extreme) extreme=fabs(fYa[i]);
750 fYamin=fYa[0]/extreme;
751 fYamax=fYa[0]/extreme;
752 for (Int_t j=0; j<fNa; j++) // Normalise the area function
754 fYa[j]=fYa[j]/extreme;
755 if (fYa[j] < fYamin) fYamin=fYa[j];
756 if (fYa[j] > fYamax) fYamax=fYa[j];
759 ///////////////////////////////////////////////////////////////////////////
760 void AliRandom::SetUser(Float_t* x,Float_t* y,Int_t n)
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
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
774 // Order input data with increasing x
777 for (Int_t i=1; i<fNa; i++) // Loop over x[]
779 for (Int_t j=0; j<i; j++) // Loop over xa[]
783 for (Int_t k=i; k>=j; k--) // Create empty position
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
792 if (j == i-1) // This x[] value is the largest so far
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[]
797 } // End loop over fXa[]
798 } // End loop over x[]
801 for (Int_t l=0; l<fNa; l++) // Fill area function
803 if (l > 0) fYa[l]+=fYa[l-1];
804 if (fabs(fYa[l]) > extreme) extreme=fabs(fYa[l]);
806 fYamin=fYa[0]/extreme;
807 fYamax=fYa[0]/extreme;
808 for (Int_t m=0; m<fNa; m++) // Normalise the area function
810 fYa[m]=fYa[m]/extreme;
811 if (fYa[m] < fYamin) fYamin=fYa[m];
812 if (fYa[m] > fYamax) fYamax=fYa[m];
815 ///////////////////////////////////////////////////////////////////////////
816 Float_t AliRandom::User()
818 // Provide a random number according to the user defined distribution
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.
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
834 for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
836 dist=fabs(ra-fYa[i]);
837 if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
840 if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
846 Int_t jbin=0; // The bin number to hold the required x value
847 if (ncand == 1) jbin=fIbins[0];
849 if (ncand > 1) // Multiple x value candidates --> pick one randomly
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];
857 if (jbin > 0) // Pick randomly a value in this x-bin
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);
866 if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl;
870 ///////////////////////////////////////////////////////////////////////////
871 void AliRandom::User(Float_t* vec,Int_t n)
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.
877 // n = The number of random numbers to be generated
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.
887 Float_t unirnf,ra,dmin,dist;
889 for (Int_t jvec=0; jvec<n; jvec++)
892 ra=Uniform(fYamin,fYamax); // Random value of the area function
893 dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered
895 for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
897 dist=fabs(ra-fYa[i]);
898 if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
901 if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
907 jbin=0; // The bin number to hold the required x value
908 if (ncand == 1) jbin=fIbins[0];
910 if (ncand > 1) // Multiple x value candidates --> pick one randomly
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];
918 if (jbin > 0) // Pick randomly a value in this x-bin
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);
927 if (ncand == 0) cout << " *AliRandom::User* No candidate found." << endl;
933 ///////////////////////////////////////////////////////////////////////////