From a0ee861a739e736f16ca92c59eca736a679eb160 Mon Sep 17 00:00:00 2001 From: Jorge-Alda Date: Tue, 26 Sep 2023 12:08:28 +0200 Subject: [PATCH 1/2] Angular observables in b->clnu and LFU differences in B0_>D*lnu --- flavio/physics/bdecays/bvlnu.py | 245 +++++++++++++++++++++++++++++++- 1 file changed, 244 insertions(+), 1 deletion(-) diff --git a/flavio/physics/bdecays/bvlnu.py b/flavio/physics/bdecays/bvlnu.py index 410546f4b..9efa31ecb 100644 --- a/flavio/physics/bdecays/bvlnu.py +++ b/flavio/physics/bdecays/bvlnu.py @@ -1,7 +1,8 @@ r"""Functions for exclusive $B\to V\ell\nu$ decays.""" -from math import sqrt, pi, cos, sin +from math import sqrt, pi, cos, sin, log import flavio +from flavio.physics.common import conjugate_par, conjugate_wc from flavio.physics.bdecays.common import lambda_K, meson_quark, meson_ff from flavio.physics.bdecays.wilsoncoefficients import wctot_dict from flavio.physics import ckm @@ -78,6 +79,54 @@ def _get_angularcoeff(q2, wc_obj, par, B, V, lep, nu): J = angular.angularcoeffs_general_v(h, q2, mB, mV, mb, mlight, ml, 0) return J + +def get_angularcoeff_cp(q2, wc_obj, par, B, V, lep): + Jlist = [_get_angularcoeff_cp(q2, wc_obj, par, B, V, lep, nu) + for nu in ['e', 'mu', 'tau']] + J = {} + J['1s'] = sum([JJ['1s'] for JJ in Jlist]) + J['1c'] = sum([JJ['1c'] for JJ in Jlist]) + J['2s'] = sum([JJ['2s'] for JJ in Jlist]) + J['2c'] = sum([JJ['2c'] for JJ in Jlist]) + J['6s'] = sum([JJ['6s'] for JJ in Jlist]) + J['6c'] = sum([JJ['6c'] for JJ in Jlist]) + J[3] = sum([JJ[3] for JJ in Jlist]) + J[4] = sum([JJ[4] for JJ in Jlist]) + J[5] = sum([JJ[5] for JJ in Jlist]) + J[7] = sum([JJ[7] for JJ in Jlist]) + J[8] = sum([JJ[8] for JJ in Jlist]) + J[9] = sum([JJ[9] for JJ in Jlist]) + return J + + +def _get_angularcoeff_cp(q2, wc_obj, par, B, V, lep, nu): + scale = config['renormalization scale']['bvll'] + par_conjugate = conjugate_par(par) + mb = running.get_mb(par_conjugate, scale) + wc = get_wceff_fccc(wc_obj, par_conjugate, + meson_quark[(B, V)], lep, nu, mb, scale, nf=5) + wc = {k: v.conjugate() for k, v in wc.items()} + if lep != nu and all(C == 0 for C in wc.values()): + # if all WCs vanish, so does the AC! + return {k: 0 for k in + ['1s', '1c', '2s', '2c', '6s', '6c', 3, 4, 5, 7, 8, 9]} + ml = par_conjugate['m_'+lep] + mB = par_conjugate['m_'+B] + mV = par_conjugate['m_'+V] + qi_qj = meson_quark[(B, V)] + if qi_qj == 'bu': + mlight = 0. # neglecting the up quark mass + if qi_qj == 'bc': + # this is needed for scalar contributions + mlight = running.get_mc(par_conjugate, scale) + N = prefactor(q2, par_conjugate, B, V, lep) + ff = get_ff(q2, par_conjugate, B, V) + h = angular.helicity_amps_v( + q2, mB, mV, mb, mlight, ml, 0, ff, wc, N) + J = angular.angularcoeffs_general_v( + h, q2, mB, mV, mb, mlight, ml, 0) + return J + def dGdq2(J): r"""$q^2$-differential branching ratio in terms of angular coefficients.""" return 3/4. * (2 * J['1s'] + J['1c']) - 1/4. * (2 * J['2s'] + J['2c']) @@ -453,6 +502,33 @@ def BR_binned_leptonflavour_function(B, V, lnum, lden, A): return lambda wc_obj, par, q2min, q2max: BR_binned_leptonflavour(q2min, q2max, wc_obj, par, B, V, lnum, lden, A) +def S_theory_num(J, J_bar, i): + return J[i] + J_bar[i] + + +def A_theory_num(J, J_bar, i): + return J[i] - J_bar[i] + + +def S_experiment_num(J, J_bar, i): + if i in [4, '6s', '6c', 7, 9]: + return -S_theory_num(J, J_bar, i) + return S_theory_num(J, J_bar, i) + + +def A_experiment_num(J, J_bar, i): + if i in [4, '6s', '6c', 7, 9]: + return -A_theory_num(J, J_bar, i) + return A_theory_num(J, J_bar, i) + + +def AFB_experiment_num(J, J_bar): + return 3/4.*S_experiment_num(J, J_bar, '6s') + + +def SA_den(J, J_bar): + return dGdq2(J) + dGdq2(J_bar) + # Observable and Prediction instances _tex = {'e': 'e', 'mu': '\mu', 'tau': r'\tau', 'l': r'\ell'} @@ -602,3 +678,170 @@ def BR_binned_leptonflavour_function(B, V, lnum, lden, A): _process_tex = _hadr[M]['tex'] + r"\tau^+\nu" _obs.add_taxonomy(_process_taxonomy + _process_tex + r"$") Prediction(_obs_name, BR_binned_tot_function('B0', 'D*+', 'tau', A=None)) + +# Angular observables +_angular_obs = {'S3': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 3), 'tex': r'S_3', 'desc': 'CP-averaged angular observable'}, + 'S4': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 4), 'tex': r'S_4', 'desc': 'CP-averaged angular observable'}, + 'S5': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 5), 'tex': r'S_5', 'desc': 'CP-averaged angular observable'}, + 'S6c': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, '6c'), 'tex': r'S_6^c', 'desc': 'CP-averaged angular observable'}, + 'S7': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 7), 'tex': r'S_7', 'desc': 'CP-averaged angular observable'}, + 'S8': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 8), 'tex': r'S_8', 'desc': 'CP-averaged angular observable'}, + 'S9': {'func_num': lambda J, J_bar: S_experiment_num(J, J_bar, 9), 'tex': r'S_9', 'desc': 'CP-averaged angular observable'}, + 'A3': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 3), 'tex': r'A_3', 'desc': 'Angular CP asymmetry'}, + 'A4': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 4), 'tex': r'A_4', 'desc': 'Angular CP asymmetry'}, + 'A5': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 5), 'tex': r'A_5', 'desc': 'Angular CP asymmetry'}, + 'A6s': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, '6s'), 'tex': r'A_6^s', 'desc': 'Angular CP asymmetry'}, + 'A7': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 7), 'tex': r'A_7', 'desc': 'Angular CP asymmetry'}, + 'A8': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 8), 'tex': r'A_8', 'desc': 'Angular CP asymmetry'}, + 'A9': {'func_num': lambda J, J_bar: A_experiment_num(J, J_bar, 9), 'tex': r'A_9', 'desc': 'Angular CP asymmetry'}, + 'AFB': {'func_num': AFB_experiment_num, 'tex': r'A_\text{FB}', 'desc': 'forward-backward asymmetry'}, + } + + +def make_metadata_binned(M, l, obs, obsdict): + _process_tex = _hadr[M]['tex'] + _tex[l]+r"^+\nu_"+_tex[l] + B = _hadr[M]['B'] + V = _hadr[M]['V'] + _obs_name = "<" + obs + ">("+M+l+"nu)" + _obs = Observable(name=_obs_name, arguments=['q2min', 'q2max']) + _obs.set_description( + 'Binned ' + obsdict['desc'] + r" in $" + _process_tex + r"$") + _obs.tex = r"$\langle " + obsdict['tex'] + \ + r"\rangle(" + _process_tex + r")$" + _obs.add_taxonomy(_process_taxonomy) + return _obs + + +def make_metadata_differential(M, l, obs, obsdict): + _process_tex = _hadr[M]['tex'] + _tex[l]+r"^+\nu_"+_tex[l] + B = _hadr[M]['B'] + V = _hadr[M]['V'] + _obs_name = obs + "("+M+l+"nu)" + _obs = Observable(name=_obs_name, arguments=['q2']) + _obs.set_description(obsdict['desc'][0].capitalize( + ) + obsdict['desc'][1:] + r" in $" + _process_tex + r"$") + _obs.tex = r"$" + obsdict['tex'] + r"(" + _process_tex + r")$" + _obs.add_taxonomy(_process_taxonomy) + return _obs + + +def nintegrate_pole(function, q2min, q2max, epsrel=0.005): + # this is a special integration function to treat the presence of the + # photon pole at low q^2. If q2min is below 0.1 GeV^2, it adds and subtracts + # the 1/q^2-enhanced pole part to split the integral into a well-behaved part + # and one that is trivially solved analytically. + # This leads to a huge speed-up. + if q2min <= 0.1 and q2min > 0: + q20 = q2min + f_q20 = function(q20) + int_a = flavio.math.integrate.nintegrate( + lambda q2: function(q2)-f_q20*q20/q2, q2min, q2max, epsrel=epsrel) + int_b = f_q20*q20 * log(q2max/q2min) + return int_a + int_b + else: + return flavio.math.integrate.nintegrate(function, q2min, q2max) + + +def BVlnu_int_ratio(func_num, func_den, q2min, q2max, B, V, l, wc_obj, par): + epsrel = 0.005 + ml = par['m_'+l] + mB = par['m_'+B] + mV = par['m_'+V] + q2min_allowed = max(ml**2, q2min) + q2max_allowed = min((mB - mV)**2, q2max) + if q2max_allowed <= q2min_allowed: + return 0 + + def J(q2): + return get_angularcoeff(q2, wc_obj, par, B, V, l) + + def Jbar(q2): + return get_angularcoeff_cp(q2, wc_obj, par, B, V, l) + + def obs_num(q2): + return func_num(J(q2), Jbar(q2)) + + def obs_den(q2): + return func_den(J(q2), Jbar(q2)) + + num = nintegrate_pole(obs_num, q2min_allowed, q2max_allowed, epsrel=epsrel) + if num == 0: + return 0 + den = nintegrate_pole(obs_den, q2min_allowed, q2max_allowed, epsrel=epsrel) + return num / den + + +def BVlnu_ratio(func_num, func_den, q2, B, V, l, wc_obj, par): + ml = par['m_'+l] + mB = par['m_'+B] + mV = par['m_'+V] + if q2 < ml**2 or q2 > (mB - mV)**2: + return 0 + J = get_angularcoeff(q2, wc_obj, par, B, V, l) + Jbar = get_angularcoeff_cp(q2, wc_obj, par, B, V, l) + + num = func_num(J, Jbar) + if num == 0: + return 0 + + den = func_den(J, Jbar) + return num / den + + +def make_obs(M, l, obs, obsdict): + B = _hadr[M]['B'] + V = _hadr[M]['V'] + func_num = obsdict['func_num'] + + # binned angular observables + _obs = make_metadata_binned(M, l, obs, obsdict) + + def func(wc_obj, par, q2min, q2max): return BVlnu_int_ratio( + func_num, SA_den, q2min, q2max, B, V, l, wc_obj, par) + Prediction(_obs.name, func) + + # differential angular observables + _obs = make_metadata_differential(M, l, obs, obsdict) + + def func(wc_obj, par, q2): return BVlnu_ratio( + func_num, SA_den, q2, B, V, l, wc_obj, par) + Prediction(_obs.name, func) + + +for l in ['e', 'mu', 'tau']: + for M in _hadr.keys(): + for obs, obsdict in _angular_obs.items(): + # angular obs + make_obs(M, l, obs, obsdict) + + +def diff(x, y): + return x-y + + +for D_obs in ['S3', 'S5', 'S7', 'S9', 'AFB']: + tex = _angular_obs[D_obs]['tex'] + obs = Observable.from_function('Dmue_{}(B0->D*lnu)'.format(D_obs), + ['{}(B0->D*munu)'.format(D_obs), + '{}(B0->D*enu)'.format(D_obs)], + diff) + obs.set_description( + r"Difference of angular observable ${}$ in $B^0\to D^{{\ast -}}\mu^+\nu$ and $B^0\to D^{{\ast -}}e^+\nu$".format(tex)) + obs.tex = r"$D_{{{}}}^{{\mu e}}(B^0\to D^{{\ast -}}\ell^+\nu)$".format(tex) + obs.add_taxonomy( + r"Process :: $b$ hadron decays :: Semi-leptonic tree-level decays :: $B\to V\ell\nu$ :: $B^0\to K^{\ast -}e^+\nu_e$") + obs.add_taxonomy( + r"Process :: $b$ hadron decays :: Semi-leptonic tree-level decays :: $B\to V\ell\nu$ :: $B^0\to K^{\ast -}\mu^+\nu_\mu$") + + obs = Observable.from_function('(B0->D*lnu)'.format(D_obs), + ['<{}>(B0->D*munu)'.format(D_obs), + '<{}>(B0->D*enu)'.format(D_obs)], + diff) + obs.set_description( + r"Binned difference of angular observable ${}$ in $B^0\to D^{{\ast -}}\mu^+\nu$ and $B^0\to K^{{\ast -}}e^+\nu$".format(tex)) + obs.tex = r"$\langle D_{{{}}}^{{\mu e}} \rangle(B^0\to D^{{\ast -}}\ell^+\nu)$".format( + tex) + obs.add_taxonomy( + r"Process :: $b$ hadron decays :: Semi-leptonic tree-level decays :: $B\to V\ell\nu$ :: $B^0\to K^{\ast -}e^+\nu_e$") + obs.add_taxonomy( + r"Process :: $b$ hadron decays :: Semi-leptonic tree-level decays :: $B\to V\ell\nu$ :: $B^0\to K^{\ast -}\mu^+\nu_\mu$") From d725ed88bf66f86bed1bfe1800a1cec5985cb1e5 Mon Sep 17 00:00:00 2001 From: AlejandroMirRamos <83365676+AlejandroMirRamos@users.noreply.github.com> Date: Thu, 28 Sep 2023 16:21:21 +0200 Subject: [PATCH 2/2] Belle-II measurements for B0->D*lnu observables --- flavio/data/measurements.yml | 65 ++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/flavio/data/measurements.yml b/flavio/data/measurements.yml index 4e7ff38ba..a58e6eadb 100644 --- a/flavio/data/measurements.yml +++ b/flavio/data/measurements.yml @@ -22430,3 +22430,68 @@ LHCb Bs->K*mumu 2018: inspire: LHCb:2018rym values: BR_LHCb(Bs->K*0mumu): 2.9 ± 1.0 ± 0.2 ± 0.3 1e-8 + +Belle-II B0->D*lnu angular 2023: + inspire: Belle-II:2023kkm + experiment: Belle-II + values: + - name: (B0->D*lnu) + q2min: 4.85 + q2max: 10.689 + value: 0.099 ± 0.056 ± 0.020 + - name: (B0->D*lnu) + q2min: 0 + q2max: 4.85 + value: -0.168 ± 0.068 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 10.689 + value: -0.024 ± 0.043 ± 0.016 + - name: (B0->D*lnu) + q2min: 4.85 + q2max: 10.689 + value: -0.026 ± 0.068 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 4.85 + value: -0.101 ± 0.069 ± 0.025 + - name: (B0->D*lnu) + q2min: 0 + q2max: 10.689 + value: -0.062 ± 0.047 ± 0.017 + - name: (B0->D*lnu) + q2min: 4.85 + q2max: 10.689 + value: -0.019 ± 0.068 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 4.85 + value: -0.055 ± 0.065 ± 0.023 + - name: (B0->D*lnu) + q2min: 0 + q2max: 10.689 + value: -0.035 ± 0.046 ± 0.016 + - name: (B0->D*lnu) + q2min: 4.85 + q2max: 10.689 + value: 0.028 ± 0.067 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 4.85 + value: -0.066 ± 0.065 ± 0.022 + - name: (B0->D*lnu) + q2min: 0 + q2max: 10.689 + value: -0.026 ± 0.046 ± 0.016 + - name: (B0->D*lnu) + q2min: 4.85 + q2max: 10.689 + value: 0.032 ± 0.067 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 4.85 + value: 0.020 ± 0.068 ± 0.024 + - name: (B0->D*lnu) + q2min: 0 + q2max: 10.689 + value: 0.020 ± 0.047 ± 0.017 \ No newline at end of file