From 385ecfccc91fb2b46e96bb3af1f7ef35df6cc6f4 Mon Sep 17 00:00:00 2001 From: jiangyi15 Date: Sun, 14 Jul 2024 18:33:38 +0800 Subject: [PATCH 1/6] feat: dispersion integral models --- tf_pwa/amp/__init__.py | 1 + tf_pwa/amp/dispersion_integral.py | 362 ++++++++++++++++++++++++++++++ tf_pwa/utils.py | 19 +- 3 files changed, 376 insertions(+), 6 deletions(-) create mode 100644 tf_pwa/amp/dispersion_integral.py diff --git a/tf_pwa/amp/__init__.py b/tf_pwa/amp/__init__.py index 645d3ebe..7e5ebc7d 100644 --- a/tf_pwa/amp/__init__.py +++ b/tf_pwa/amp/__init__.py @@ -12,6 +12,7 @@ from .amp import AmplitudeModel, create_amplitude from .base import * from .core import * +from .dispersion_integral import DispersionIntegralParticle from .flatte import ParticleFlatte from .Kmatrix import KmatrixSingleChannelParticle from .kmatrix_simple import KmatrixSimple diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py new file mode 100644 index 00000000..971b23cf --- /dev/null +++ b/tf_pwa/amp/dispersion_integral.py @@ -0,0 +1,362 @@ +import numpy as np +import tensorflow as tf + +from tf_pwa.amp.core import Particle, register_particle +from tf_pwa.data import data_to_numpy + + +def build_integral_tail(fun, x_center, tail, s_th, N=1001, _epsilon=1e-9): + """ + + Integration of the tail parts using tan transfrom. + + .. math:: + \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' = \\int_{\\arctan s_{max}}^{\\frac{\\pi}{2}} \\frac{f(\\tan x)}{\\tan x-s} \\frac{\\mathrm{d} \\tan x}{\\mathrm{d} x} \\mathrm{d} x + + """ + x = np.linspace(np.arctan(tail), np.pi / 2 - _epsilon, N) + tanx = tf.tan(x) + dtanxdx = 1 / tf.cos(x) ** 2 + fx = fun(tanx) / (tanx - x_center[:, None]) * dtanxdx + int_f = tf.reduce_mean(fx, axis=-1) * ( + np.pi / 2 - _epsilon - np.arctan(tail) + ) + return int_f + + +def build_integral(fun, s_range, s_th, N=1001, add_tail=True, method="tf"): + """ + + .. math:: + F(s) = P\\int_{s_{th}}^{\\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + = \\int_{s_{th}}^{s- \\epsilon} \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s + \\epsilon}^{s_{max} } \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + + It require same :math:`\\epsilon` for :math:`s- \\epsilon` and :math:`s- \\epsilon` to get the Cauchy Principal Value. We used bin center to to keep the same :math:`\\epsilon` from left and right bound. + + """ + if method == "scipy": + return build_integral_scipy(fun, s_range, s_th, N=N, add_tail=add_tail) + else: + return build_integral_tf(fun, s_range, s_th, N=N, add_tail=add_tail) + + +def build_integral_scipy( + fun, s_range, s_th, N=1001, add_tail=True, _epsilon=1e-6 +): + """ + + .. math:: + F(s) = P\\int_{s_{th}}^{\\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + = \\int_{s_{th}}^{s- \\epsilon} \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s + \\epsilon}^{s_{max} } \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + + It require same :math:`\\epsilon` for :math:`s- \\epsilon` and :math:`s- \\epsilon` to get the Cauchy Principal Value. We used bin center to to keep the same :math:`\\epsilon` from left and right bound. + + """ + from scipy.integrate import quad + + s_min, s_max = s_range + x = np.linspace(s_min, s_max, N) + + ret = [] + + def f(s): + with tf.device("CPU"): + y = data_to_numpy(fun(s)) + return y + + for xi in x: + if xi < s_th: + y, e = quad(lambda s: f(s) / (s - xi), s_th + _epsilon, np.inf) + else: + y1, e1 = quad(lambda s: f(s) / (s - xi), s_th, xi - _epsilon) + y2, e2 = quad( + lambda s: f(s) / (s - xi), xi + _epsilon, s_max + 0.1 + ) + y3, e2 = quad(lambda s: f(s) / (s - xi), s_max + 0.1, np.inf) + y = y1 + y2 + y3 + ret.append(y) + return x, np.stack(ret) + + +def build_integral_tf(fun, s_range, s_th, N=1001, add_tail=True): + """ + + .. math:: + I(s) = P\\int_{s_{th}}^{\\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + = \\int_{s_{th}}^{s- \\epsilon} \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s + \\epsilon}^{s_{max} } \\frac{f(s')}{s'-s} \\mathrm{d} s' + \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' + + It require same :math:`\\epsilon` for :math:`s- \\epsilon` and :math:`s- \\epsilon` to get the Cauchy Principal Value. We used bin center to to keep the same :math:`\\epsilon` from left and right bound. + + """ + s_min, s_max = s_range + delta = (s_min - s_max) / (N - 1) + x = np.linspace(s_min - delta / 2, s_max + delta / 2, N + 1) + x_center = (x[1:] + x[:-1]) / 2 + int_x = x[x > s_th] + fx = fun(int_x) / (int_x - x_center[:, None]) + int_f = tf.reduce_mean(fx, axis=-1) * (s_max - s_th) + if add_tail: + int_f = int_f + build_integral_tail(fun, x_center, s_max, s_th, N) + return x_center, int_f + + +class LinearInterpFunction: + def __init__(self, x_range, y): + x_min, x_max = x_range + N = y.shape[-1] + self.x_min = x_min + self.x_max = x_max + self.delta = (self.x_max - self.x_min) / (N - 1) + self.N = N + self.y = y + + def __call__(self, x): + diff = (x - self.x_min) / self.delta + idx0 = diff // 1.0 + idx = tf.cast(idx0, tf.int32) + left = tf.gather(self.y, idx) + right = tf.gather(self.y, idx + 1) + k = diff - idx0 + return (1 - k) * left + k * right + + +class DispersionIntegralFunction(LinearInterpFunction): + def __init__(self, fun, s_range, s_th, N=1001, method="tf"): + self.fun = fun + x_center, y = build_integral(fun, s_range, s_th, N, method=method) + super().__init__((x_center[0], x_center[-1]), y) + + +@register_particle("DI") +class DispersionIntegralParticle(Particle): + """ + + "DI" (Dispersion Integral) model is the model used in `PRD78,074023(2008) `_ . In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. + + .. math:: + f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} [Re \\Pi_i(s) - Re\\Pi_i(m_0^2)] - i \\sum_{i} \\rho'_i(s) } + + where :math:`\\rho'_i(s) = g_i^2 \\rho_i(s) F_i^2(s)` is the phase space with barrier factor :math:`F_i^2(s)=\\exp(-\\alpha k_i^2)`. + + The real parts of :math:`\Pi(s)` is defined using the dispersion intergral + + .. math:: + + Re \\Pi_i(s) = \\frac{1}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' = \\lim_{\\epsilon \\rightarrow 0} \\left[ \\int_{s_{th,i}}^{s-\\epsilon} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' +\\int_{s+\\epsilon}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s'\\right] + + The reprodution of the Fig1 in `PRD78,074023(2008) `_ . + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle + >>> plt.clf() + >>> m = np.linspace(0.6, 1.6, 1001) + >>> s = m * m + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> p = DispersionIntegralParticle("A0_980_DI", mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) + >>> p.init_params() + >>> y1 = p.rho_prime(s, *p.mass_list[0]) + >>> scale1 = 1/np.max(y1) + >>> x1 = p.int_f[0](s)/np.pi + >>> p.alpha = 2.5 + >>> p.init_integral() + >>> y2 = p.rho_prime(s, *p.mass_list[0]) + >>> scale2 = 1/np.max(y2) + >>> x2 = p.int_f[0](s)/np.pi + >>> p_ref = DispersionIntegralParticle("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") + >>> p_ref.init_params() + >>> s_ref = np.linspace(0.6**2, 1.6**2, 11) + >>> x_ref = p_ref.int_f[0](s_ref)/np.pi + >>> plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") + >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") + >>> plt.plot(m, y2* scale2, linestyle="--") + >>> plt.plot(m, x2* scale2, linestyle="--") + >>> plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> plt.legend() + + The Argand plot is + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> plt.clf() + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> from tf_pwa.utils import plot_particle_model + >>> plot_particle_model("DI", dict(mass=0.98, mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + + """ + + def __init__( + self, + *args, + mass_range=(0, 5), + mass_list=[[0.493677, 0.493677]], + l_list=None, + int_N=1001, + int_method="tf", + **kwargs + ): + super().__init__(*args, **kwargs) + self.mass_range = mass_range + self.srange = (mass_range[0] ** 2, mass_range[1] ** 2) + self.mass_list = mass_list + self.int_N = int_N + self.int_method = int_method + if l_list is None: + l_list = [0] * len(mass_list) + self.l_list = l_list + + def init_params(self): + super().init_params() + self.alpha = 2.0 + self.gi = [] + for idx, _ in enumerate(self.mass_list): + name = f"g_{idx}" + self.gi.append(self.add_var(name)) + self.init_integral() + + def init_integral(self): + self.int_f = [] + for idx, ((m1, m2), l) in enumerate(zip(self.mass_list, self.l_list)): + fi = lambda s: self.rho_prime(s, m1, m2, l) + int_fi = DispersionIntegralFunction( + fi, + self.srange, + (m1 + m2) ** 2, + N=self.int_N, + method=self.int_method, + ) + self.int_f.append(int_fi) + + def q2_ch(self, s, m1, m2): + return (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / s / 4 + + def rho_prime(self, s, m1, m2, l=0): + q2 = self.q2_ch(s, m1, m2) + q2 = tf.where(s > (m1 + m2) ** 2, q2, tf.zeros_like(q2)) + rho = 2 * tf.sqrt(q2 / s) + F = tf.exp(-self.alpha * q2) + return rho * F**2 + + def __call__(self, m): + s = m * m + m0 = self.get_mass() + s0 = m0 * m0 + gi = tf.stack([var() ** 2 for var in self.gi], axis=-1) + im = tf.stack( + [self.rho_prime(s, m1, m2) for m1, m2 in self.mass_list], axis=-1 + ) + re = tf.stack([f(s) - f(s0) for f in self.int_f], axis=-1) + real = s0 - s - tf.reduce_sum(gi * re, axis=-1) / np.pi + imag = tf.reduce_sum(gi * im, axis=-1) + dom = real**2 + imag**2 + return tf.complex(real / dom, imag / dom) + + def get_amp(self, *args, **kwargs): + return self(args[0]["m"]) + + +@register_particle("DI2") +class DispersionIntegralParticle2(DispersionIntegralParticle): + """ + + Dispersion Integral model. In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. + + .. math:: + f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } + + where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i(s)^2`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. + + The real parts of :math:`\Pi(s)` is defined using the dispersion intergral + + .. math:: + + Re \\Pi_i(s) = \\frac{(s-s_{th,i})}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)(s'-s_{th})} \\mathrm{d} s' + + The shape of Chew-Mandelstam function + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle2 + >>> plt.clf() + >>> m = np.linspace(0.6, 1.6, 1001) + >>> s = m * m + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> p = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=600) + >>> p.init_params() + >>> y1 = p.rho_prime(s, *p.mass_list[0])* (s-M_K**2*4) + >>> scale1 = 1/np.max(y1) + >>> x1 = p.int_f[0](s) * (s-M_K**2*4)/np.pi + >>> p_ref = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") + >>> p_ref.init_params() + >>> s_ref = np.linspace(0.6**2, 1.6**2, 11) + >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref-M_K**2*4)/np.pi + >>> plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") + >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") + >>> plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> plt.legend() + + The Argand plot is + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> plt.clf() + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> from tf_pwa.utils import plot_particle_model + >>> plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=601), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + + """ + + def rho_prime(self, s, m1, m2, l=0): + q2 = self.q2_ch(s, m1, m2) + q2 = tf.where(s > (m1 + m2) ** 2, q2, tf.zeros_like(q2)) + rho = 2 * tf.sqrt(q2 / s) + from tf_pwa.breit_wigner import Bprime_q2 + + rhop = rho * q2**l * Bprime_q2(l, q2, 1 / self.d**2, self.d) + return rhop / (s - (m1 + m2) ** 2) + + def __call__(self, m): + s = m * m + m0 = self.get_mass() + s0 = m0 * m0 + gi = tf.stack([var() ** 2 for var in self.gi], axis=-1) + im = tf.stack( + [ + self.rho_prime(s, m1, m2) * (s - (m1 + m2) ** 2) + for m1, m2 in self.mass_list + ], + axis=-1, + ) + re = tf.stack( + [ + f(s) * (s - (m1 + m2) ** 2) + for f, (m1, m2) in zip(self.int_f, self.mass_list) + ], + axis=-1, + ) + re0 = tf.stack( + [ + f(s0) * (s0 - (m1 + m2) ** 2) + for f, (m1, m2) in zip(self.int_f, self.mass_list) + ], + axis=-1, + ) + real = s0 - s - tf.reduce_sum(gi * (re - re0), axis=-1) / np.pi + imag = tf.reduce_sum(gi * im, axis=-1) + dom = real**2 + imag**2 + return tf.complex(real / dom, imag / dom) + + def get_amp(self, *args, **kwargs): + return self(args[0]["m"]) diff --git a/tf_pwa/utils.py b/tf_pwa/utils.py index 1d951c7e..2165c080 100644 --- a/tf_pwa/utils.py +++ b/tf_pwa/utils.py @@ -341,14 +341,14 @@ def fit_normal(data, weights=None): return np.array([mu, sigma]), np.array([mu_error, sigma_error]) -def create_test_config(model_name, params={}, plot_params={}): +def create_test_config(model_name, params={}, plot_params={}, top_mass=1.0): from tf_pwa.config_loader import ConfigLoader config_dic = { "data": {"dat_order": ["B", "C", "D"]}, "decay": {"A": [["R_BC", "D"]], "R_BC": ["B", "C"]}, "particle": { - "$top": {"A": {"J": 0, "P": -1, "mass": 1.0}}, + "$top": {"A": {"J": 0, "P": -1, "mass": top_mass}}, "$finals": { "B": {"J": 0, "P": -1, "mass": 0.1}, "C": {"J": 0, "P": -1, "mass": 0.1}, @@ -372,14 +372,21 @@ def create_test_config(model_name, params={}, plot_params={}): def plot_particle_model( - model_name, params={}, plot_params={}, axis=None, special_points=None + model_name, + params={}, + plot_params={}, + axis=None, + special_points=None, + mrange=(0.2, 0.9 - 1e-12), ): import matplotlib.pyplot as plt import numpy as np - config = create_test_config(model_name, params, plot_params) + config = create_test_config( + model_name, params, plot_params, top_mass=mrange[1] + 0.1 + 1e-12 + ) f = config.get_particle_function("R_BC") - m = np.linspace(0.2, 0.9 - 1e-12, 2000) + m = np.linspace(mrange[0], mrange[1], 2000) a = f(m).numpy() if special_points is not None: special_points = np.array(special_points) @@ -401,7 +408,7 @@ def plot_particle_model( ax2.scatter(special_points, np.abs(at) ** 2) ax2.set_ylabel("$|A|^2$") ax2.set_ylim((0, None)) - ax2.axvline(x=0.2, linestyle="--") + ax2.axvline(x=mrange[0], linestyle="--") ax2.yaxis.set_label_position("right") ax2.yaxis.tick_right() ax1.plot(np.real(a), m) From ee0993de58b383e46c295461842b4d603692116e Mon Sep 17 00:00:00 2001 From: jiangyi15 Date: Sun, 14 Jul 2024 19:21:22 +0800 Subject: [PATCH 2/6] ci: fixed ci error --- tf_pwa/amp/dispersion_integral.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py index 971b23cf..0fdd7af8 100644 --- a/tf_pwa/amp/dispersion_integral.py +++ b/tf_pwa/amp/dispersion_integral.py @@ -168,7 +168,7 @@ class DispersionIntegralParticle(Particle): >>> x2 = p.int_f[0](s)/np.pi >>> p_ref = DispersionIntegralParticle("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") >>> p_ref.init_params() - >>> s_ref = np.linspace(0.6**2, 1.6**2, 11) + >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) >>> x_ref = p_ref.int_f[0](s_ref)/np.pi >>> plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") @@ -297,7 +297,7 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): >>> x1 = p.int_f[0](s) * (s-M_K**2*4)/np.pi >>> p_ref = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") >>> p_ref.init_params() - >>> s_ref = np.linspace(0.6**2, 1.6**2, 11) + >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref-M_K**2*4)/np.pi >>> plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") From a86db60f93281680ff7232da9219398136442391 Mon Sep 17 00:00:00 2001 From: jiangyi15 Date: Sun, 14 Jul 2024 19:45:50 +0800 Subject: [PATCH 3/6] ci: fix ci error --- tf_pwa/amp/dispersion_integral.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py index 0fdd7af8..328b5c37 100644 --- a/tf_pwa/amp/dispersion_integral.py +++ b/tf_pwa/amp/dispersion_integral.py @@ -170,12 +170,12 @@ class DispersionIntegralParticle(Particle): >>> p_ref.init_params() >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) >>> x_ref = p_ref.int_f[0](s_ref)/np.pi - >>> plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") - >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") - >>> plt.plot(m, y2* scale2, linestyle="--") - >>> plt.plot(m, x2* scale2, linestyle="--") - >>> plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") - >>> plt.legend() + >>> _ = plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") + >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") + >>> _ = plt.plot(m, y2* scale2, linestyle="--") + >>> _ = plt.plot(m, x2* scale2, linestyle="--") + >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> _ = plt.legend() The Argand plot is @@ -187,7 +187,7 @@ class DispersionIntegralParticle(Particle): >>> M_eta = 0.547862 >>> M_pi = 0.1349768 >>> from tf_pwa.utils import plot_particle_model - >>> plot_particle_model("DI", dict(mass=0.98, mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> _ = plot_particle_model("DI", dict(mass=0.98, mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) """ @@ -299,10 +299,10 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): >>> p_ref.init_params() >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref-M_K**2*4)/np.pi - >>> plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") - >>> plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") - >>> plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") - >>> plt.legend() + >>> _ = plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") + >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") + >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> _ = plt.legend() The Argand plot is @@ -314,7 +314,7 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): >>> M_eta = 0.547862 >>> M_pi = 0.1349768 >>> from tf_pwa.utils import plot_particle_model - >>> plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=601), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> _ = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=601), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) """ From 951a67995b8c3d87c25decb9e7280bcb473b78d7 Mon Sep 17 00:00:00 2001 From: jiangyi15 Date: Thu, 5 Sep 2024 20:42:31 +0800 Subject: [PATCH 4/6] refactor: combine diffrent type of dispersion intergral; fixed bugs of sym ls2hel --- tf_pwa/amp/core.py | 10 +- tf_pwa/amp/dispersion_integral.py | 246 +++++++++++++++++++++--------- 2 files changed, 181 insertions(+), 75 deletions(-) diff --git a/tf_pwa/amp/core.py b/tf_pwa/amp/core.py index 43f9a147..ad303c63 100644 --- a/tf_pwa/amp/core.py +++ b/tf_pwa/amp/core.py @@ -890,11 +890,15 @@ def _get_cg_matrix( if out_sym: from sympy.physics.quantum.cg import CG - sint = lambda x: sym.simplify(_spin_int(x * 2)) / 2 + sint = ( + lambda x: sym.simplify(sym.sign(x) * _spin_int(abs(x) * 2)) / 2 + ) sqrt = lambda x: sym.sqrt(sint(x)) - def my_cg_coef(a, b, c, d, e, f): - return CG(*list(map(sint, [a, c, b, d, e, f]))).doit() + def my_cg_coef(j1, j2, m1, m2, j3, m3): + ret = CG(*list(map(sint, [j1, m1, j2, m2, j3, m3]))).doit() + print(j1, j2, m1, m2, j3, m3, ret) + return ret ret = ret.tolist() for i, ls_i in enumerate(ls): diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py index 328b5c37..4f7aa35d 100644 --- a/tf_pwa/amp/dispersion_integral.py +++ b/tf_pwa/amp/dispersion_integral.py @@ -5,25 +5,6 @@ from tf_pwa.data import data_to_numpy -def build_integral_tail(fun, x_center, tail, s_th, N=1001, _epsilon=1e-9): - """ - - Integration of the tail parts using tan transfrom. - - .. math:: - \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' = \\int_{\\arctan s_{max}}^{\\frac{\\pi}{2}} \\frac{f(\\tan x)}{\\tan x-s} \\frac{\\mathrm{d} \\tan x}{\\mathrm{d} x} \\mathrm{d} x - - """ - x = np.linspace(np.arctan(tail), np.pi / 2 - _epsilon, N) - tanx = tf.tan(x) - dtanxdx = 1 / tf.cos(x) ** 2 - fx = fun(tanx) / (tanx - x_center[:, None]) * dtanxdx - int_f = tf.reduce_mean(fx, axis=-1) * ( - np.pi / 2 - _epsilon - np.arctan(tail) - ) - return int_f - - def build_integral(fun, s_range, s_th, N=1001, add_tail=True, method="tf"): """ @@ -91,15 +72,65 @@ def build_integral_tf(fun, s_range, s_th, N=1001, add_tail=True): s_min, s_max = s_range delta = (s_min - s_max) / (N - 1) x = np.linspace(s_min - delta / 2, s_max + delta / 2, N + 1) + shift = (s_th - s_min) % delta + x = x + shift * delta x_center = (x[1:] + x[:-1]) / 2 - int_x = x[x > s_th] + int_x = x[x > s_th + delta / 4] fx = fun(int_x) / (int_x - x_center[:, None]) - int_f = tf.reduce_mean(fx, axis=-1) * (s_max - s_th) + int_f = tf.reduce_mean(fx, axis=-1) * (int_x[-1] - int_x[0] + delta) if add_tail: - int_f = int_f + build_integral_tail(fun, x_center, s_max, s_th, N) + int_f = int_f + build_integral_tail( + fun, x_center, x[-1] + delta / 2, s_th, N + ) return x_center, int_f +def build_integral_tail(fun, x_center, tail, s_th, N=1001, _epsilon=1e-9): + """ + + Integration of the tail parts using tan transfrom. + + .. math:: + \\int_{s_{max}}^{\infty} \\frac{f(s')}{s'-s} \\mathrm{d} s' = \\int_{\\arctan s_{max}}^{\\frac{\\pi}{2}} \\frac{f(\\tan x)}{\\tan x-s} \\frac{\\mathrm{d} \\tan x}{\\mathrm{d} x} \\mathrm{d} x + + """ + x_min, x_max = np.arctan(tail), np.pi / 2 + delta = (x_min - x_max) / N + x = np.linspace(x_min + delta / 2, x_max - delta / 2, N) + tanx = tf.tan(x) + dtanxdx = 1 / tf.cos(x) ** 2 + fx = fun(tanx) / (tanx - x_center[:, None]) * dtanxdx + int_f = tf.reduce_mean(fx, axis=-1) * (np.pi / 2 - np.arctan(tail)) + return int_f + + +def complex_q(s, m1, m2): + q2 = (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / (4 * s) + q = tf.sqrt(tf.complex(q2, tf.zeros_like(q2))) + return q + + +def chew_mandelstam(m, m1, m2): + """ + Chew-Mandelstam function + """ + s = m * m + C = lambda x: tf.complex(x, tf.zeros_like(x)) + m1 = tf.cast(m1, s.dtype) + m2 = tf.cast(m2, s.dtype) + q = complex_q(s, m1, m2) + s1 = m1 * m1 + s2 = m2 * m2 + a = ( + C(2 / m) + * q + * tf.math.log((C(s1 + s2 - s) + C(2 * m) * q) / C(2 * m1 * m2)) + ) + b = (s1 - s2) * (1 / s - 1 / (m1 + m2) ** 2) * tf.math.log(m1 / m2) + ret = a - C(b) + return ret / (16 * np.pi**2) + + class LinearInterpFunction: def __init__(self, x_range, y): x_min, x_max = x_range @@ -131,7 +162,7 @@ def __init__(self, fun, s_range, s_th, N=1001, method="tf"): class DispersionIntegralParticle(Particle): """ - "DI" (Dispersion Integral) model is the model used in `PRD78,074023(2008) `_ . In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. + "DI" (Dispersion Integral) model is the model used in `PRD78,074023(2008) `_ . In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allowed in the integration, unless `dyn_int=True`. .. math:: f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} [Re \\Pi_i(s) - Re\\Pi_i(m_0^2)] - i \\sum_{i} \\rho'_i(s) } @@ -177,7 +208,7 @@ class DispersionIntegralParticle(Particle): >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") >>> _ = plt.legend() - The Argand plot is + The Argand plot .. plot:: @@ -199,6 +230,7 @@ def __init__( l_list=None, int_N=1001, int_method="tf", + dyn_int=False, **kwargs ): super().__init__(*args, **kwargs) @@ -210,6 +242,7 @@ def __init__( if l_list is None: l_list = [0] * len(mass_list) self.l_list = l_list + self.dyn_int = dyn_int def init_params(self): super().init_params() @@ -236,22 +269,34 @@ def init_integral(self): def q2_ch(self, s, m1, m2): return (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / s / 4 + def im_weight(self, s, m1, m2, l=0): + return tf.ones_like(s) + def rho_prime(self, s, m1, m2, l=0): q2 = self.q2_ch(s, m1, m2) q2 = tf.where(s > (m1 + m2) ** 2, q2, tf.zeros_like(q2)) rho = 2 * tf.sqrt(q2 / s) F = tf.exp(-self.alpha * q2) - return rho * F**2 + return rho * F**2 / self.im_weight(s, m1, m2, l) def __call__(self, m): + if self.dyn_int: + self.init_integral() s = m * m m0 = self.get_mass() s0 = m0 * m0 gi = tf.stack([var() ** 2 for var in self.gi], axis=-1) - im = tf.stack( - [self.rho_prime(s, m1, m2) for m1, m2 in self.mass_list], axis=-1 - ) - re = tf.stack([f(s) - f(s0) for f in self.int_f], axis=-1) + ims = [] + res = [] + for (m1, m2), li, f in zip(self.mass_list, self.l_list, self.int_f): + w = self.im_weight(s, m1, m2, li) + w0 = self.im_weight(s0, m1, m2, li) + tmp_i = self.rho_prime(s, m1, m2, li) * w + tmp_r = f(s) * w - f(s0) * w0 + ims.append(tmp_i) + res.append(tmp_r) + im = tf.stack(ims, axis=-1) + re = tf.stack(res, axis=-1) real = s0 - s - tf.reduce_sum(gi * re, axis=-1) / np.pi imag = tf.reduce_sum(gi * im, axis=-1) dom = real**2 + imag**2 @@ -270,41 +315,53 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): .. math:: f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } - where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i(s)^2`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. + where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i^2(s)`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. The real parts of :math:`\Pi(s)` is defined using the dispersion intergral .. math:: - Re \\Pi_i(s) = \\frac{(s-s_{th,i})}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)(s'-s_{th})} \\mathrm{d} s' + Re \\Pi_i(s) = \\frac{\\color{red}(s-s_{th,i})}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)({\\color{red} s'-s_{th,i} })} \\mathrm{d} s' + + .. note:: - The shape of Chew-Mandelstam function + Small `int_N` will have bad precision. + + The shape of :math:`\\Pi(s)` and comparing to Chew-Mandelstam function :math:`\\Sigma(s)` .. plot:: >>> import matplotlib.pyplot as plt - >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle2 + >>> import tensorflow as tf + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle2, chew_mandelstam >>> plt.clf() >>> m = np.linspace(0.6, 1.6, 1001) >>> s = m * m >>> M_K = 0.493677 >>> M_eta = 0.547862 >>> M_pi = 0.1349768 - >>> p = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=600) + >>> p = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) >>> p.init_params() - >>> y1 = p.rho_prime(s, *p.mass_list[0])* (s-M_K**2*4) + >>> p2 = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=1001) + >>> p2.init_params() + >>> y1 = p.rho_prime(s, M_K, M_K)* (s-M_K**2*4) >>> scale1 = 1/np.max(y1) >>> x1 = p.int_f[0](s) * (s-M_K**2*4)/np.pi + >>> x2 = p2.int_f[0](s) * (s-M_K**2*4)/np.pi >>> p_ref = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") >>> p_ref.init_params() >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref-M_K**2*4)/np.pi + >>> x_chew = chew_mandelstam(s, M_K, M_K) * np.pi**2 * 4 + >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s), N=101$") + >>> _ = plt.plot(m, x2* scale1, label="Re $\\\\Pi (s), N=1001$") + >>> _ = plt.plot(m, tf.math.real(x_chew).numpy() * scale1, label="$Re\\\\Sigma(s)$") + >>> _ = plt.plot(m, tf.math.imag(x_chew).numpy() * scale1, label="$Im \\\\Sigma(s)$") >>> _ = plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") - >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") >>> _ = plt.legend() - The Argand plot is + The Argand plot .. plot:: @@ -314,7 +371,10 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): >>> M_eta = 0.547862 >>> M_pi = 0.1349768 >>> from tf_pwa.utils import plot_particle_model - >>> _ = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=601), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> axis = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> _ = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=1001), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05], axis=axis) + >>> _ = axis[3].legend(["N=101", "N=1001"]) + """ @@ -324,39 +384,81 @@ def rho_prime(self, s, m1, m2, l=0): rho = 2 * tf.sqrt(q2 / s) from tf_pwa.breit_wigner import Bprime_q2 - rhop = rho * q2**l * Bprime_q2(l, q2, 1 / self.d**2, self.d) - return rhop / (s - (m1 + m2) ** 2) + rhop = rho * q2**l * Bprime_q2(l, q2, 1 / self.d**2, self.d) ** 2 + return rhop / self.im_weight(s, m1, m2, l) - def __call__(self, m): - s = m * m - m0 = self.get_mass() - s0 = m0 * m0 - gi = tf.stack([var() ** 2 for var in self.gi], axis=-1) - im = tf.stack( - [ - self.rho_prime(s, m1, m2) * (s - (m1 + m2) ** 2) - for m1, m2 in self.mass_list - ], - axis=-1, - ) - re = tf.stack( - [ - f(s) * (s - (m1 + m2) ** 2) - for f, (m1, m2) in zip(self.int_f, self.mass_list) - ], - axis=-1, - ) - re0 = tf.stack( - [ - f(s0) * (s0 - (m1 + m2) ** 2) - for f, (m1, m2) in zip(self.int_f, self.mass_list) - ], - axis=-1, - ) - real = s0 - s - tf.reduce_sum(gi * (re - re0), axis=-1) / np.pi - imag = tf.reduce_sum(gi * im, axis=-1) - dom = real**2 + imag**2 - return tf.complex(real / dom, imag / dom) + def im_weight(self, s, m1, m2, l=0): + return s - (m1 + m2) ** 2 - def get_amp(self, *args, **kwargs): - return self(args[0]["m"]) + +@register_particle("DI3") +class DispersionIntegralParticle3(DispersionIntegralParticle2): + """ + + Dispersion Integral model. In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. + + .. math:: + f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } + + where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i^2(s)`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. + + The real parts of :math:`\Pi(s)` is defined using the dispersion intergral + + .. math:: + + Re \\Pi_i(s) = \\frac{\\color{red}s}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)({\\color{red}s'})} \\mathrm{d} s' + + .. note:: + + Small `int_N` will have bad precision. + + The shape of :math:`\\Pi(s)` + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> import tensorflow as tf + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle3, chew_mandelstam + >>> plt.clf() + >>> m = np.linspace(0.6, 1.6, 1001) + >>> s = m * m + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> p = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) + >>> p.init_params() + >>> p2 = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=1001) + >>> p2.init_params() + >>> y1 = p.rho_prime(s, M_K, M_K)* (s) + >>> scale1 = 1/np.max(y1) + >>> x1 = p.int_f[0](s) * (s)/np.pi + >>> x2 = p2.int_f[0](s) * (s)/np.pi + >>> p_ref = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") + >>> p_ref.init_params() + >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) + >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref)/np.pi + >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s), N=101$") + >>> _ = plt.plot(m, x2* scale1, label="Re $\\\\Pi (s), N=1001$") + >>> _ = plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") + >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> _ = plt.legend() + + The Argand plot + + .. plot:: + + >>> import matplotlib.pyplot as plt + >>> plt.clf() + >>> M_K = 0.493677 + >>> M_eta = 0.547862 + >>> M_pi = 0.1349768 + >>> from tf_pwa.utils import plot_particle_model + >>> axis = plot_particle_model("DI3", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=101, dyn_int=True), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> _ = plot_particle_model("DI3", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=1001), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05], axis=axis) + >>> _ = axis[3].legend(["N=101", "N=1001"]) + + + """ + + def im_weight(self, s, m1, m2, l=0): + return s From 6dd0d0bd36e99ace78cbc5161a80cde2328c5cb7 Mon Sep 17 00:00:00 2001 From: jiangyi15 Date: Thu, 31 Oct 2024 00:51:53 +0800 Subject: [PATCH 5/6] refctor: rewrite DI model --- tf_pwa/amp/dispersion_integral.py | 308 +++++++++++++----------------- tf_pwa/breit_wigner.py | 175 +++++++++++++++++ 2 files changed, 310 insertions(+), 173 deletions(-) diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py index 4f7aa35d..16f5e010 100644 --- a/tf_pwa/amp/dispersion_integral.py +++ b/tf_pwa/amp/dispersion_integral.py @@ -1,7 +1,52 @@ +""" + +The dispersion relation is come from the theory for functions of complex variables. +When a function :math:`f(z)` is analytic in the real axis and vanish when :math:`z \\rightarrow \\infty`, +The function has the property that the integration over the edge of the upper half plane is zero, + +.. math:: + \\int_C \\frac{f(z)}{z -z_0 } \\mathrm{d} z = \\lim_{\\epsilon\\rightarrow 0}\\left[\\int_{-\\infty}^{z_0-\\epsilon}\\frac{f(z)}{z - z_0} \\mathrm{d} z + \\int_{z_0+\\epsilon}^{\infty} \\frac{f(z)}{z -z_0 } \\mathrm{d} z + \\int_{z_0-\\epsilon}^{z_0+\\epsilon}\\frac{f(z)}{z -z_0 } \\mathrm{d} z \\right] = 0 + +.. plot:: + + >>> import matplotlib.pyplot as plt + >>> import numpy as np + >>> _ = plt.plot([1.0], [0.0], marker="o") + >>> phi = np.linspace(0,np.pi) + >>> _ = plt.plot(1.0 + 0.1*np.cos(phi), 0.1*np.sin(phi), c="C1") + >>> _ = plt.plot(-3 * np.cos(phi), 3 *np.sin(phi), c="C2") + >>> _ = plt.plot([-3, 0.9], [0,0], c="C3") + >>> _ = plt.plot([1.1, 3], [0,0], c="C3") + >>> _ = plt.xticks([-3, 0, 1, 3], labels=["$-\\infty$", "0", "z0", "$-\\infty$"]) + >>> _ = plt.yticks([]) + + +The integration suround the pole :math:`z_0` contibute a residue term. + +.. math:: + i \\pi f(z_0) = P\\int_{-\\infty}^{+\\infty}\\frac{f(z)}{z - z_0} \\mathrm{d} z = \\lim_{\\epsilon \\rightarrow 0}\\left[\\int_{-\\infty}^{z_0-\\epsilon}\\frac{f(z)}{z - z_0} \\mathrm{d} z + \\int_{z_0+\\epsilon}^{\infty} \\frac{f(z)}{z -z_0 } \\mathrm{d} z\\right] + + +The physical amplitude have the same property after some subtraction of infinity. + +.. math:: + Re f(s) - Re f(s_0) = \\frac{1}{\\pi}P\\int_{s_{th}}^{\\infty} \left[\\frac{Im f(s')}{s' - s} - \\frac{Im f(s')}{s' - s_0}\\right]\mathrm{d} s' = \\frac{(s-s_0)}{\\pi}P\\int_{s_{th}}^{\\infty} \\frac{Im f(s')}{(s' - s)(s' - s_0)}\mathrm{d} s' + +Sometimes, additional substrction will be used to make sure that the intergration is finity. + +.. math:: + Re f(s) - Re \\left[\\sum_{k=0}^{n-1}\\frac{(s-s_0)^k f^{(k)}(s_0)}{k!}\\right] = \\frac{(s-s_0)^n}{\\pi}P\\int_{s_{th}}^{\\infty} \\frac{Im f(s')}{(s' - s)(s' - s_0)^n}\mathrm{d} s' + +More detials can be found in `Dispersion Relation Integrals `_. + +""" + + import numpy as np import tensorflow as tf from tf_pwa.amp.core import Particle, register_particle +from tf_pwa.breit_wigner import chew_mandelstam, complex_q from tf_pwa.data import data_to_numpy @@ -104,34 +149,9 @@ def build_integral_tail(fun, x_center, tail, s_th, N=1001, _epsilon=1e-9): return int_f -def complex_q(s, m1, m2): - q2 = (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / (4 * s) - q = tf.sqrt(tf.complex(q2, tf.zeros_like(q2))) - return q - - -def chew_mandelstam(m, m1, m2): - """ - Chew-Mandelstam function - """ - s = m * m - C = lambda x: tf.complex(x, tf.zeros_like(x)) - m1 = tf.cast(m1, s.dtype) - m2 = tf.cast(m2, s.dtype) - q = complex_q(s, m1, m2) - s1 = m1 * m1 - s2 = m2 * m2 - a = ( - C(2 / m) - * q - * tf.math.log((C(s1 + s2 - s) + C(2 * m) * q) / C(2 * m1 * m2)) - ) - b = (s1 - s2) * (1 / s - 1 / (m1 + m2) ** 2) * tf.math.log(m1 / m2) - ret = a - C(b) - return ret / (16 * np.pi**2) - - class LinearInterpFunction: + "class for linear interpolation" + def __init__(self, x_range, y): x_min, x_max = x_range N = y.shape[-1] @@ -152,6 +172,8 @@ def __call__(self, x): class DispersionIntegralFunction(LinearInterpFunction): + "class for interpolation of dispersion integral." + def __init__(self, fun, s_range, s_th, N=1001, method="tf"): self.fun = fun x_center, y = build_integral(fun, s_range, s_th, N, method=method) @@ -162,50 +184,58 @@ def __init__(self, fun, s_range, s_th, N=1001, method="tf"): class DispersionIntegralParticle(Particle): """ - "DI" (Dispersion Integral) model is the model used in `PRD78,074023(2008) `_ . In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allowed in the integration, unless `dyn_int=True`. + Dispersion Integral model. In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. .. math:: - f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} [Re \\Pi_i(s) - Re\\Pi_i(m_0^2)] - i \\sum_{i} \\rho'_i(s) } + f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } - where :math:`\\rho'_i(s) = g_i^2 \\rho_i(s) F_i^2(s)` is the phase space with barrier factor :math:`F_i^2(s)=\\exp(-\\alpha k_i^2)`. + where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i^2(s)`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. The real parts of :math:`\Pi(s)` is defined using the dispersion intergral .. math:: - Re \\Pi_i(s) = \\frac{1}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' = \\lim_{\\epsilon \\rightarrow 0} \\left[ \\int_{s_{th,i}}^{s-\\epsilon} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' +\\int_{s+\\epsilon}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s'\\right] + Re \\Pi_i(s) = \\frac{\\color{red}(s-s_{0,i})}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s){\\color{red} (s'-s_{0,i})}} \\mathrm{d} s' - The reprodution of the Fig1 in `PRD78,074023(2008) `_ . + By default, :math:`s_{0,i}=0`, it can be change to other value though option `s0: value`. `value=sth` for :math:`s_{th,i}`. + + .. note:: + + Small `int_N` will have bad precision. + + The shape of :math:`\\Pi(s)` and comparing to Chew-Mandelstam function :math:`\\Sigma(s)` .. plot:: >>> import matplotlib.pyplot as plt - >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle + >>> import tensorflow as tf + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle, chew_mandelstam, complex_q >>> plt.clf() >>> m = np.linspace(0.6, 1.6, 1001) >>> s = m * m >>> M_K = 0.493677 >>> M_eta = 0.547862 >>> M_pi = 0.1349768 - >>> p = DispersionIntegralParticle("A0_980_DI", mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) + >>> p = DispersionIntegralParticle("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) >>> p.init_params() - >>> y1 = p.rho_prime(s, *p.mass_list[0]) - >>> scale1 = 1/np.max(y1) - >>> x1 = p.int_f[0](s)/np.pi - >>> p.alpha = 2.5 - >>> p.init_integral() - >>> y2 = p.rho_prime(s, *p.mass_list[0]) - >>> scale2 = 1/np.max(y2) - >>> x2 = p.int_f[0](s)/np.pi + >>> p2 = DispersionIntegralParticle("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=1001) + >>> p2.init_params() + >>> y1 = p.rho_prime(s, M_K, M_K)* (s) + >>> x1 = p.int_f[0](s) * (s)/np.pi + >>> x2 = p2.int_f[0](s) * (s)/np.pi >>> p_ref = DispersionIntegralParticle("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") >>> p_ref.init_params() >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) - >>> x_ref = p_ref.int_f[0](s_ref)/np.pi - >>> _ = plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") - >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") - >>> _ = plt.plot(m, y2* scale2, linestyle="--") - >>> _ = plt.plot(m, x2* scale2, linestyle="--") - >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") + >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref)/np.pi + >>> x_chew = chew_mandelstam(m, M_K, M_K) + >>> x_qm = 2 * complex_q(s, M_K, M_K).numpy() / m + >>> _ = plt.plot(m, (x1 - np.max(x2)), label="Re $\\\\Pi (s) - \\\\Pi(s_{th}), N=101$") + >>> _ = plt.plot(m, (x2 - np.max(x2)), label="Re $\\\\Pi (s) - \\\\Pi(s_{th}), N=1001$") + >>> _ = plt.plot(m, y1, label="Im $\\\\Pi (s)$") + >>> _ = plt.plot(m, tf.math.real(x_chew).numpy(), label="$Re\\\\Sigma(s)$", ls="--") + >>> _ = plt.plot(m, tf.math.imag(x_chew).numpy(), label="$Im \\\\Sigma(s)$", ls="--") + >>> _ = plt.plot(m, np.real(x_qm), label="2q/m", ls=":") + >>> _ = plt.scatter(np.sqrt(s_ref), (x_ref - np.max(x2)), label="scipy integration") >>> _ = plt.legend() The Argand plot @@ -218,7 +248,10 @@ class DispersionIntegralParticle(Particle): >>> M_eta = 0.547862 >>> M_pi = 0.1349768 >>> from tf_pwa.utils import plot_particle_model - >>> _ = plot_particle_model("DI", dict(mass=0.98, mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> axis = plot_particle_model("DI", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,0], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) + >>> _ = plot_particle_model("DI", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,0], int_N=1001), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05], axis=axis) + >>> _ = axis[3].legend(["N=101", "N=1001"]) + """ @@ -231,6 +264,7 @@ def __init__( int_N=1001, int_method="tf", dyn_int=False, + s0=0.0, **kwargs ): super().__init__(*args, **kwargs) @@ -243,6 +277,9 @@ def __init__( l_list = [0] * len(mass_list) self.l_list = l_list self.dyn_int = dyn_int + if not isinstance(s0, (list, tuple)): + s0 = [s0] * len(mass_list) + self.s0 = s0 def init_params(self): super().init_params() @@ -255,8 +292,10 @@ def init_params(self): def init_integral(self): self.int_f = [] - for idx, ((m1, m2), l) in enumerate(zip(self.mass_list, self.l_list)): - fi = lambda s: self.rho_prime(s, m1, m2, l) + for idx, ((m1, m2), l, sA) in enumerate( + zip(self.mass_list, self.l_list, self.s0) + ): + fi = lambda s: self.rho_prime(s, m1, m2, l, sA) int_fi = DispersionIntegralFunction( fi, self.srange, @@ -269,15 +308,19 @@ def init_integral(self): def q2_ch(self, s, m1, m2): return (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / s / 4 - def im_weight(self, s, m1, m2, l=0): - return tf.ones_like(s) - - def rho_prime(self, s, m1, m2, l=0): + def rho_prime(self, s, m1, m2, l=0, s0=0.0): q2 = self.q2_ch(s, m1, m2) q2 = tf.where(s > (m1 + m2) ** 2, q2, tf.zeros_like(q2)) rho = 2 * tf.sqrt(q2 / s) - F = tf.exp(-self.alpha * q2) - return rho * F**2 / self.im_weight(s, m1, m2, l) + from tf_pwa.breit_wigner import Bprime_q2 + + rhop = rho * q2**l * Bprime_q2(l, q2, 1 / self.d**2, self.d) ** 2 + return rhop / self.im_weight(s, m1, m2, l, s0) + + def im_weight(self, s, m1, m2, l=0, s0=0.0): + if s0 is None or s0 == "sth": + s0 = (m1 + m2) ** 2 + return s - s0 def __call__(self, m): if self.dyn_int: @@ -288,10 +331,12 @@ def __call__(self, m): gi = tf.stack([var() ** 2 for var in self.gi], axis=-1) ims = [] res = [] - for (m1, m2), li, f in zip(self.mass_list, self.l_list, self.int_f): - w = self.im_weight(s, m1, m2, li) - w0 = self.im_weight(s0, m1, m2, li) - tmp_i = self.rho_prime(s, m1, m2, li) * w + for (m1, m2), li, f, sA in zip( + self.mass_list, self.l_list, self.int_f, self.s0 + ): + w = self.im_weight(s, m1, m2, li, sA) + w0 = self.im_weight(s0, m1, m2, li, sA) + tmp_i = self.rho_prime(s, m1, m2, li, sA) * w tmp_r = f(s) * w - f(s0) * w0 ims.append(tmp_i) res.append(tmp_r) @@ -306,58 +351,53 @@ def get_amp(self, *args, **kwargs): return self(args[0]["m"]) -@register_particle("DI2") -class DispersionIntegralParticle2(DispersionIntegralParticle): +@register_particle("DI_a0") +class DispersionIntegralParticleA0(DispersionIntegralParticle): """ - Dispersion Integral model. In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. + "DI_a0" model is the model used in `PRD78,074023(2008) `_ . In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allowed in the integration, unless `dyn_int=True`. .. math:: - f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } + f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} [Re \\Pi_i(s) - Re\\Pi_i(m_0^2)] - i \\sum_{i} \\rho'_i(s) } - where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i^2(s)`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. + where :math:`\\rho'_i(s) = g_i^2 \\rho_i(s) F_i^2(s)` is the phase space with barrier factor :math:`F_i^2(s)=\\exp(-\\alpha k_i^2)`. - The real parts of :math:`\Pi(s)` is defined using the dispersion intergral + The real parts of :math:`\\Pi(s)` is defined using the dispersion intergral .. math:: - Re \\Pi_i(s) = \\frac{\\color{red}(s-s_{th,i})}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)({\\color{red} s'-s_{th,i} })} \\mathrm{d} s' - - .. note:: - - Small `int_N` will have bad precision. + Re \\Pi_i(s) = \\frac{1}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' = \\lim_{\\epsilon \\rightarrow 0} \\left[ \\int_{s_{th,i}}^{s-\\epsilon} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s' +\\int_{s+\\epsilon}^{\\infty} \\frac{\\rho'_i(s')}{s' - s} \\mathrm{d} s'\\right] - The shape of :math:`\\Pi(s)` and comparing to Chew-Mandelstam function :math:`\\Sigma(s)` + The reprodution of the Fig1 in `PRD78,074023(2008) `_ . .. plot:: >>> import matplotlib.pyplot as plt - >>> import tensorflow as tf - >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle2, chew_mandelstam + >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticleA0 >>> plt.clf() >>> m = np.linspace(0.6, 1.6, 1001) >>> s = m * m >>> M_K = 0.493677 >>> M_eta = 0.547862 >>> M_pi = 0.1349768 - >>> p = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) + >>> p = DispersionIntegralParticleA0("A0_980_DI", mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) >>> p.init_params() - >>> p2 = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=1001) - >>> p2.init_params() - >>> y1 = p.rho_prime(s, M_K, M_K)* (s-M_K**2*4) + >>> y1 = p.rho_prime(s, *p.mass_list[0]) >>> scale1 = 1/np.max(y1) - >>> x1 = p.int_f[0](s) * (s-M_K**2*4)/np.pi - >>> x2 = p2.int_f[0](s) * (s-M_K**2*4)/np.pi - >>> p_ref = DispersionIntegralParticle2("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") + >>> x1 = p.int_f[0](s)/np.pi + >>> p.alpha = 2.5 + >>> p.init_integral() + >>> y2 = p.rho_prime(s, *p.mass_list[0]) + >>> scale2 = 1/np.max(y2) + >>> x2 = p.int_f[0](s)/np.pi + >>> p_ref = DispersionIntegralParticleA0("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") >>> p_ref.init_params() >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) - >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref-M_K**2*4)/np.pi - >>> x_chew = chew_mandelstam(s, M_K, M_K) * np.pi**2 * 4 - >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s), N=101$") - >>> _ = plt.plot(m, x2* scale1, label="Re $\\\\Pi (s), N=1001$") - >>> _ = plt.plot(m, tf.math.real(x_chew).numpy() * scale1, label="$Re\\\\Sigma(s)$") - >>> _ = plt.plot(m, tf.math.imag(x_chew).numpy() * scale1, label="$Im \\\\Sigma(s)$") - >>> _ = plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") + >>> x_ref = p_ref.int_f[0](s_ref)/np.pi + >>> _ = plt.plot(m, y1* scale1, label="$\\\\rho'(s)$") + >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s)$") + >>> _ = plt.plot(m, y2* scale2, linestyle="--", label="$\\\\rho'(s),\\\\alpha=2.5$") + >>> _ = plt.plot(m, x2* scale2, linestyle="--", label="Re $\\\\Pi (s),\\\\alpha=2.5$") >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") >>> _ = plt.legend() @@ -371,94 +411,16 @@ class DispersionIntegralParticle2(DispersionIntegralParticle): >>> M_eta = 0.547862 >>> M_pi = 0.1349768 >>> from tf_pwa.utils import plot_particle_model - >>> axis = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) - >>> _ = plot_particle_model("DI2", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=1001), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05], axis=axis) - >>> _ = axis[3].legend(["N=101", "N=1001"]) - + >>> _ = plot_particle_model("DI_a0", dict(mass=0.98, mass_range=(0,2.0), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) """ - def rho_prime(self, s, m1, m2, l=0): + def rho_prime(self, s, m1, m2, l=0, s0=None): q2 = self.q2_ch(s, m1, m2) q2 = tf.where(s > (m1 + m2) ** 2, q2, tf.zeros_like(q2)) rho = 2 * tf.sqrt(q2 / s) - from tf_pwa.breit_wigner import Bprime_q2 - - rhop = rho * q2**l * Bprime_q2(l, q2, 1 / self.d**2, self.d) ** 2 - return rhop / self.im_weight(s, m1, m2, l) - - def im_weight(self, s, m1, m2, l=0): - return s - (m1 + m2) ** 2 - - -@register_particle("DI3") -class DispersionIntegralParticle3(DispersionIntegralParticle2): - """ - - Dispersion Integral model. In the model a linear interpolation is used to avoid integration every times in fitting. No paramters are allow in the integration. - - .. math:: - f(s) = \\frac{1}{m_0^2 - s - \\sum_{i} g_i^2 [Re\\Pi_i(s) -Re\\Pi_i(m_0^2) + i Im \\Pi_i(s)] } - - where :math:`Im \\Pi_i(s)=\\rho_i(s)n_i^2(s)`, :math:`n_i(s)={q}^{l} {B_l'}(q,1/d, d)`. - - The real parts of :math:`\Pi(s)` is defined using the dispersion intergral - - .. math:: - - Re \\Pi_i(s) = \\frac{\\color{red}s}{\\pi} P \\int_{s_{th,i}}^{\\infty} \\frac{Im \\Pi_i(s')}{(s' - s)({\\color{red}s'})} \\mathrm{d} s' - - .. note:: - - Small `int_N` will have bad precision. - - The shape of :math:`\\Pi(s)` - - .. plot:: - - >>> import matplotlib.pyplot as plt - >>> import tensorflow as tf - >>> from tf_pwa.amp.dispersion_integral import DispersionIntegralParticle3, chew_mandelstam - >>> plt.clf() - >>> m = np.linspace(0.6, 1.6, 1001) - >>> s = m * m - >>> M_K = 0.493677 - >>> M_eta = 0.547862 - >>> M_pi = 0.1349768 - >>> p = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=101) - >>> p.init_params() - >>> p2 = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=1001) - >>> p2.init_params() - >>> y1 = p.rho_prime(s, M_K, M_K)* (s) - >>> scale1 = 1/np.max(y1) - >>> x1 = p.int_f[0](s) * (s)/np.pi - >>> x2 = p2.int_f[0](s) * (s)/np.pi - >>> p_ref = DispersionIntegralParticle3("A0_980_DI", mass_range=(0.6,1.6), mass_list=[[M_K, M_K],[M_eta, M_pi]], int_N=11, int_method="scipy") - >>> p_ref.init_params() - >>> s_ref = np.linspace(0.6**2, 1.6**2-1e-6, 11) - >>> x_ref = p_ref.int_f[0](s_ref) * (s_ref)/np.pi - >>> _ = plt.plot(m, x1* scale1, label="Re $\\\\Pi (s), N=101$") - >>> _ = plt.plot(m, x2* scale1, label="Re $\\\\Pi (s), N=1001$") - >>> _ = plt.plot(m, y1* scale1, label="Im $\\\\Pi (s)$") - >>> _ = plt.scatter(np.sqrt(s_ref), x_ref* scale1, label="scipy integration") - >>> _ = plt.legend() - - The Argand plot - - .. plot:: - - >>> import matplotlib.pyplot as plt - >>> plt.clf() - >>> M_K = 0.493677 - >>> M_eta = 0.547862 - >>> M_pi = 0.1349768 - >>> from tf_pwa.utils import plot_particle_model - >>> axis = plot_particle_model("DI3", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=101, dyn_int=True), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05]) - >>> _ = plot_particle_model("DI3", dict(mass=0.98, mass_range=(0.5,1.7), mass_list=[[M_K, M_K],[M_eta, M_pi]], l_list=[0,1], int_N=1001), {"R_BC_g_0": 0.415,"R_BC_g_1": 0.405}, mrange=[0.93, 1.05], axis=axis) - >>> _ = axis[3].legend(["N=101", "N=1001"]) - - - """ + F = tf.exp(-self.alpha * q2) + return rho * F**2 / self.im_weight(s, m1, m2, l, s0) - def im_weight(self, s, m1, m2, l=0): - return s + def im_weight(self, s, m1, m2, l=0, s0=None): + return tf.ones_like(s) diff --git a/tf_pwa/breit_wigner.py b/tf_pwa/breit_wigner.py index e3f54362..fffaf094 100644 --- a/tf_pwa/breit_wigner.py +++ b/tf_pwa/breit_wigner.py @@ -388,3 +388,178 @@ def get_bprime_coeff(l): coeffs = Hjw2.as_dict() ret = [coeffs.get((2 * l - 2 * i,), 0) for i in range(l + 1)] return ret + + +def complex_q(s, m1, m2): + q2 = (s - (m1 + m2) ** 2) * (s - (m1 - m2) ** 2) / (4 * s) + q = tf.sqrt(tf.complex(q2, tf.zeros_like(q2))) + return q + + +def chew_mandelstam(m, m1, m2): + """ + Chew-Mandelstam function in `PDG 2024 Eq 50.44 `_ multiply :math:`16\\pi` factor. + + .. math:: + \\Sigma(m) = \\frac{1}{\\pi}\\left[ + \\frac{2q}{m} \\ln \\left(\\frac{ m_1^2 + m_2^2 - m^2 + 2mq }{ 2 m_1 m_2}\\right) + - (m_1^2 - m_2^2) (\\frac{1}{m^2} - \\frac{1}{(m_1+m_2)^2}) \\ln \\frac{m_1}{m_2} + \\right] + + for :math:`m>(m_1+m_2)` + + .. math:: + Im\\Sigma(m) = \\frac{1}{i}\\frac{1}{\\pi} \\frac{2q}{m} \\ln (-1) = \\frac{2q}{m} + + .. math:: + Re\\Sigma(m) = \\frac{1}{\\pi}\\left[ + \\frac{2q}{m} \\ln \\left( \\frac{ m^2 - m_1^2 - m_2^2 - 2mq }{ 2 m_1 m_2}\\right) + - (m_1^2 - m_2^2) (\\frac{1}{m^2} - \\frac{1}{(m_1+m_2)^2}) \\ln \\frac{m_1}{m_2} + \\right] + + + """ + s = m * m + C = lambda x: tf.complex(x, tf.zeros_like(x)) + m1 = tf.cast(m1, s.dtype) + m2 = tf.cast(m2, s.dtype) + q = complex_q(s, m1, m2) + s1 = m1 * m1 + s2 = m2 * m2 + a = ( + C(2 / m) + * q + * tf.math.log((C(s1 + s2 - s) + C(2 * m) * q) / C(2 * m1 * m2)) + ) + b = (s1 - s2) * (1 / s - 1 / (m1 + m2) ** 2) * tf.math.log(m1 / m2) + ret = a - C(b) + return ret / math.pi + + +def chew_mandelstam_l(m, m1, m2, l): + """ + + TODO. Function from `J.Math.Phys. 25 (1984) 3540 `_ , with some modifies to be same as `chew_mandelstam` function. + + compare with `chew_mandelstam` function. + + + :math:`m_1=0.4, m_2=0.1` + + .. plot:: + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> from tf_pwa.breit_wigner import chew_mandelstam, chew_mandelstam_l + >>> m = np.linspace(0.3, 2.0) + >>> y1 = chew_mandelstam_l(m, 0.4, 0.1, l=0) + >>> y2 = chew_mandelstam(m, 0.4, 0.1) + >>> _ = plt.plot(m, np.real(y1), label="Re $C_0(m)$") + >>> _ = plt.plot(m, np.imag(y1), label="Im $C_0(m)$") + >>> _ = plt.plot(m, np.real(y1)-np.mean(np.real(y1-y2)), label="Re $C_0(m) - \\delta$") + >>> _ = plt.plot(m, np.real(y2), ls="--", label="Re $\\Sigma(m)$") + >>> _ = plt.plot(m, np.imag(y2), ls="--", label="Im $\\Sigma(m)$") + >>> _ = plt.legend() + + :math:`m_1=m_2=0.4` + + .. plot:: + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> from tf_pwa.breit_wigner import chew_mandelstam, chew_mandelstam_l + >>> m = np.linspace(0.3, 2.0) + >>> y1 = chew_mandelstam_l(m, 0.4, 0.4, l=0) + >>> y2 = chew_mandelstam(m, 0.4, 0.4) + >>> _ = plt.plot(m, np.real(y1), label="Re $C_0(m)$") + >>> _ = plt.plot(m, np.imag(y1), label="Im $C_0(m)$") + >>> _ = plt.plot(m, np.real(y1)-np.mean(np.real(y1-y2)), label="Re $C_0(m)- \\delta$") + >>> _ = plt.plot(m, np.real(y2), ls="--", label="Re $\\Sigma(m)$") + >>> _ = plt.plot(m, np.imag(y2), ls="--", label="Im $\\Sigma(m)$") + >>> _ = plt.legend() + + .. plot:: + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> from tf_pwa.breit_wigner import chew_mandelstam, chew_mandelstam_l + >>> m = np.linspace(0.3, 2.0) + >>> y1 = chew_mandelstam_l(m, 0.4, 0.1, l=0) + >>> y2 = chew_mandelstam_l(m, 0.4, 0.1, l=1) + >>> y3 = chew_mandelstam_l(m, 0.4, 0.1, l=2) + >>> _ = plt.plot(m, np.real(y1), label="Re $C_0(m)$") + >>> _ = plt.plot(m, np.imag(y1), label="Im $C_0(m)$") + >>> _ = plt.plot(m, np.real(y2), ls="--", label=" Re $C_1(m)$") + >>> _ = plt.plot(m, np.imag(y2), ls="--", label="Im $C_1(m)$") + >>> _ = plt.plot(m, np.real(y3), ls=":", label="Re $C_2(m)$") + >>> _ = plt.plot(m, np.imag(y3), ls=":", label="Im $C_2(m)$") + >>> _ = plt.legend() + + """ + s = m * m + C = lambda x: tf.complex(x, tf.zeros_like(x)) + same_mass = abs(m1 - m2) < 1e-6 + m1 = tf.cast(m1, s.dtype) + m2 = tf.cast(m2, s.dtype) + s1 = m1 * m1 + s2 = m2 * m2 + a = (m1 + m2) ** 2 + b = (m1 - m2) ** 2 + lam = (s - a) * (s - b) / s / s + + lam_n = [lam**i for i in range(l + 1)] + + ret_1_1 = 0 + for r in range(l + 1): + ret_1_1 = ret_1_1 + lam_n[l - r] / (2 * r + 1) + ret_1 = ( + -( + C(lam) ** (l + 0.5) + * ( + tf.math.log( + -( + ( + (tf.sqrt(C(s - a)) + tf.sqrt(C(s - b))) + / C(2 * tf.sqrt(m1 * m2)) + ) + ** (2) + ) + ) + ) + + C(ret_1_1) + ) + / math.pi + ) + # print(ret_1) + # return ret_1 + if same_mass: # b = 0 + return ret_1 + + mu = (a - b) ** 2 / (16 * a * b) + nu = (a + b) / (2 * tf.sqrt(a * b)) + omega = nu - tf.sqrt(a * b) / s + mu_n = [mu**i for i in range(l + 1)] + + f1, f2 = [1], [1] + for i in range(1, l + 1): + f1.append(f1[-1] * i) + f2.append(f2[-1] * (2 * i - 1)) # TODO: (2n-1)!! for n=0? + Omega_n = [(-2) ** n * f2[n] / f1[n] for n in range(l + 1)] + + ret_2_1 = omega * tf.math.log(m1 / m2) / math.pi + ret_2_2 = Omega_n[0] * lam_n[l] + for i in range(1, l + 1): + ret_2_2 = ret_2_2 + Omega_n[i] * lam_n[l - i] * mu_n[i] + ret_2 = ret_2_1 * ret_2_2 + ret_3_1 = omega * nu / 2 / math.pi + ret_3_2 = 0 + for p in range(1, l + 1): + ret_3_2_1 = 0 + for q in range(l - p + 1): + ret_3_2_1 = ret_3_2_1 + Omega_n[q] * lam_n[l - p - q] * mu_n[q] + ret_3_2_2 = 0 + for r in range(p - 1): + ret_3_2_2 = ret_3_2_2 + Omega_n[r] * mu_n[r] + ret_3_2 = ret_3_2 + ret_3_2_1 * ret_3_2_2 / p + ret_3 = ret_3_1 * ret_3_2 + return ret_1 + C(ret_2 + ret_3) From 81896cdaed9b735abf57a24e880b942cebdcd19c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 30 Oct 2024 16:53:28 +0000 Subject: [PATCH 6/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tf_pwa/amp/dispersion_integral.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tf_pwa/amp/dispersion_integral.py b/tf_pwa/amp/dispersion_integral.py index 16f5e010..29486245 100644 --- a/tf_pwa/amp/dispersion_integral.py +++ b/tf_pwa/amp/dispersion_integral.py @@ -41,7 +41,6 @@ """ - import numpy as np import tensorflow as tf