Amplitude model with `sympy`#

This section is a follow-up of the previous chapter, Reaction and Models, to formulate the amplitude model for the \(\gamma p \to \eta\pi^0 p\) channel example symbolically. See TR‑033 for a purely numerical tutorial.

The model we want to implement is

\[\begin{split} \begin{array}{rcl} I &=& \left|A^{12} + A^{23} + A^{31}\right|^2 \\ A^{12} &=& \frac{\sum a_m Y_2^m (\Omega_1)}{s_{12}-m^2_{a_2}+im_{a_2} \Gamma_{a_2}} \\ A^{23} &=& \frac{\sum b_m Y_1^m (\Omega_2)}{s_{23}-m^2_{\Delta^+}+im_{\Delta^+} \Gamma_{\Delta^+}} \\ A^{31} &=& \frac{c_0}{s_{31}-m^2_{N^*}+im_{N^*} \Gamma_{N^*}} \,, \end{array} \end{split}\]

where \(1=\eta\), \(2=\pi^0\), and \(3=p\).

Model implementation#

l_max = 2

\(A^{12}\)#

s12, m_a2, gamma_a2, l12 = sp.symbols(r"s_{12} m_{a_2} \Gamma_{a_2} l_{12}")
theta1, phi1 = sp.symbols("theta_1 phi_1")
a = sp.IndexedBase("a")
m = sp.symbols("m", cls=sp.Idx)
A12 = sp.Sum(a[m] * sp.Ynm(l12, m, theta1, phi1), (m, -l12, l12)) / (
    s12 - m_a2**2 + sp.I * m_a2 * gamma_a2
)
A12

\[\displaystyle \frac{\sum_{m=- l_{12}}^{l_{12}} Y_{l_{12}}^{m}\left(\theta_{1},\phi_{1}\right) {a}_{m}}{i \Gamma_{a_2} m_{a_2} - m_{a_2}^{2} + s_{12}}\]

A12_funcs = [
    sp.lambdify(
        [s12, *(a[j] for j in range(-l_max, l_max + 1)), m_a2, gamma_a2, theta1, phi1],
        expr=A12.subs(l12, i).doit().expand(func=True),
    )
    for i in range(l_max + 1)
]

\(A^{23}\)#

s23, m_delta, gamma_delta, l23 = sp.symbols(
    r"s_{23} m_{\Delta^+} \Gamma_{\Delta^+} l_{23}"
)
b = sp.IndexedBase("b")
m = sp.symbols("m", cls=sp.Idx)
theta2, phi2 = sp.symbols("theta_2 phi_2")
A23 = sp.Sum(b[m] * sp.Ynm(l23, m, theta2, phi2), (m, -l23, l23)) / (
    s23 - m_delta**2 + sp.I * m_delta * gamma_delta
)
A23

\[\displaystyle \frac{\sum_{m=- l_{23}}^{l_{23}} Y_{l_{23}}^{m}\left(\theta_{2},\phi_{2}\right) {b}_{m}}{i \Gamma_{\Delta^+} m_{\Delta^+} - \left(m_{\Delta^+}\right)^{2} + s_{23}}\]

A23_funcs = [
    sp.lambdify(
        [
            s23,
            *(b[j] for j in range(-l_max, l_max + 1)),
            m_delta,
            gamma_delta,
            theta2,
            phi2,
        ],
        A23.subs(l23, i).doit().expand(func=True),
    )
    for i in range(l_max + 1)
]

\(A^{31}\)#

c = sp.IndexedBase("c")
s31, m_nstar, gamma_nstar = sp.symbols(r"s_{31} m_{N^*} \Gamma_{N^*}")
theta3, phi3, l31 = sp.symbols("theta_3 phi_3 l_{31}")
A31 = sp.Sum(c[m] * sp.Ynm(l31, m, theta3, phi3), (m, -l31, l31)) / (
    s31 - m_nstar**2 + sp.I * m_nstar * gamma_nstar
)
A31

\[\displaystyle \frac{\sum_{m=- l_{31}}^{l_{31}} Y_{l_{31}}^{m}\left(\theta_{3},\phi_{3}\right) {c}_{m}}{i \Gamma_{N^*} m_{N^*} - \left(m_{N^*}\right)^{2} + s_{31}}\]

A31_funcs = [
    sp.lambdify(
        [
            s31,
            *(c[j] for j in range(-l_max, l_max + 1)),
            m_nstar,
            gamma_nstar,
            theta3,
            phi3,
        ],
        A31.subs(l31, i).doit().expand(func=True),
    )
    for i in range(l_max + 1)
]

\(I = |A|^2 = |A^{12}+A^{23}+A^{31}|^2\)#

intensity_expr = sp.Abs(A12 + A23 + A31) ** 2
intensity_expr

\[\displaystyle \left|{\frac{\sum_{m=- l_{12}}^{l_{12}} Y_{l_{12}}^{m}\left(\theta_{1},\phi_{1}\right) {a}_{m}}{i \Gamma_{a_2} m_{a_2} - m_{a_2}^{2} + s_{12}} + \frac{\sum_{m=- l_{23}}^{l_{23}} Y_{l_{23}}^{m}\left(\theta_{2},\phi_{2}\right) {b}_{m}}{i \Gamma_{\Delta^+} m_{\Delta^+} - \left(m_{\Delta^+}\right)^{2} + s_{23}} + \frac{\sum_{m=- l_{31}}^{l_{31}} Y_{l_{31}}^{m}\left(\theta_{3},\phi_{3}\right) {c}_{m}}{i \Gamma_{N^*} m_{N^*} - \left(m_{N^*}\right)^{2} + s_{31}}}\right|^{2}\]

Phase Space Generation#

Mass for \(p\gamma\) system

E_lab_gamma = 8.5
m_proton = 0.938
m_0 = np.sqrt(2 * E_lab_gamma * m_proton + m_proton**2)
m_eta = 0.548
m_pi = 0.135
m_0

np.float64(4.101931740046389)

rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=m_0,
    final_state_masses={1: m_eta, 2: m_pi, 3: m_proton},
)
phsp_momenta = phsp_generator.generate(500_000, rng)

Kinematic variables#

p1 = FourMomentumSymbol("p1", shape=[])
p2 = FourMomentumSymbol("p2", shape=[])
p3 = FourMomentumSymbol("p3", shape=[])
p12 = ArraySum(p1, p2)
p23 = ArraySum(p2, p3)
p31 = ArraySum(p3, p1)

theta1_expr, phi1_expr = formulate_helicity_angles(p1, p2)
theta2_expr, phi2_expr = formulate_helicity_angles(p2, p3)
theta3_expr, phi3_expr = formulate_helicity_angles(p3, p1)

s12_expr = SquaredInvariantMass(p12)
s23_expr = SquaredInvariantMass(p23)
s31_expr = SquaredInvariantMass(p31)

\[\begin{split}\begin{array}{rcl} \theta_{1} &=& \theta\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{12}}\right|}{E\left({p}_{12}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{12}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{12}\right)\right) p_{1}\right) \\ \theta_{2} &=& \theta\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{23}}\right|}{E\left({p}_{23}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{23}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{23}\right)\right) p_{2}\right) \\ \theta_{3} &=& \theta\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{31}}\right|}{E\left({p}_{31}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{31}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{31}\right)\right) p_{3}\right) \\ \phi_{1} &=& \phi\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{12}}\right|}{E\left({p}_{12}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{12}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{12}\right)\right) p_{1}\right) \\ \phi_{2} &=& \phi\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{23}}\right|}{E\left({p}_{23}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{23}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{23}\right)\right) p_{2}\right) \\ \phi_{3} &=& \phi\left(\boldsymbol{B_z}\left(\frac{\left|\vec{{p}_{31}}\right|}{E\left({p}_{31}\right)}\right) \boldsymbol{R_y}\left(- \theta\left({p}_{31}\right)\right) \boldsymbol{R_z}\left(- \phi\left({p}_{31}\right)\right) p_{3}\right) \\ s_{12} &=& m_{{p}_{12}}^2 \\ s_{23} &=& m_{{p}_{23}}^2 \\ s_{31} &=& m_{{p}_{31}}^2 \\ \end{array}\end{split}\]

helicity_transformer = SympyDataTransformer.from_sympy(
    kinematic_variables, backend="jax"
)

phsp = helicity_transformer(phsp_momenta)
list(phsp)

['theta_1',
 'theta_2',
 'theta_3',
 'phi_1',
 'phi_2',
 'phi_3',
 's_{12}',
 's_{23}',
 's_{31}']

Parameters#

Note

The mass and width of each resonance is customized to make the resonance bands in the Dalitz plot more visible.

Visualization#

Model components#

phi = np.pi / 4
theta = np.pi / 4
s_values = np.linspace(0, 5, num=500)
A12_values = A12_funcs[2](s_values, *a_vals, m_a2_val, gamma_a2_val, theta, phi)
A23_values = A23_funcs[1](s_values, *b_vals, m_delta_val, gamma_delta_val, theta, phi)
A31_values = A31_funcs[0](s_values, *c_vals, m_nstar_val, gamma_nstar_val, theta, phi)

../_images/a9cdcb632d95275c6cc7f2dbe0c6e3ce03bb18e800d142f8336315d5a1b1b932.svg

Unitarity is preserved in each of the subsystems (assuming fixed \(\phi,\theta\)), because we assume there is only one resonance in the subsystem.

../_images/8d2911398f652cf5bf87b8d1559d6368d010f817e14802a5b73b37f5f951c98e.svg

Dalitz Plot#

intensity_funcs = np.array([
    [
        [
            create_parametrized_function(
                expression=intensity_expr.xreplace({l12: i, l23: j, l31: k})
                .doit()
                .expand(func=True),
                parameters=parameters_default,
                backend="jax",
            )
            for i in range(l_max + 1)
        ]
        for j in range(l_max + 1)
    ]
    for k in range(l_max + 1)
])
intensity_funcs.shape

(3, 3, 3)

../_images/2d071cf6b87a17e6d17f287bd85f62daabfbe2d005ee518d67a327d8443f4fc6.png

Projection#

../_images/9d3bdacd2ad0f2ae760b100429da9bc8749716b1de72ab342c36f55ea1c4ec83.svg

Amplitude model with sympy

Contents

Amplitude model with `sympy`#

Model implementation#

\(A^{12}\)#

\(A^{23}\)#

\(A^{31}\)#

\(I = |A|^2 = |A^{12}+A^{23}+A^{31}|^2\)#

Phase Space Generation#

Kinematic variables#

Parameters#

Visualization#

Model components#

Dalitz Plot#

Projection#

Amplitude model with sympy

Contents

Amplitude model with sympy#

Model implementation#

\(A^{12}\)#

\(A^{23}\)#

\(A^{31}\)#

\(I = |A|^2 = |A^{12}+A^{23}+A^{31}|^2\)#

Phase Space Generation#

Kinematic variables#

Parameters#

Visualization#

Model components#

Dalitz Plot#

Projection#

Amplitude model with `sympy`#