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

# parameters taken from PDG and table I of 1311.0903
c.set_constraint('alpha_s', '0.1184(7)')
c.set_constraint('f_Bs', '0.2277(45)')
c.set_constraint('f_B0', '0.1905(42)')
c.set_constraint('Vcb', 4.24e-2)
c.set_constraint('Vub', 3.82e-3)
c.set_constraint('delta', radians(73.))
c.set_constraint('DeltaGamma/Gamma_Bs', 0.1226)

par = c.get_central_all()


wc_obj = WilsonCoefficients()
wc = wctot_dict(wc_obj, 'bsmumu', 4.8, par)
wc_e = wctot_dict(wc_obj, 'bsee', 4.8, par)
wc_tau = wctot_dict(wc_obj, 'bstautau', 4.8, par)


class TestBll(unittest.TestCase):
    def test_bsll(self):
        # just some trivial tests to see if calling the functions raises an error
        self.assertGreater(br_lifetime_corr(0.08, -1), 0)
        self.assertEqual(len(amplitudes(par, wc, 'Bs', 'mu', 'mu')), 2)
        # ADeltaGamma should be +1.0 in the SM
        self.assertEqual(ADeltaGamma(par, wc, 'Bs', 'mu'), 1.0)
        self.assertEqual(flavio.sm_prediction('ADeltaGamma(Bs->mumu)'), 1.0)
        # BR should be around 3.5e-9
        self.assertAlmostEqual(br_inst(par, wc, 'Bs', 'mu', 'mu')*1e9, 3.5, places=0)
        # correction factor should enhance the BR by roughly 7%
        self.assertAlmostEqual(br_timeint(par, wc, 'Bs', 'mu', 'mu')/br_inst(par, wc, 'Bs', 'mu', 'mu'), 1.07, places=2)
        # ratio of Bs->mumu and Bs->ee BRs should be roughly given by ratio of squared masses
        self.assertAlmostEqual(

c = copy.deepcopy(default_parameters)
# parameters taken from PDG and table I of 1311.0903
c.set_constraint('alpha_s', '0.1184(7)')
c.set_constraint('f_Bs', '0.2277(45)')
c.set_constraint('f_B0', '0.1905(42)')
c.set_constraint('Vcb', 4.24e-2)
c.set_constraint('Vub', 3.82e-3)
c.set_constraint('delta', radians(73.))
c.set_constraint('DeltaGamma/Gamma_Bs', 0.1226)

par = c.get_central_all()


wc_obj = WilsonCoefficients()
wc = wctot_dict(wc_obj, 'bsmumu', 4.8, par)
wc_e = wctot_dict(wc_obj, 'bsee', 4.8, par)
wc_tau = wctot_dict(wc_obj, 'bstautau', 4.8, par)


class TestBll(unittest.TestCase):
    def test_bsll(self):
        # just some trivial tests to see if calling the functions raises an error
        self.assertGreater(br_lifetime_corr(0.08, -1), 0)
        self.assertEqual(len(amplitudes(par, wc, 'Bs', 'mu', 'mu')), 2)
        # ADeltaGamma should be +1.0 in the SM
        self.assertEqual(ADeltaGamma(par, wc, 'Bs', 'mu'), 1.0)
        self.assertEqual(flavio.sm_prediction('ADeltaGamma(Bs->mumu)'), 1.0)
        # BR should be around 3.5e-9
        self.assertAlmostEqual(br_inst(par, wc, 'Bs', 'mu', 'mu')*1e9, 3.5, places=0)
        # correction factor should enhance the BR by roughly 7%
        self.assertAlmostEqual(br_timeint(par, wc, 'Bs', 'mu', 'mu')/br_inst(par, wc, 'Bs', 'mu', 'mu'), 1.07, places=2)

def test_wctot(self):
        wc_low_correct = np.array([ -2.93671059e-01,   1.01676402e+00,  -5.87762813e-03,
        -8.70666812e-02,   4.11098919e-04,   1.10641294e-03,
        -2.95662859e-01,  -1.63048361e-01,   4.11363023e+00,
        -4.19345312e+00,   3.43507549e-03,   1.22202095e-03,
        -1.03192325e-03,  -1.00703396e-04,  -3.17810374e-03])
        wc_obj = eft.WilsonCoefficients()
        wc_low = wilsoncoefficients.wctot_dict(wc_obj, 'bsmumu', 4.2, par)
        wc_names = ['C1_bs', 'C2_bs', 'C3_bs', 'C4_bs', 'C5_bs', 'C6_bs', 'C7_bs', 'C8_bs', 'C9_bsmumu', 'C10_bsmumu', 'C3Q_bs', 'C4Q_bs', 'C5Q_bs', 'C6Q_bs', 'Cb_bs', 'C1p_bs', 'C2p_bs', 'C3p_bs', 'C4p_bs', 'C5p_bs', 'C6p_bs', 'C7p_bs', 'C8p_bs', 'C9p_bsmumu', 'C10p_bsmumu', 'C3Qp_bs', 'C4Qp_bs', 'C5Qp_bs', 'C6Qp_bs', 'Cbp_bs', 'CS_bsmumu', 'CP_bsmumu', 'CSp_bsmumu', 'CPp_bsmumu']
        wc_low_array = np.asarray([wc_low[key] for key in wc_names])
        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_low_array[6] = wc_low_array[6] + np.dot(yi, wc_low_array[:6]) # c7eff
        wc_low_array[7] = wc_low_array[7] + np.dot(zi, wc_low_array[:6]) # c8eff
        np.testing.assert_almost_equal(wc_low_array[:15], wc_low_correct, decimal=2)

def test_acp(self):
        # check that the SM central values for the individual B->Xsgamma
        # and B->Xdgamma (ignoring long distance contributions) roughly
        # agree with the values quoted in hep-ph/0312260
        wc_sm_s = flavio.physics.bdecays.wilsoncoefficients.wctot_dict(wc_obj, 'bsee', scale=2, par=par_dict, nf_out=5)
        wc_sm_d = flavio.physics.bdecays.wilsoncoefficients.wctot_dict(wc_obj, 'bdee', scale=2, par=par_dict, nf_out=5)
        p_ave_s = flavio.physics.bdecays.bxgamma.PE0_BR_BXgamma(wc_sm_s, par_dict, 's', 1.6)
        p_asy_s = flavio.physics.bdecays.bxgamma.PE0_ACP_BXgamma(wc_sm_s, par_dict, 's', 1.6)
        acp_s = p_asy_s/p_ave_s
        self.assertAlmostEqual(100*acp_s, 0.44, delta=0.5)
        p_ave_d = flavio.physics.bdecays.bxgamma.PE0_BR_BXgamma(wc_sm_d, par_dict, 'd', 1.6)
        p_asy_d = flavio.physics.bdecays.bxgamma.PE0_ACP_BXgamma(wc_sm_d, par_dict, 'd', 1.6)
        acp_d = p_asy_d/p_ave_d
        # check that the s+d CP asymmetry vanishes
        self.assertAlmostEqual(flavio.sm_prediction('ACP(B->Xgamma)'), 0, delta=1e-9)

def function(wc_obj, par_dict, q2, cp_conjugate):
        par = par_dict.copy()
        if cp_conjugate:
            par = conjugate_par(par)
        wc = wctot_dict(wc_obj, label, scale, par)
        if cp_conjugate:
            wc = conjugate_wc(wc)
        return flavio.physics.bdecays.bvll.qcdf.helicity_amps_qcdf(q2, wc, par, B, V)
    return function

def __init__(self, B, V, lep, wc_obj, par):
        """Initialize the class and cache results needed more often."""
        self.B = B
        self.V = V
        self.lep = lep
        self.wc_obj = wc_obj
        self.par = par
        self.par_conjugate = conjugate_par(par)
        self.prefactor = prefactor(None, self.par, B, V)
        self.prefactor_conjugate = prefactor(None, self.par_conjugate, B, V)
        self.scale = config['renormalization scale']['bvll']
        self.label = meson_quark[(B,V)] + lep + lep # e.g. bsmumu, bdtautau
        self.wctot_dict = wctot_dict(wc_obj, self.label, self.scale, par)
        self._ff = {}
        self._wceff = {}
        self._wceff_bar = {}
        self._ha = {}
        self._ha_bar = {}
        self._j = {}
        self._j_bar = {}
        self.ml = par['m_'+lep]
        self.mB = par['m_'+B]
        self.mV = par['m_'+V]
        self.mb = running.get_mb(par, self.scale)

def ADG_func(wc_obj, par):
        scale = config['renormalization scale']['bll']
        label = meson_quark[B]+lep+lep
        wc = wctot_dict(wc_obj, label, scale, par)
        return ADeltaGamma(par, wc, B, lep)
    return ADG_func

def helicity_amps_ff(q2, wc_obj, par_dict, B, P, lep, cp_conjugate):
    par = par_dict.copy()
    if cp_conjugate:
        par = conjugate_par(par)
    scale = config['renormalization scale']['bpll']
    label = meson_quark[(B,P)] + lep + lep # e.g. bsmumu, bdtautau
    wc = wctot_dict(wc_obj, label, scale, par)
    if cp_conjugate:
        wc = conjugate_wc(wc)
    wc_eff = get_wceff(q2, wc, par, B, P, lep, scale)
    ml = par['m_'+lep]
    mB = par['m_'+B]
    mP = par['m_'+P]
    mb = running.get_mb(par, scale)
    N = prefactor(q2, par, B, P)
    ff = get_ff(q2, par, B, P)
    h = angular.helicity_amps_p(q2, mB, mP, mb, 0, ml, ml, ff, wc_eff, N)
    return h

def ACPBXgamma(wc_obj, par, E0):
    r"""Direct CP asymmetry of $B\to X_{s+d}\gamma$. `E0` is the photon energy
    cutoff $E_0$ in GeV (currently works only for `E0=1.6`)."""
    scale = flavio.config['renormalization scale']['bxgamma']
    # these are the b->qee Wilson coefficients - they contain the b->qgamma ones as a subset
    xi_t_d = flavio.physics.ckm.xi('t', 'bd')(par)
    xi_t_s = flavio.physics.ckm.xi('t', 'bs')(par)
    wc_d = flavio.physics.bdecays.wilsoncoefficients.wctot_dict(wc_obj, 'bdee', scale, par, nf_out=5)
    wc_s = flavio.physics.bdecays.wilsoncoefficients.wctot_dict(wc_obj, 'bsee', scale, par, nf_out=5)
    br_d = abs(xi_t_d)**2 * PE0_BR_BXgamma(wc_d, par, 'd', E0)
    br_s = abs(xi_t_s)**2 * PE0_BR_BXgamma(wc_s, par, 's', E0)
    as_d = abs(xi_t_d)**2 * PE0_ACP_BXgamma(wc_d, par, 'd', E0)
    as_s = abs(xi_t_s)**2 * PE0_ACP_BXgamma(wc_s, par, 's', E0)
    # return (as_s)/(br_s + br_d)
    return (as_s + as_d)/(br_s + br_d)

def get_wceff(q2, wc_obj, par, B, V, lep):
    scale = config['bdecays']['scale_bvll']
    wc = wctot_dict(wc_obj, meson_quark[(B,V)]+lep+lep, scale, par)
    xi_u = ckm.xi('u',meson_quark[(B,V)])(par)
    xi_t = ckm.xi('t',meson_quark[(B,V)])(par)
    qiqj=meson_quark[(B,V)]
    Yq2 = matrixelements.Y(q2, wc, par, scale, qiqj) + (xi_u/xi_t)*matrixelements.Yu(q2, wc, par, scale, qiqj)
        #   b) NNLO Q1,2
    delta_C7 = matrixelements.delta_C7(par=par, wc=wc, q2=q2, scale=scale, qiqj=qiqj)
    delta_C9 = matrixelements.delta_C9(par=par, wc=wc, q2=q2, scale=scale, qiqj=qiqj)
    mb = running.get_mb(par, scale)
    ll = lep + lep
    c = {}
    c['7']  = wc['C7eff_'+qiqj]      + delta_C7
    c['7p'] = wc['C7effp_'+qiqj]
    c['v']  = wc['C9_'+qiqj+ll]      + delta_C9 + Yq2
    c['vp'] = wc['C9p_'+qiqj+ll]
    c['a']  = wc['C10_'+qiqj+ll]
    c['ap'] = wc['C10p_'+qiqj+ll]