hggamma <- function(beta, sigma, lambda, timet)
{

s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)

b <- lambda/sigma
haz <- (lambda>0)*((b*(timet^(b-1))*dgamma(timet^b,shape=s,rate=r))/(1-pgamma(timet^b,shape=s,rate=r)))
haz <- haz+(lambda<0)*((-b*(timet^(b-1))*dgamma(timet^b,shape=s,rate=r))/pgamma(timet^b,shape=s,rate=r))
haz
}