How to use the flavio.physics.bdecays function in flavio

def test_bxclnu_lfu(self):
        self.assertAlmostEqual(
            flavio.sm_prediction('Rmue(B->Xclnu)'), 1, delta=0.01)
        self.assertAlmostEqual(
            flavio.sm_prediction('Rtaumu(B->Xclnu)'),
            flavio.sm_prediction('Rtaul(B->Xclnu)'),
            delta=0.01)
        self.assertEqual(
            flavio.physics.bdecays.bxlnu.BR_tot_leptonflavour(wc_obj, par, 'e', 'e'), 1)
        self.assertAlmostEqual(
            flavio.physics.bdecays.bxlnu.BR_tot_leptonflavour(wc_obj, par, 'tau', 'e'),
            1/flavio.physics.bdecays.bxlnu.BR_tot_leptonflavour(wc_obj, par, 'e', 'tau'))

def test_qcdf_contributions(self):
        # check that the individual contributions sum to the total contribution
        q2 = 4.679
        B = 'B0'
        V = 'K*0'
        lep = 'mu'
        amps_all = flavio.physics.bdecays.bvll.qcdf_interpolate.helicity_amps_qcdf(q2, par, B, V)
        amps_WA = flavio.physics.bdecays.bvll.qcdf_interpolate.helicity_amps_qcdf(q2, par, B, V, contribution='WA')
        amps_O8 = flavio.physics.bdecays.bvll.qcdf_interpolate.helicity_amps_qcdf(q2, par, B, V, contribution='O8')
        amps_QSS = flavio.physics.bdecays.bvll.qcdf_interpolate.helicity_amps_qcdf(q2, par, B, V, contribution='QSS')
        for i in amps_all.keys():
            if not amps_all[i] == 0:
                self.assertAlmostEqual(amps_all[i]/(amps_WA[i]+amps_O8[i]+amps_QSS[i]), 1, places=5)

def helicity_amps_qcdf(q2, wc, par, B, V, **kwargs):
    if q2 > 6:
        warnings.warn("The QCDF corrections should not be trusted for q2 above 6 GeV^2")
    mB = par['m_'+B]
    mV = par['m_'+V]
    X = sqrt(lambda_K(mB**2,q2,mV**2))/2.
    ta = transversity_amps_qcdf(q2, wc, par, B, V, **kwargs)
    h = flavio.physics.bdecays.angular.transversity_to_helicity(ta)
    return h

def _br_blnu(wc_obj, par, B, lep, nu):
    r"""Branching ratio of $B_q\to\ell^+\nu_\ell$."""
    bq = meson_quark[B]
    # CKM element
    if bq == 'bc':
        Vxb = flavio.physics.ckm.get_ckm(par)[1,2]
    elif bq == 'bu':
        Vxb = flavio.physics.ckm.get_ckm(par)[0,2]
    # renormalization scale
    scale = flavio.config['renormalization scale']['bll']
    # Wilson coefficients
    wc = wc_obj.get_wc(bq + lep + 'nu' + nu, scale, par)
    # add SM contribution to Wilson coefficient
    if lep == nu:
        wc['CVL_'+bq+lep+'nu'+nu] += flavio.physics.bdecays.wilsoncoefficients.get_CVLSM(par, scale, nf=5)
    mb = flavio.physics.running.running.get_mb(par, scale)
    if B == 'B+':
        mq = 0  # neglecting up quark mass
    elif B == 'Bc':
        mq = flavio.physics.running.running.get_mc(par, scale)
    return br_plnu_general(wc, par, Vxb, B, bq, lep, nu, mb, mq, delta=0)

def bvlilj_obs(function, q2, wc, par, B, V, l1, l2):
    ml1 = par['m_'+l1]
    ml2 = par['m_'+l2]
    mB = par['m_'+B]
    mV = par['m_'+V]
    if q2 < (ml1+ml2)**2 or q2 > (mB-mV)**2:
        return 0
    scale = flavio.config['renormalization scale']['bvll']
    mb = flavio.physics.running.running.get_mb(par, scale)
    h = helicity_amps(q2, wc, par, B, V, l1, l2)
    J = flavio.physics.bdecays.angular.angularcoeffs_general_v(h, q2, mB, mV, mb, 0, ml1, ml2)
    return function(J)

self.__init__()

        # Read the parameter metadata from the default YAML data file
        _read_yaml_object_metadata(pkgutil.get_data('flavio', 'data/parameters_metadata.yml'), self)

        # Read the uncorrelated parameter values from the default YAML data file
        _read_yaml_object_values(pkgutil.get_data('flavio', 'data/parameters_uncorrelated.yml'), self)

        # Read the correlated parameter values from the default YAML data file
        _read_yaml_object_values_correlated(pkgutil.get_data('flavio', 'data/parameters_correlated.yml'), self)

        # Read the parameters from the default PDG data file
        read_pdg(2018, self)

        # Read default parameters for B->V form factors
        flavio.physics.bdecays.formfactors.b_v.bsz_parameters.bsz_load('v2', 'LCSR', ('B->omega', 'B->rho'), self)
        flavio.physics.bdecays.formfactors.b_v.bsz_parameters.bsz_load('v2', 'LCSR-Lattice', ('B->K*', 'Bs->phi', 'Bs->K*'), self)

        # Read default parameters for B->P form factors
        flavio.physics.bdecays.formfactors.b_p.bsz_parameters.gkvd_load('v1', 'LCSR-Lattice', ('B->K', 'B->pi'), self)

        # Read default parameters for Lambdab->Lambda form factors
        flavio.physics.bdecays.formfactors.lambdab_12.lattice_parameters.lattice_load_ho(self)

def prefactor(q2, par, scale):
    xi_t = flavio.physics.ckm.xi('t','bs')(par)
    alphaem = flavio.physics.running.running.get_alpha(par, scale)['alpha_e']
    mLb = par['m_Lambdab']
    mL = par['m_Lambda']
    la_K = flavio.physics.bdecays.common.lambda_K(mLb**2, mL**2, q2)
    return par['GF'] * xi_t * alphaem * sqrt(q2) * la_K**(1/4.) / sqrt(3 * 2 * mLb**3 * pi**5) / 32.

def transversity_amps_qcdf(q2, wc, par, B, V, **kwargs):
    """QCD factorization corrections to B->Vll transversity amplitudes."""
    mB = par['m_'+B]
    mV = par['m_'+V]
    scale = config['renormalization scale']['bvll']
    # using the b quark pole mass here!
    mb = running.get_mb_pole(par)
    N = flavio.physics.bdecays.bvll.amplitudes.prefactor(q2, par, B, V)/4
    T_perp_ = T_perp(q2, par, wc, B, V, scale, **kwargs)
    T_para_ = T_para(q2, par, wc, B, V, scale, **kwargs)
    ta = {}
    ta['perp_L'] = N * sqrt(2)*2 * (mB**2-q2) * mb / q2 * T_perp_
    ta['perp_R'] =  ta['perp_L']
    ta['para_L'] = -ta['perp_L']
    ta['para_R'] =  ta['para_L']
    ta['0_L'] = ( N * mb * (mB**2 - q2)**2 )/(mB**2 * mV * sqrt(q2)) * T_para_
    ta['0_R'] = ta['0_L']
    ta['t'] = 0
    ta['S'] = 0
    return ta

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

def helicity_amps(q2, wc_obj, par, B, V, nu1, nu2):
    scale = flavio.config['renormalization scale']['bvll']
    label = flavio.physics.bdecays.common.meson_quark[(B,V)] + nu1 + nu2 # e.g. bsnuenue, bdnutaunumu
    wc = wc_obj.get_wc(label, scale, par)
    if nu1 == nu2: # add the SM contribution if neutrino flavours coincide
        wc['CL_'+label] += flavio.physics.bdecays.wilsoncoefficients.CL_SM(par)
    mB = par['m_'+B]
    mV = par['m_'+V]
    N = prefactor(q2, par, B, V)
    ff = get_ff(q2, par, B, V)
    wc_eff = flavio.physics.bdecays.wilsoncoefficients.get_wceff_nunu(q2, wc, par, B, V, nu1, nu2, scale)
    # below, mb is set to 4 just to save time. It doesn't enter at all as there
    # are no dipole operators.
    h = flavio.physics.bdecays.angular.helicity_amps_v(q2, mB, mV, 4, 0, 0, 0, ff, wc_eff, N)
    return h