Hi,
I have found what I think to be a bug in the random number generation of RooGamma. The routine does not work for small gamma (called alpha elsewhere and from now on in this post to avoid confusion with the name of the distribution) values. The original paper describing the routine suggests that for values of alpha less than one, a random number can be sampled by first sampling from a gamma distribution of alpha+1 and multiplying by a uniform random number raised to the power of 1/alpha.
I have applied what I think is a reasonable fix:
void RooGamma::generateEvent(Int_t code)
{
assert(code==1) ;
//algorithm adapted from code example in:
//Marsaglia, G. and Tsang, W. W.
//A Simple Method for Generating GammPdf Variables
//ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
//
//The speed of this algorithm depends on the speed of generating normal variates.
//The algorithm is limited to $\gamma \geq 0$ !
while(1) {
double d = 0;
double c = 0;
double xgen = 0;
double v = 0;
double u = 0;
d = gamma -1./3.;
//Edit from Luke Winstrom
if(gamma<1) d++;
c = 1./TMath::Sqrt(9.*d);
while(v <= 0.){
xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + cxgen;
}
v = vvv; u = RooRandom::randomGenerator()->Uniform();
if( u < 1.-.0331(xgenxgen)(xgenxgen) ) {
if ( (((dv)* beta + mu ) < x.max()) && (((dv) beta + mu) > x.min()) ) {
x = ((dv) beta + mu) ;
break;
}
}
if( TMath::Log(u) < 0.5xgenxgen + d*(1.-v + TMath::Log(v)) ) {
if ( (((dv) beta + mu ) < x.max()) && (((dv) beta + mu) > x.min()) ) {
x = ((dv) beta + mu) ;
break;
}
}
}
//Edit from Luke Winstrom
if(gamma<1) x=x*pow(RooRandom::randomGenerator()->Uniform(),1/gamma);
return;
}
Cheers,
Luke