def prob_imf(m1, m2, s1z, s2z, **kwargs):
    ''' Return probability density for power-law
        Parameters
        ----------
        m1: array
            Component masses 1
        m2: array
            Component masses 2
        s1z: array
            Aligned spin 1(Not in use currently)
        s2z:
            Aligned spin 2(Not in use currently)
        **kwargs: string
            Keyword arguments as model parameters

        Returns
        -------
        p_m1_m2: array
           the probability density for m1, m2 pair
    '''

    min_mass = kwargs.get('min_mass', 5.)
    max_mass = kwargs.get('max_mass', 95.)
    alpha = kwargs.get('alpha', -2.35)
    max_mtotal = min_mass + max_mass
    m1, m2 = np.array(m1), np.array(m2)

    C_imf = max_mass**(alpha + 1)/(alpha + 1)
    C_imf -= min_mass**(alpha + 1)/(alpha + 1)

    xx = np.minimum(m1, m2)
    m1 = np.maximum(m1, m2)
    m2 = xx

    bound = np.sign(max_mtotal - m1 - m2)
    bound += np.sign(max_mass - m1) * np.sign(m2 - min_mass)
    idx = np.where(bound != 2)

    p_m1_m2 = np.zeros_like(m1)
    idx = np.where(m1 <= max_mtotal/2.)
    p_m1_m2[idx] = (1./C_imf) * m1[idx]**alpha /(m1[idx] - min_mass)
    idx = np.where(m1 > max_mtotal/2.)
    p_m1_m2[idx] = (1./C_imf) * m1[idx]**alpha /(max_mass - m1[idx])
    p_m1_m2[idx] = 0

    return p_m1_m2/2.
Return probability density for power-law
        Parameters
        ----------
        m1: array
            Component masses 1
        m2: array
            Component masses 2
        s1z: array
            Aligned spin 1(Not in use currently)
        s2z:
            Aligned spin 2(Not in use currently)
        **kwargs: string
            Keyword arguments as model parameters

        Returns
        -------
        p_m1_m2: array
           the probability density for m1, m2 pair