How to use the flavio.classes function in flavio

def test_lfv(self):
        obs_1 = flavio.classes.Observable["BR(B0->K*emu)"]
        obs_2 = flavio.classes.Observable["BR(B0->K*mue)"]
        self.assertEqual(obs_1.prediction_central(flavio.default_parameters, wc_sm), 0)
        # BR(B->K*emu) should be 4 times larger as Wilson coeff is 2x the mue one
        self.assertAlmostEqual(
            obs_1.prediction_central(flavio.default_parameters, wc_np)
            /obs_2.prediction_central(flavio.default_parameters, wc_np),
            4.,  places=10)
        # test for errors
        flavio.sm_prediction("BR(B+->K*mue)")
        flavio.sm_prediction("BR(B0->rhotaue)")
        flavio.sm_prediction("BR(B+->rhotaumu)")
        flavio.sm_prediction("BR(Bs->phimutau)")

def test_lambdablambdall_subleading(self):
        ta_high = flavio.classes.AuxiliaryQuantity(
        'Lambdab->Lambdall subleading effects at high q2'
        ).prediction_central(par_nominal, wc_sm, q2=15, cp_conjugate=False)
        ta_low = flavio.classes.AuxiliaryQuantity(
        'Lambdab->Lambdall subleading effects at low q2'
        ).prediction_central(par_nominal, wc_sm, q2=1, cp_conjugate=False)
        # check that the keys contain all the transversity amps
        allamps = {('para0','L'), ('para1','L'), ('perp0','L'), ('perp1','L'),
                   ('para0','R'), ('para1','R'), ('perp0','R'), ('perp1','R')}
        self.assertEqual(set(ta_high.keys()), allamps)
        self.assertEqual(set(ta_low.keys()), allamps)
        # check that the central values are actually all zero

Defaults to 100. Less is faster but less precise.
    - `threads` (optional): if bigger than 1, number of threads to use for
      parallel computation of uncertainties
    - `plot_args` (optional): dictionary with keyword arguments to be passed
      to the matplotlib plot function, e.g. 'c' for colour.
    - `fill_args` (optional): dictionary with keyword arguments to be passed
      to the matplotlib fill_between function, e.g. 'facecolor'

    A word of caution regarding the `q2steps_err` option. By default, the
    uncertainty is only computed at 10 steps in $q^2$ and is interpolated in
    between. This can be enough if the uncertainty does not vary strongly
    with $q^2$. However, when the starting point or end point of the plot range
    is outside the physical phase space, the uncertainty will vanish at that
    point and the interpolation might be inaccurate.
    """
    obs = flavio.classes.Observable[obs_name]
    if obs.arguments != ['q2']:
        raise ValueError(r"Only observables that depend on $q^2$ (and nothing else) are allowed")
    step = (q2max-q2min)/(q2steps-1)
    q2_arr = np.arange(q2min, q2max+step, step)
    step = (q2max-q2min)/(q2steps_err-1)
    q2_err_arr = np.arange(q2min, q2max+step, step)
    # fix to avoid bounds_error in interp1d due to lack of numerical precision
    q2_err_arr[-1] = q2_arr[-1]
    if wc is None:
        wc = flavio.physics.eft._wc_sm # SM Wilson coefficients
        obs_err_arr = [flavio.sm_uncertainty(obs_name, q2, threads=threads) for q2 in q2_err_arr]
        obs_arr = [flavio.sm_prediction(obs_name, q2) for q2 in q2_arr]
    else:
        obs_err_arr = [flavio.np_uncertainty(obs_name, wc_obj=wc, q2=q2, threads=threads) for q2 in q2_err_arr]
        obs_arr = [flavio.np_prediction(obs_name, wc, q2) for q2 in q2_arr]
    ax = plt.gca()

# Observable and Prediction instances

_had = {'rho0': r'\rho^0', 'phi': r'\phi'}
_shortname = {'rho0': 'rho', 'phi': 'phi'}
_lep = {'e': ' e', 'mu': r'\mu',}

for V in _had:
    for l in _lep:
        _obs_name = "BR(tau->" + _shortname[V] + l + r")"
        _obs = flavio.classes.Observable(_obs_name)
        _process_tex = r"\tau^+\to " + _had[V] + _lep[l] + r"^+"
        _process_taxonomy = r'Process :: $\tau$ lepton decays :: LFV decays :: $\tau\to V\ell$ :: $' + _process_tex + r"$"
        _obs.add_taxonomy(_process_taxonomy)
        _obs.set_description(r"Branching ratio of $" + _process_tex + r"$")
        _obs.tex = r"$\text{BR}(" + _process_tex + r")$"
        flavio.classes.Prediction(_obs_name, br_tauvl_fct(V, l))

def _warn_meas_corr(self):
        """Warn the user if the fit contains multiple correlated measurements of
        an observable that is not included in the fit parameters, as this will
        lead to inconsistent results."""
        corr_with = {}
        # iterate over all measurements constraining at least one fit obs.
        for name in self.get_measurements:
            m = flavio.classes.Measurement[name]
            # iterate over all fit obs. constrained by this measurement
            for obs in set(self.observables) & set(m.all_parameters):
                # the constraint on this fit obs.
                constraint = m._parameters[obs][1]
                # find all the other obs. constrained by this constraint
                for c, p in m._constraints:
                    if c == constraint:
                        par = p
                        break
                for p in par:
                    # if the other obs. are not fit obs., append them to the list
                    if p not in self.observables:
                        if p not in corr_with:
                            corr_with[p] = [obs]
                        else:
                            corr_with[p].append(obs)

# Observable and Prediction instances

_lep = {'e': 'e', 'mu': r'\mu', 'tau': r'\tau'}
_had = {'D+': r'D^+', 'Ds': r'D_s',}

for D in _had:
  for l in _lep:
    _process_tex = _had[D]+r"\to "+_lep[l]+r"^+\nu_"+_lep[l]
    _process_taxonomy = r'Process :: $c$ hadron decays :: Leptonic tree-level decays :: $D\to \ell\nu$ :: $'

    _obs_name = "BR("+D+"->"+l+"nu)"
    _obs = flavio.classes.Observable(_obs_name)
    _obs.set_description(r"Branching ratio of $"+_process_tex+r"$")
    _obs.tex = r"$\text{BR}("+_process_tex+r")$"
    _obs.add_taxonomy(_process_taxonomy + _process_tex + r'$')
    flavio.classes.Prediction(_obs_name, br_Dlnu_fct(D, l))

_hadr = {
'B0->K': {'tex': r"B^0\to K^0", 'B': 'B0', 'P': 'K0', },
'B+->K': {'tex': r"B^+\to K^+", 'B': 'B+', 'P': 'K+', },
'B0->pi': {'tex': r"B^0\to \pi^0", 'B': 'B0', 'P': 'pi0', },
'B+->pi': {'tex': r"B^+\to \pi^+", 'B': 'B+', 'P': 'pi+', },
}

for M in _hadr.keys():

    _process_tex = _hadr[M]['tex'] + r"\nu\bar\nu"
    _process_taxonomy = r'Process :: $b$ hadron decays :: FCNC decays :: $B\to P\nu\bar\nu$ :: $'  + _process_tex + r'$'

    # binned branching ratio
    _obs_name = "("+M+"nunu)"
    _obs = flavio.classes.Observable(name=_obs_name, arguments=['q2min', 'q2max'])
    _obs.set_description(r"Binned differential branching ratio of $" + _process_tex + r"$")
    _obs.tex = r"$\langle \frac{d\text{BR}}{dq^2} \rangle(" + _process_tex + r")$"
    _obs.add_taxonomy(_process_taxonomy)
    flavio.classes.Prediction(_obs_name, dbrdq2_int_summed(_hadr[M]['B'], _hadr[M]['P']))

    # differential branching ratio
    _obs_name = "dBR/dq2("+M+"nunu)"
    _obs = flavio.classes.Observable(name=_obs_name, arguments=['q2'])
    _process_tex = _hadr[M]['tex'] + r"\nu\bar\nu"
    _obs.set_description(r"Differential branching ratio of $" + _process_tex + r"$")
    _obs.tex = r"$\frac{d\text{BR}}{dq^2}(" + _process_tex + r")$"
    _obs.add_taxonomy(_process_taxonomy)
    flavio.classes.Prediction(_obs_name, dbrdq2_summed(_hadr[M]['B'], _hadr[M]['P']))

    # total branching ratio
    _obs_name = "BR("+M+"nunu)"

def make_observable(l, ltex, lname, lfunc):
    """Instantiate the observable and add metadata."""
    obs_name = "a_" + l
    obs = flavio.classes.Observable(obs_name)
    obs.set_description(r"Anomalous magnetic moment of the " + lname)
    obs.tex = r"$a_{}$".format(ltex)
    obs.add_taxonomy(r"Process :: Dipole moments :: Lepton anomalous magnetic moments :: $a_{}$".format(ltex))
    flavio.classes.Prediction(obs_name, lfunc)

_obs.set_description(r"Binned differential branching ratio of $" + _process_tex + r"$")
    _obs.tex = r"$\langle \frac{d\text{BR}}{dq^2} \rangle(" + _process_tex + r")$"
    _obs.add_taxonomy(_process_taxonomy + _process_tex + r"$")
    flavio.classes.Prediction(_obs_name, dbrdq2_int_summed(_hadr[M]['B'], _hadr[M]['V']))

    # differential branching ratio
    _obs_name = "dBR/dq2("+M+"nunu)"
    _obs = flavio.classes.Observable(name=_obs_name, arguments=['q2'])
    _obs.set_description(r"Differential branching ratio of $" + _process_tex + r"$")
    _obs.tex = r"$\frac{d\text{BR}}{dq^2}(" + _process_tex + r")$"
    _obs.add_taxonomy(_process_taxonomy + _process_tex + r"$")
    flavio.classes.Prediction(_obs_name, dbrdq2_summed(_hadr[M]['B'], _hadr[M]['V']))

    # total branching ratio
    _obs_name = "BR("+M+"nunu)"
    _obs = flavio.classes.Observable(name=_obs_name)
    _obs.set_description(r"Branching ratio of $" + _process_tex + r"$")
    _obs.tex = r"$\text{BR}(" + _process_tex + r")$"
    _obs.add_taxonomy(_process_taxonomy + _process_tex + r"$")
    flavio.classes.Prediction(_obs_name, BRtot_summed(_hadr[M]['B'], _hadr[M]['V']))

    if 'K*' in M: # FL only implemented for K*

        # binned FL
        _obs_name = "("+M+"nunu)"
        _obs = flavio.classes.Observable(name=_obs_name, arguments=['q2min', 'q2max'])
        _obs.set_description(r"Binned longitudinal polarization fraction of $" + _process_tex + r"$")
        _obs.tex = r"$\langle F_L \rangle(" + _process_tex + r")$"
        _obs.add_taxonomy(_process_taxonomy + _process_tex + r"$")
        flavio.classes.Prediction(_obs_name, fl_int_summed(_hadr[M]['B'], _hadr[M]['V']))

        # differential FL

_obs.set_description(r"Ratio of $Z^0$ partial widths to ${}$ pairs vs. all hadrons".format(tex))
            flavio.classes.Prediction(_obs_name, Rq(f))

_obs_name = "R_uc"
_obs = flavio.classes.Observable(_obs_name)
_obs.tex = r"$R_{uc}^0$"
_obs.add_taxonomy(r'Process :: $Z^0$ decays :: Flavour conserving decays :: $' + _process_tex + r"$")
_obs.set_description(r"Average ratio of $Z^0$ partial widths to $u$ or $c$ pairs vs. all hadrons")
flavio.classes.Prediction(_obs_name, Rq1q2('u', 'c'))

# LFV Z decays
for (f1, f2) in [('e', 'mu'), ('e', 'tau'), ('mu', 'tau'), ]:
    tex1 = _leptons[f1]
    tex2 = _leptons[f2]
    _obs_name = "BR(Z->{}{})".format(f1, f2)
    _obs = flavio.classes.Observable(_obs_name)
    _process_tex = r"Z^0\to {}^\pm{}^\mp".format(tex1, tex2)
    _obs.tex = r"$\text{BR}(" + _process_tex + r")$"
    _obs.set_description(r"Branching ratio of $" + _process_tex + r"$")
    _obs.add_taxonomy(r'Process :: $Z^0$ decays :: FCNC decays :: $' + _process_tex + r"$")
    flavio.classes.Prediction(_obs_name, BRZ_fct_av(f1, f2))


_obs_name = "GammaZ"
_obs = flavio.classes.Observable(_obs_name)
_obs.tex = r"$\Gamma_Z$"
_obs.set_description(r"Total width of the $Z^0$ boson")
for _f in (_leptons, _uquarks, _dquarks):
    for f, tex in _f.items():
        if f in _leptons:
            _process_tex = r"Z^0\to {}^+{}^-".format(tex, tex)
        else: