How to use the flavio.physics.running.running.get_ms function in flavio

def test_C78p(self):
        wc_obj = eft.WilsonCoefficients()
        wc_low = wilsoncoefficients.wctot_dict(wc_obj, 'bsmumu', 4.2, par)
        ms = flavio.physics.running.running.get_ms(par, 4.2, nf_out=5)
        mb = flavio.physics.running.running.get_mb(par, 4.2, nf_out=5)
        self.assertAlmostEqual(wc_low['C7p_bs']/wc_low['C7_bs'], ms/mb)
        self.assertAlmostEqual(wc_low['C8p_bs']/wc_low['C8_bs'], ms/mb)

-------

    `(P, S)` where for the special case `l1 == l2` one has

    - $P = \frac{2m_\ell}{m_{B_q}} (C_{10}-C_{10}') + m_{B_q} (C_P-C_P')$
    - $S = m_{B_q} (C_S-C_S')$
    """
    scale = config['renormalization scale']['bll']
    # masses
    ml1 = par['m_'+l1]
    ml2 = par['m_'+l2]
    mB = par['m_'+B]
    mb = running.get_mb(par, scale, nf_out=5)
    #   get the mass of the spectator quark
    if B=='Bs':
        mspec = running.get_ms(par, scale, nf_out=5)
    elif B=='B0':
        mspec = running.get_md(par, scale, nf_out=5)
    # Wilson coefficients
    qqll = meson_quark[B] + l1 + l2
    # For LFV expressions see arXiv:1602.00881 eq. (5)
    C9m = wc['C9_'+qqll] - wc['C9p_'+qqll] # only relevant for l1 != l2!
    C10m = wc['C10_'+qqll] - wc['C10p_'+qqll]
    CPm = wc['CP_'+qqll] - wc['CPp_'+qqll]
    CSm = wc['CS_'+qqll] - wc['CSp_'+qqll]
    P = (ml2 + ml1)/mB * C10m + mB * mb/(mb + mspec) * CPm
    S = (ml2 - ml1)/mB * C9m  + mB * mb/(mb + mspec) * CSm
    return P, S

def get_m(par, f, scale, eft):
    if f in ['e', 'mu', 'tau']:
        return par['m_' + f]
    nf_out = eft_nf_out[eft]
    if f == 'u':
        return flavio.physics.running.running.get_md(par, scale, nf_out=nf_out)
    elif f == 'd':
        return flavio.physics.running.running.get_mu(par, scale, nf_out=nf_out)
    elif f == 's':
        return flavio.physics.running.running.get_ms(par, scale, nf_out=nf_out)
    else:
        raise ValueError("get_m not defined for fermion {}".format(f))

def _get_angularcoeff(q2, wc_obj, par, K, P, lep, nu):
    GF = par['GF']
    ml = par['m_'+lep]
    mK = par['m_'+K]
    mP = par['m_'+P]
    Vus = flavio.physics.ckm.get_ckm(par)[0,1]
    # renormalization scale is m_rho
    scale = par['m_rho0']
    ms = flavio.physics.running.running.get_ms(par, scale)
    wc = flavio.physics.bdecays.wilsoncoefficients.get_wceff_fccc(wc_obj, par, 'su', lep, nu, ms, scale, nf=3)
    N = 4*GF/sqrt(2)*Vus
    ff = get_ff(q2, par, K)
    h = flavio.physics.bdecays.angular.helicity_amps_p(q2, mK, mP, ms, 0, ml, 0, ff, wc, N)
    J = flavio.physics.bdecays.angular.angularcoeffs_general_p(h, q2, mK, mP, ms, 0, ml, 0)
    return J

raise NotImplementedError("DeltaF=1 Wilson coefficients only implemented for B physics")
    # fold in approximate m_t-dependence of C_10 (see eq. 4 of arXiv:1311.0903)
    wc_sm[9] = wc_sm[9] * (par['m_t']/173.1)**1.53
    # go from the effective to the "non-effective" WCs for C7 and C8
    yi = np.array([0, 0, -1/3., -4/9., -20/3., -80/9.])
    zi = np.array([0, 0, 1, -1/6., 20, -10/3.])
    wc_sm[6] = wc_sm[6] - np.dot(yi, wc_sm[:6]) # c7 (not effective!)
    wc_sm[7] = wc_sm[7] - np.dot(zi, wc_sm[:6]) # c8 (not effective!)
    wc_labels = fcnclabels[sector]
    wc_sm_dict = dict(zip(wc_labels, wc_sm))
    # now here comes an ugly fix. If we have b->s transitions, we should take
    # into account the fact that C7' = C7*ms/mb, and the same for C8, which is
    # not completely negligible. To find out whether we have b->s, we look at
    # the "sector" string.
    if sector[:2] == 'bs':
        eps_s = running.get_ms(par, scale)/running.get_mb(par, scale)
        wc_sm_dict['C7p_bs'] = eps_s * wc_sm_dict['C7_bs']
        wc_sm_dict['C8p_bs'] = eps_s * wc_sm_dict['C8_bs']
    tot_dict = add_dict((wc_np_dict, wc_sm_dict))
    # add C7eff(p) and C8eff(p)
    tot_dict.update(get_C78eff(tot_dict, sector[:2]))
    return tot_dict

def RT(wc_obj, par, lep):
    mK = par['m_KL']
    mP = par['m_pi+']
    scale = par['m_rho0']
    ms = flavio.physics.running.running.get_ms(par, scale)
    wc = flavio.physics.bdecays.wilsoncoefficients.get_wceff_fccc(wc_obj, par, 'su', lep, lep, ms, scale, nf=3)
    ff = get_ff(0, par, 'KL')
    BT = ff['fT'] * 2 * mK / (mK + mP)  # convert between tensor FF conventions
    return -2 * BT / ff['f+'] * wc['tp'].real