docs: def m_b(mbmb, scale, f, alphasMZ=0.1185, loop=3):
      r"""Get running b quark mass in the MSbar scheme at the scale `scale`
      in the theory with `f` dynamical quark flavours starting from $m_b(m_b)$"""
      if scale == mbmb and f == 5:
          return mbmb  # nothing to do
      _sane(scale, f)
      alphas_mb = alpha_s(mbmb, 5, alphasMZ=alphasMZ, loop=loop)
      crd = rundec.CRunDec()
      if f == 5:
          alphas_scale = alpha_s(scale, f, alphasMZ=alphasMZ, loop=loop)
          return crd.mMS2mMS(mbmb, alphas_mb, alphas_scale, f, loop)
      elif f == 4:
          crd.nfMmu.Mth = 4.8
          crd.nfMmu.muth = 4.8
          crd.nfMmu.nf = 5
          return crd.mH2mL(mbmb, alphas_mb, mbmb, crd.nfMmu, scale, loop)
      elif f == 3:
          mc = 1.3
          crd.nfMmu.Mth = 4.8
          crd.nfMmu.muth = 4.8
          crd.nfMmu.nf = 5
          mbmc = crd.mH2mL(mbmb, alphas_mb, mbmb, crd.nfMmu, mc, loop)
          crd.nfMmu.Mth = mc
          crd.nfMmu.muth = mc
          crd.nfMmu.nf = 4
          alphas_mc = alpha_s(mc, 4, alphasMZ=alphasMZ, loop=loop)
          return crd.mH2mL(mbmc, alphas_mc, mc, crd.nfMmu, scale, loop)
      elif f == 6:
          crd.nfMmu.Mth = 170
          crd.nfMmu.muth = 170
          crd.nfMmu.nf = 6
          return crd.mL2mH(mbmb, alphas_mb, mbmb, crd.nfMmu, scale, loop)
      else:
          raise ValueError("Invalid input: f={}, scale={}".format(f, scale))