def METHOD_NAME(data, sum_w, sum_w2, a=1, b=0):
    """
    Log-likelihood based on the poisson-gamma mixture. This is a Poisson likelihood using a Gamma prior.
    This implementation is based on the implementation of Austin Schneider (aschneider@icecube.wisc.edu)
    -- Input variables --
    data = data histogram
    sum_w = MC histogram
    sum_w2 = Uncertainty map (sum of weights squared in each bin)
    a, b = hyperparameters of gamma prior for MC counts; default values of a = 1 and b = 0 corresponds to LEff (eq 3.16) https://doi.org/10.1007/JHEP06(2019)030
           a = 0 and b = 0 corresponds to LMean (Table 2) https://doi.org/10.1007/JHEP06(2019)030
    -- Output --
    llh = LLH values in each bin
    -- Notes --
    Shape of data, sum_w, sum_w2 and llh are identical
    """
    llh = np.ones(data.shape) * -np.inf # Binwise LLH values
    bad_bins = np.logical_or(sum_w <= 0, sum_w2 < 0) # Bins where the MC is 0 or less than 0
    # LLH would be 0 for the bad bins if the data is 0
    zero_llh = np.logical_and(data == 0, bad_bins)
    llh[zero_llh] = 0 # Zero LLH for these bins if data is also 0
    good_bins = ~bad_bins
    poisson_bins = np.logical_and(sum_w2 == 0, good_bins) # In the limit that sum_w2 == 0, the llh converges to poisson
    llh[poisson_bins] = poissonLLH(data[poisson_bins], sum_w[poisson_bins]) # Poisson LLH since limiting case
   # Calculate hyperparameters for the gamma posterior for MC counts
    regular_bins = np.logical_and(good_bins, ~poisson_bins) # Bins on which the poisson_gamma LLH would be evaluated
    alpha = sum_w[regular_bins]**2./sum_w2[regular_bins] + a
    beta = sum_w[regular_bins]/sum_w2[regular_bins] + b
    k = data[regular_bins]
    # Poisson-gamma LLH
    L = alpha*np.log(beta) + special.loggamma(k+alpha).real - special.loggamma(k+1.0).real - (k+alpha)*np.log1p(beta) - special.loggamma(alpha).real
    llh[regular_bins] = L
    return llh