def METHOD_NAME(data, unweighted_mc, weights):
    """
    Barlow-Beeston log-likelihood (constant terms not omitted)
    Link to paper: https://doi.org/10.1016/0010-4655(93)90005-W
    -- Input variables --
    data = data histogram
    mc = weighted MC histogram
    unweighted_mc = unweighted MC histogream
    weights = weight of each bin
    -- Output --
    llh = LLH values in each bin
    -- Notes --
    Shape of data, mc, unweighted_mc, weights and llh must be identical
    """
    # The actual barlow LLH
    def llh(A_, k, w, a):
        SMALL_VAL = 1.e-10
        f = w*A_
        # Takes care of log(0) problems
        if not len(A_) > 1:
            f = max((f, SMALL_VAL))
            A_ = max((A_, SMALL_VAL))
        # The loggamma() terms takes care of the log(value!) for non-integer values
        return -1.*(k*np.log(f) - f + a*np.log(A_) - A_ - special.loggamma(k+1) - special.loggamma(a+1))
    A = np.array(unweighted_mc) # Expected unweighted counts in a bin
    # For each bin the appropriate 'A' will now be found using the current counts as a seed
    # This will be done by using a minimiser to ensure that the value of 'A' found minimises the -LLH
    # Find values of A such that llh in each bin is a maximum
    for i, val in enumerate(A):
        # If the unweighted MC counts in the bin is 0, A = 0
        if val == 0:
            continue
        # Otherwise, find the value of A
        arg = (data[i], weights[i], unweighted_mc[i])
        # Powell works fast and is fine for our purposes
        result = optimize.minimize(fun=llh, x0=val, args=arg, method='Powell')
        # Check that the minimisation ran properly
        if result.success:
            A[i] = result.x
        else:
            print("Something went wrong...")
            print("Minimiser message: ")
            print("------------------")
            print(result.message)
            return -np.inf
    LLH = llh(A, data, weights, unweighted_mc)
    return -1*LLH # Return LLH (not negative LLH)