Polarimeter vector field

Polarimeter vector field#

Amplitude model#

The helicity amplitude for the \(\Lambda_c \to p K \pi\) decay reads as a sum of three partial wave series, incorporating \(\Lambda^{**}\) resonances, \(\Delta^{**}\) resonances, and \(K^{**}\) resonances. The particles are ordered as \(\Lambda_c(\mathbf{0}) \to p(\mathbf{1}) \pi(\mathbf{2}) K(\mathbf{3})\).

(1)#\[\begin{split} \begin{align} \mathcal{A}_{\nu,\lambda}(m_{K\pi},m_{pK}) &= \sum_{\nu',\lambda'} \big[\\ &\qquad d_{\nu,\nu'}^{1/2}(\zeta_{1(1)}^0) \mathcal{A}_{\nu',\lambda'}^{K} d_{\lambda',\lambda}^{1/2}(\zeta_{1(1)}^1) &\mathbf{subsystem\,1\;}(\to 23)\\ &\qquad + d_{\nu,\nu'}^{1/2}(\zeta_{2(1)}^0) \mathcal{A}_{\nu',\lambda'}^{\Lambda} d_{\lambda',\lambda}^{1/2}(\zeta_{2(1)}^1) &\mathbf{subsystem\,2\;}(\to 31)\\ &\qquad + d_{\nu,\nu'}^{1/2}(\zeta_{3(1)}^0) \mathcal{A}_{\nu',\lambda'}^{\Delta} d_{\lambda',\lambda}^{1/2}(\zeta_{3(1)}^1)\big]\,. &\mathbf{subsystem\,3\;}(\to 12)\\ \end{align} \end{split}\]

where \(\zeta^{i}_{j(k)}\) is the Wigner rotation for particle \(i,k\) and chain \(j\). The number in brackets indicates the overall definition of the helicity states \(|1/2,\nu\rangle\) and \(|1/2,\lambda\rangle\) for \(\Lambda_c\) and proton, respectively.

We use the particle-2 convention for the helicity couplings, which leads to the phase factor in the transitions \(\Lambda_c\to K^{**}p\), and \(\Lambda\to K p\):

\[\begin{split} \begin{align} \mathcal{A}^{K}_{\nu,\lambda} &= \sum_{j,\tau} \delta_{\nu,\tau - \lambda}\mathcal{H}^{\Lambda_c \to K^{**} p}_{\tau,\lambda} (-)^{1/2 - \lambda} \,d^{j}_{\lambda,0} (\theta_{23}) \, \mathcal{H}^{K^{**} \to \pi K}_{0,0}\\ % \mathcal{A}^{\Lambda}_{\nu,\lambda} &= \sum_{j,\tau} \delta_{\nu,\tau} \mathcal{H}^{\Lambda_c \to \Lambda^{**} \pi}_{\tau,0} d^{j}_{\tau,-\lambda} (\theta_{31}) \mathcal{H}^{\Lambda^{**} \to K p}_{0,\lambda} (-)^{1/2-\lambda} \\ % \mathcal{A}^{\Delta}_{\nu,\lambda} &= \sum_{j,\tau} \delta_{\nu,\tau} \mathcal{H}^{\Lambda_c \to \Delta^{**} K}_{\tau,0} d^{j}_{\tau,\lambda}(\theta_{12}) \mathcal{H}^{\Delta^{**} \to p\pi}_{\lambda,0} \,. \end{align} \end{split}\]

The helicity couplings in the particle-2 convention obey simple properties with respect to the parity transformation:

\[ \begin{align} \mathcal{H}^{A\to BC}_{-\lambda,-\lambda'} = P_A P_B P_C (-)^{j_A-j_B-j_C} \mathcal{H}^{A\to BC}_{\lambda,\lambda'} \end{align} \]

It reduced amount of the couplings in the strong decay of isobars. Moreover the magnitude of the couplings cannot be determined separately, therefore, it is set to 1:

\[\begin{split} \begin{align} \mathcal{H}^{\Lambda^{**} \to K p}_{0,1/2} &= 1\,, & \mathcal{H}^{\Delta^{**} \to p\pi}_{1/2,0} &= 1\,,& \mathcal{H}^{K^{**} \to \pi K}_{0,0} &= 1, \\ \mathcal{H}^{\Lambda^{**} \to K p}_{0,-1/2} &= -P_\Lambda (-)^{j-1/2}\,, & \mathcal{H}^{\Delta^{**} \to p\pi}_{-1/2,0} &= -P_\Delta (-)^{j-1/2}\,, && \end{align} \end{split}\]

The helicity couplings for the \(\Lambda_c^+\) decay are fit parameters. There are four of them for the \(K^{**}\) chain, and two for both the \(\Delta^{**}\) and \(\Lambda^{**}\) chains.

Resonances and LS-scheme#

Resonance choices and their \(LS\)-couplings are defined as follows:

resonance_names = {
    1: ["K*(892)0"],
    2: ["Lambda(1520)", "Lambda(1670)"],
    3: ["Delta(1232)++"],
}

Compute LS-couplings Hide code cell source

@frozen
class Resonance:
    particle: Particle
    l_R: int
    l_Λc: int

    @staticmethod
    def generate_ls(particle: Particle, chain_id: int) -> Resonance:
        LS_prod = generate_ls(Λc, particle, siblings[chain_id], strong=False)
        LS_prod = [L for L, S in LS_prod]
        LS_dec = generate_ls(particle, *decay_products[chain_id])
        LS_dec = [L for L, S in LS_dec]
        return Resonance(particle, l_R=min(LS_dec), l_Λc=min(LS_prod))


def generate_ls(
    parent: Particle,
    child1: Particle,
    child2: Particle,
    strong: bool = True,
    max_L: int = 3,
):
    s1 = child1.spin
    s2 = child2.spin
    s_values = arange(abs(s1 - s2), s1 + s2)
    LS_values = set()
    for S in s_values:
        for L in arange(0, max_L):
            if not abs(L - S) <= parent.spin <= L + S:
                continue
            η0, η1, η2 = [
                int(parent.parity),
                int(child1.parity),
                int(child2.parity),
            ]
            if strong and η0 != η1 * η2 * (-1) ** L:
                continue
            LS_values.add((L, S))
    return sorted(LS_values)


def arange(x1, x2):
    spin_range = np.arange(float(x1), +float(x2) + 0.5)
    return list(map(sp.Rational, spin_range))


resonance_particles = {
    chain_id: [PDG[name] for name in names]
    for chain_id, names in resonance_names.items()
}
ls_resonances = {
    chain_id: [Resonance.generate_ls(particle, chain_id) for particle in particles]
    for chain_id, particles in resonance_particles.items()
}


def jp(particle: Particle):
    p = "+" if particle.parity > 0 else "-"
    j = sp.Rational(particle.spin)
    return Rf"\({j}^{p}\)"


def create_html_table_row(*items, typ="td"):
    items = (f"<{typ}>{i}</{typ}>" for i in items)
    return "<tr>" + "".join(items) + "</tr>\n"


column_names = [
    "resonance",
    R"\(j^P\)",
    R"\(m\) (MeV)",
    R"\(\Gamma_0\) (MeV)",
    R"\(l_R\)",
    R"\(l_{\Lambda_c}^\mathrm{min}\)",
]
src = "<table>\n"
src += create_html_table_row(*column_names, typ="th")
for chain_id, resonance_list in ls_resonances.items():
    child1, child2 = decay_products[chain_id]
    for resonance in resonance_list:
        src += create_html_table_row(
            Rf"\({resonance.particle.latex} \to {child1.latex} {child2.latex}\)",
            jp(resonance.particle),
            int(1e3 * resonance.particle.mass),
            int(1e3 * resonance.particle.width),
            resonance.l_R,
            resonance.l_Λc,
        )
src += "</table>\n"
HTML(src)

resonance	\(j^P\)	\(m\) (MeV)	\(\Gamma_0\) (MeV)	\(l_R\)	\(l_{\Lambda_c}^\mathrm{min}\)
\(K^{*}(892)^{0} \to \pi^{+} K^{-}\)	\(1^-\)	895	47	1	0
\(\Lambda(1520) \to p K^{-}\)	\(3/2^-\)	1519	16	2	1
\(\Lambda(1670) \to p K^{-}\)	\(1/2^-\)	1674	30	0	0
\(\Delta(1232)^{++} \to p \pi^{+}\)	\(3/2^+\)	1232	117	1	1

Aligned amplitude#

\[\displaystyle \sum_{\lambda^{\prime}=-1/2}^{1/2} \sum_{\nu^{\prime}=-1/2}^{1/2}{A^{K}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{1(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{1(1)}\right) + A^{\Delta}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{3(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{3(1)}\right) + A^{\Lambda}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{2(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{2(1)}\right)}\]

Dynamics#

The lineshape function is factored out of the \(\Lambda_c\) helicity coupling:

\[ \begin{align} \mathcal{H}^{\Lambda_c \to R x}_{\lambda,\lambda'} = \hat{\mathcal{H}}^{\Lambda_c \to R x}_{\lambda,\lambda'}\,\mathcal{R}(s)\,. \end{align} \]

The relativistic Breit-Wigner parametrization reads:

\[ \begin{align} \mathcal{R}(s) = \left(\frac{q}{q_0}\right)^{l_{\Lambda_c}^\text{min}} \frac{F_{l_{\Lambda_c}^\text{min}}(q R)}{F_{l_{\Lambda_c}^\text{min}}(q_0 R)}\, \frac{1}{m^2-s-im\Gamma(s)} \left(\frac{p}{p_0}\right)^{l_R} \frac{F_{l_R}(pR)}{F_{l_R}(p_0R)}, \end{align} \]

with energy-dependent width given by

\[ \begin{align} \Gamma(s) = \Gamma_0 \left(\frac{p}{p_0}\right)^{2l_R+1} \frac{m}{\sqrt{s}} \, \frac{F_{l_R}^2(pR)}{F_{l_R}^2(p_0R)}\,, \end{align} \]

The form-factor \(F\) is the Blatt-Weisskopf factor with the length factor \(R=5\,\)GeV\(^{-1}\):

\[ \begin{align} F_0(pR) &= 1\,,& F_1(pR) &= \sqrt{\frac{1}{1+(pR)^2}}\,,& F_2(pR) &= \sqrt{\frac{1}{9+3(pR)^2+(pR)^4}}\,. \end{align} \]

The break-up momenta is calculated for every decay chain separately. Using the notations \(0->R(\to ij) k\), one writes:

\[ \begin{align} p &= \lambda^{1/2}(s,m_i^2,m_j^2)/(2\sqrt{s})\,, & q &= \lambda^{1/2}(s,m_0^2,m_k^2)/(2m_0)\,. \end{align} \]

The momenta with subindex zero are computed for nominal mass of the resonance, \(s=m^2\). The three-argument Källén function reads:

(2)#\[ \begin{align} \lambda(x,y,z) = x^2+y^2+z^2 - 2xy-2yz-2zx\,. \end{align} \]

Formulation with SymPy

\[\begin{split}\displaystyle \begin{align*} F_{L}\left(z\right) = & \begin{cases} 1 & \text{for}\: L = 0 \\\frac{1}{\sqrt{z^{2} + 1}} & \text{for}\: L = 1 \\\frac{1}{\sqrt{z^{4} + 3 z^{2} + 9}} & \text{for}\: L = 2 \end{cases} \end{align*}\end{split}\]

\[\displaystyle \begin{eqnarray} \lambda\left(x, y, z\right) & = & x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \end{eqnarray}\]

\[\displaystyle \begin{eqnarray} p_{s} & = & \frac{\sqrt{\lambda\left(s, m_{i}^{2}, m_{j}^{2}\right)}}{2 \sqrt{s}} \end{eqnarray}\]

\[\displaystyle \begin{eqnarray} q_{s} & = & \frac{\sqrt{\lambda\left(s, m_{0}^{2}, m_{k}^{2}\right)}}{2 m_{0}} \end{eqnarray}\]

\[\displaystyle \begin{eqnarray} \Gamma\left(s\right) & = & Γ_{0} \frac{m}{\sqrt{s}} \frac{F_{l_{R}}\left(R p_{s}\right)^{2}}{F_{l_{R}}\left(R p_{m^{2}}\right)^{2}} \left(\frac{p_{s}}{p_{m^{2}}}\right)^{2 l_{R} + 1} \end{eqnarray}\]

\[\displaystyle \begin{eqnarray} \mathcal{R}\left(s\right) & = & \frac{\frac{F_{l_{R}}\left(R p_{s}\right)}{F_{l_{R}}\left(R p_{m^{2}}\right)} \frac{F_{l_{\Lambda_c}}\left(R q_{s}\right)}{F_{l_{\Lambda_c}}\left(R q_{m^{2}}\right)} \left(\frac{p_{s}}{p_{m^{2}}}\right)^{l_{R}} \left(\frac{q_{s}}{q_{m^{2}}}\right)^{l_{\Lambda_c}}}{m^{2} - i m \Gamma\left(s\right) - s} \end{eqnarray}\]

Decay chain amplitudes#

\[\displaystyle \sum_{\tau=-1}^{1}{\left(-1\right)^{\frac{1}{2} - \lambda} \delta_{\nu, - \lambda + \tau} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, \tau, \lambda} \mathcal{R}\left(\sigma_{1}\right) d^{1}_{\tau,0}\left(\theta_{23}\right)}\]

\[\displaystyle \sum_{\tau=-3/2}^{3/2}{\left(-1\right)^{\frac{1}{2} - \lambda} \delta_{\nu \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, \lambda} \mathcal{H}^\mathrm{production}_{\Lambda(1520), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{3}{2}}_{\tau,- \lambda}\left(\theta_{31}\right)} + \sum_{\tau=-1/2}^{1/2}{\left(-1\right)^{\frac{1}{2} - \lambda} \delta_{\nu \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, \lambda} \mathcal{H}^\mathrm{production}_{\Lambda(1670), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{1}{2}}_{\tau,- \lambda}\left(\theta_{31}\right)}\]

\[\displaystyle \sum_{\tau=-3/2}^{3/2}{\delta_{\nu \tau} \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, \lambda, 0} \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \tau, 0} \mathcal{R}\left(\sigma_{3}\right) d^{\frac{3}{2}}_{\tau,\lambda}\left(\theta_{12}\right)}\]

Angle definitions#

Angles with repeated lower indices are trivial. The other angles are computed from the invariants and these are the angles with sign

\[\begin{split} \begin{align} \zeta_{1(1)}^{0} &= \hat{\theta}_{1(1)}^{0} = 0\,, & \zeta_{1(1)}^{1} &= 0\,,\\ \zeta_{2(1)}^0 &=\hat{\theta}_{2(1)} =-\hat{\theta}_{1(2)}\,,\\ \zeta_{3(1)}^0 &= \hat{\theta}_{3(1)}\,, &\zeta_{3(1)}^1 &= -\zeta_{1(3)}^1\,, \end{align} \end{split}\]

The expressions for the cosine of the positive (anticlockwise) angles, \(\theta_{12}, \theta_{23}, \theta_{13}\) and \(\hat\theta_{1(2)}, \hat\theta_{3(1)}, \zeta^1_{1(3)}\) can be expressed in terms of Mandelstam variables \(\sigma_1, \sigma_2, \sigma_3\) using [Mikhasenko et al., 2020], Appendix A:

\[\begin{split}\displaystyle \begin{array}{rcl} \theta_{12} & = & \operatorname{acos}{\left(\frac{2 \sigma_{3} \left(- m_{K}^{2} - m_{p}^{2} + \sigma_{2}\right) - \left(- m_{K}^{2} + m_{\Lambda_c}^{2} - \sigma_{3}\right) \left(m_{p}^{2} - m_{\pi}^{2} + \sigma_{3}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{K}^{2}, \sigma_{3}\right)} \sqrt{\lambda\left(\sigma_{3}, m_{p}^{2}, m_{\pi}^{2}\right)}} \right)} \\ \theta_{23} & = & \operatorname{acos}{\left(\frac{2 \sigma_{1} \left(- m_{p}^{2} - m_{\pi}^{2} + \sigma_{3}\right) - \left(- m_{K}^{2} + m_{\pi}^{2} + \sigma_{1}\right) \left(- m_{p}^{2} + m_{\Lambda_c}^{2} - \sigma_{1}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{p}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{1}, m_{\pi}^{2}, m_{K}^{2}\right)}} \right)} \\ \theta_{31} & = & \operatorname{acos}{\left(\frac{2 \sigma_{2} \left(- m_{K}^{2} - m_{\pi}^{2} + \sigma_{1}\right) - \left(m_{K}^{2} - m_{p}^{2} + \sigma_{2}\right) \left(- m_{\pi}^{2} + m_{\Lambda_c}^{2} - \sigma_{2}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{\pi}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{K}^{2}, m_{p}^{2}\right)}} \right)} \\ \zeta^0_{1(1)} & = & 0 \\ \zeta^0_{2(1)} & = & - \operatorname{acos}{\left(\frac{- 2 m_{\Lambda_c}^{2} \left(- m_{p}^{2} - m_{\pi}^{2} + \sigma_{3}\right) + \left(m_{p}^{2} + m_{\Lambda_c}^{2} - \sigma_{1}\right) \left(m_{\pi}^{2} + m_{\Lambda_c}^{2} - \sigma_{2}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{\pi}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(m_{\Lambda_c}^{2}, \sigma_{1}, m_{p}^{2}\right)}} \right)} \\ \zeta^0_{3(1)} & = & \operatorname{acos}{\left(\frac{- 2 m_{\Lambda_c}^{2} \left(- m_{K}^{2} - m_{p}^{2} + \sigma_{2}\right) + \left(m_{K}^{2} + m_{\Lambda_c}^{2} - \sigma_{3}\right) \left(m_{p}^{2} + m_{\Lambda_c}^{2} - \sigma_{1}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{p}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(m_{\Lambda_c}^{2}, \sigma_{3}, m_{K}^{2}\right)}} \right)} \\ \zeta^1_{1(1)} & = & 0 \\ \zeta^1_{2(1)} & = & \operatorname{acos}{\left(\frac{2 m_{p}^{2} \left(- m_{K}^{2} - m_{\Lambda_c}^{2} + \sigma_{3}\right) + \left(- m_{K}^{2} - m_{p}^{2} + \sigma_{2}\right) \left(m_{p}^{2} + m_{\Lambda_c}^{2} - \sigma_{1}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{p}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{p}^{2}, m_{K}^{2}\right)}} \right)} \\ \zeta^1_{3(1)} & = & - \operatorname{acos}{\left(\frac{2 m_{p}^{2} \left(- m_{\pi}^{2} - m_{\Lambda_c}^{2} + \sigma_{2}\right) + \left(- m_{p}^{2} - m_{\pi}^{2} + \sigma_{3}\right) \left(m_{p}^{2} + m_{\Lambda_c}^{2} - \sigma_{1}\right)}{\sqrt{\lambda\left(m_{\Lambda_c}^{2}, m_{p}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{3}, m_{p}^{2}, m_{\pi}^{2}\right)}} \right)} \\ \end{array}\end{split}\]

where \(m_0\) is the mass of the initial state \(\Lambda_c\) and \(m_1, m_2, m_3\) are the masses of \(p, \pi, K\), respectively:

\[\begin{split}\displaystyle \begin{array}{rcl} m_{\Lambda_c} & = & 2.28646 \\ m_{p} & = & 0.938272081 \\ m_{\pi} & = & 0.13957039 \\ m_{K} & = & 0.493677 \\ \end{array}\end{split}\]

Helicity coupling values#

\[\begin{split}\displaystyle \begin{array}{rcl} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} & = & 1 \\ \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, - \frac{1}{2}} & = & -1 \\ \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, - \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, \frac{1}{2}, 0} & = & 1 \\ \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, - \frac{1}{2}, 0} & = & 1 \\ \end{array}\end{split}\]

\[\begin{split}\displaystyle \begin{array}{rcl} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, 0, - \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, -1, - \frac{1}{2}} & = & 1.0 - 1.0 i \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, 1, \frac{1}{2}} & = & -3.0 - 3.0 i \\ \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, 0, \frac{1}{2}} & = & -1.0 - 4.0 i \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{0}, 0, - \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{0}, -1, - \frac{1}{2}} & = & 1.0 - 1.0 i \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{0}, 1, \frac{1}{2}} & = & -3.0 - 3.0 i \\ \mathcal{H}^\mathrm{production}_{K_{0}^{*}(1430)^{0}, 0, \frac{1}{2}} & = & -1.0 - 4.0 i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{0}, 0, - \frac{1}{2}} & = & 1 \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{0}, -1, - \frac{1}{2}} & = & 1.0 - 1.0 i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{0}, 1, \frac{1}{2}} & = & -3.0 - 3.0 i \\ \mathcal{H}^\mathrm{production}_{K_{2}^{*}(1430)^{0}, 0, \frac{1}{2}} & = & -1.0 - 4.0 i \\ \mathcal{H}^\mathrm{production}_{\Lambda(1520), \frac{1}{2}, 0} & = & 1.5 \\ \mathcal{H}^\mathrm{production}_{\Lambda(1520), - \frac{1}{2}, 0} & = & 0.3 \\ \mathcal{H}^\mathrm{production}_{\Lambda(1670), \frac{1}{2}, 0} & = & -0.5 + 1.0 i \\ \mathcal{H}^\mathrm{production}_{\Lambda(1670), - \frac{1}{2}, 0} & = & -0.3 - 0.1 i \\ \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \frac{1}{2}, 0} & = & -13.0 + 5.0 i \\ \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, - \frac{1}{2}, 0} & = & -7.0 + 3.0 i \\ \end{array}\end{split}\]

Intensity expression#

Incoherent sum of the amplitudes defined by Aligned amplitude:

\[\displaystyle \sum_{\lambda=-1/2}^{1/2} \sum_{\nu=-1/2}^{1/2}{\left|{\sum_{\lambda^{\prime}=-1/2}^{1/2} \sum_{\nu^{\prime}=-1/2}^{1/2}{A^{K}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{1(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{1(1)}\right) + A^{\Delta}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{3(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{3(1)}\right) + A^{\Lambda}_{\nu^{\prime}, \lambda^{\prime}} d^{\frac{1}{2}}_{\lambda^{\prime},\lambda}\left(\zeta^1_{2(1)}\right) d^{\frac{1}{2}}_{\nu,\nu^{\prime}}\left(\zeta^0_{2(1)}\right)}}\right|^{2}}\]

Remaining free_symbols are indeed the specific amplitudes as defined by Decay chain amplitudes:

The specific amplitudes from Decay chain amplitudes need to be formulated for each value of \(\nu, \lambda\), so that they can be substituted in the top expression:

\[\begin{split}\displaystyle \begin{array}{rcl} A^{K}_{- \frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-1}^{1}{- \delta_{- \frac{1}{2}, \tau + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, \tau, - \frac{1}{2}} \mathcal{R}\left(\sigma_{1}\right) d^{1}_{\tau,0}\left(\theta_{23}\right)} \\ A^{K}_{- \frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-1}^{1}{\delta_{- \frac{1}{2}, \tau - \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, \tau, \frac{1}{2}} \mathcal{R}\left(\sigma_{1}\right) d^{1}_{\tau,0}\left(\theta_{23}\right)} \\ A^{K}_{\frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-1}^{1}{- \delta_{\frac{1}{2}, \tau + \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, \tau, - \frac{1}{2}} \mathcal{R}\left(\sigma_{1}\right) d^{1}_{\tau,0}\left(\theta_{23}\right)} \\ A^{K}_{\frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-1}^{1}{\delta_{\frac{1}{2}, \tau - \frac{1}{2}} \mathcal{H}^\mathrm{decay}_{K^{*}(892)^{0}, 0, 0} \mathcal{H}^\mathrm{production}_{K^{*}(892)^{0}, \tau, \frac{1}{2}} \mathcal{R}\left(\sigma_{1}\right) d^{1}_{\tau,0}\left(\theta_{23}\right)} \\ A^{\Lambda}_{- \frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{- \delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, - \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1520), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{3}{2}}_{\tau,\frac{1}{2}}\left(\theta_{31}\right)} + \sum_{\tau=-1/2}^{1/2}{- \delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, - \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1670), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{1}{2}}_{\tau,\frac{1}{2}}\left(\theta_{31}\right)} \\ A^{\Lambda}_{- \frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1520), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{3}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{31}\right)} + \sum_{\tau=-1/2}^{1/2}{\delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1670), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{1}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{31}\right)} \\ A^{\Lambda}_{\frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{- \delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, - \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1520), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{3}{2}}_{\tau,\frac{1}{2}}\left(\theta_{31}\right)} + \sum_{\tau=-1/2}^{1/2}{- \delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, - \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1670), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{1}{2}}_{\tau,\frac{1}{2}}\left(\theta_{31}\right)} \\ A^{\Lambda}_{\frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1520), 0, \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1520), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{3}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{31}\right)} + \sum_{\tau=-1/2}^{1/2}{\delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Lambda(1670), 0, \frac{1}{2}} \mathcal{H}^\mathrm{production}_{\Lambda(1670), \tau, 0} \mathcal{R}\left(\sigma_{2}\right) d^{\frac{1}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{31}\right)} \\ A^{\Delta}_{- \frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, - \frac{1}{2}, 0} \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \tau, 0} \mathcal{R}\left(\sigma_{3}\right) d^{\frac{3}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{12}\right)} \\ A^{\Delta}_{- \frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{- \frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, \frac{1}{2}, 0} \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \tau, 0} \mathcal{R}\left(\sigma_{3}\right) d^{\frac{3}{2}}_{\tau,\frac{1}{2}}\left(\theta_{12}\right)} \\ A^{\Delta}_{\frac{1}{2}, - \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, - \frac{1}{2}, 0} \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \tau, 0} \mathcal{R}\left(\sigma_{3}\right) d^{\frac{3}{2}}_{\tau,- \frac{1}{2}}\left(\theta_{12}\right)} \\ A^{\Delta}_{\frac{1}{2}, \frac{1}{2}} & = & \sum_{\tau=-3/2}^{3/2}{\delta_{\frac{1}{2} \tau} \mathcal{H}^\mathrm{decay}_{\Delta(1232)^{++}, \frac{1}{2}, 0} \mathcal{H}^\mathrm{production}_{\Delta(1232)^{++}, \tau, 0} \mathcal{R}\left(\sigma_{3}\right) d^{\frac{3}{2}}_{\tau,\frac{1}{2}}\left(\theta_{12}\right)} \\ \end{array}\end{split}\]

Polarization sensitivity#

We introduce the polarimeter vector field (polarization sensitivity) of the \(\Lambda_c\) decay. It is defined by three quantities \((\alpha_x,\alpha_y,\alpha_z)\) forming a three-dimensional vector \(\vec\alpha\) dependent on just two decay variables, \(\sigma_1=m_{K\pi}^2\), and \(\sigma_2=m_{pK}^2\).

The polarimeter vector field is computed by averaging the Pauli matrices \(\vec\sigma\) contracted with the \(\Lambda_c^+\) helicity indices given the transition amplitude.

(3)#\[ \begin{align} \vec\alpha(m_{K\pi},m_{pK}) = \sum_{\lambda,\nu,\nu'} A^{*}_{\nu,\lambda}\vec\sigma_{\nu,\nu'} A_{\nu',\lambda} \,\big / \sum_{\lambda,\nu} \left|A_{\nu,\lambda}\right|^2 \end{align} \]

The quantities \(\vec\alpha(m_{K\pi},m_{pK})\) give the model-independent representation of the \(\Lambda_c^+\) decay process. It can be used to study \(\Lambda_c^+\) production polarization using

(4)#\[ \begin{align} I(\alpha,\beta,\gamma,m_{K\pi},m_{pK}) = I_0(m_{K\pi},m_{pK})\, \left(1 + \sum_{i,j} P_i R_{ij}(\alpha,\beta,\gamma) \alpha_j(m_{K\pi},m_{pK}) \right)\,, \end{align} \]

where \(R_{ij}(\alpha,\beta,\gamma)\) is a three-dimensional rotation matrix:

\[ \begin{align} R(\alpha,\beta,\gamma) = R_z(\alpha)R_y(\beta)R_z(\gamma)\,, \end{align} \]

and \(I_0\) is the averaged decay rate

\[ \begin{align} I_0(m_{K\pi},m_{pK}) = \sum_{\lambda,\nu}\left|A_{\nu,\lambda}\right|^2\,. \end{align} \]

def to_index(helicity):
    """Symbolic conversion of half-value helicities to Pauli matrix indices."""
    # https://github.com/ComPWA/compwa.github.io/pull/129#issuecomment-1096599896
    return sp.Piecewise(
        (1, sp.LessThan(helicity, 0)),
        (0, True),
    )


ν_prime = sp.Symbol(R"\nu^{\prime}")
polarimetry_exprs = tuple(
    PoolSum(
        formulate_aligned_amplitude(ν, λ).conjugate()
        * msigma(i)[to_index(ν), to_index(ν_prime)]
        * formulate_aligned_amplitude(ν_prime, λ),
        (λ, [-half, +half]),
        (ν, [-half, +half]),
        (ν_prime, [-half, +half]),
    )
    / intensity_expr
    for i in (1, 2, 3)
)

Properties of the vector \(\vec\alpha\)#

The vector \(\vec \alpha\) introduced in Eq. (3) obeys the following properties:

It is a three-dimensional vector defined in the rest frame of the decaying particle. Particularly, it is transformed as a regular vector in case initial (alignment) configuration change.
The length of the vector is limited by 1: \(|\vec{\alpha}| < 1\)
\(\alpha_y=0\) for the decays of a fermion to a fermions and (pseudo)scalar

Here is the prove of the second statement:

\[\begin{split} \begin{align} I_{\nu',\nu} = \sum_{\lambda} A_{\nu',\lambda}^* A_{\nu,\lambda} = \begin{pmatrix} a & c^*\\ c & b \end{pmatrix} = \frac{a+b}{2}\left( \mathbb{I} + (\vec{\sigma} \cdot \vec{\alpha}) \right)\,, \end{align} \end{split}\]

where

\[\begin{split} \begin{align} a &= \left|A_{+,+}\right|^2+\left|A_{+,-}\right|^2\,,\\ \nonumber b &= \left|A_{-,+}\right|^2+\left|A_{-,-}\right|^2\,,\\ \nonumber c &= A_{+,+}^*A_{-,+} + A_{+,-}^*A_{-,-}\,, \end{align} \end{split}\]

and

\[ \begin{align} \alpha_x &= \frac{\text{Re}\,c}{a+b}\,,& \alpha_x &= \frac{\text{Im}\,c}{a+b}\,,& \alpha_z &= \frac{a-b}{a+b}\,, \end{align} \]

To constraint the length of the \(\vec\alpha\), one notices \(ab - c \ge 0\). Therefore,

\[ \begin{align} |\vec\alpha|^2 &= \frac{(a-b)^2+c^2}{(a+b)^2} = \frac{(a+b)^2-4ab+c^2}{(a+b)^2} \leq \frac{(a+b)^2-3ab}{(a+b)^2} \leq 1\, \end{align} \]

since \(a,b \geq 0\).

Computations with TensorWaves#

Conversion to computational backend#

The full expression tree can be converted to a computational, parametrized function as follows. Note that identify all coupling symbols are interpreted as parameters. The remaining symbols (the angles) become arguments to the function.

free_parameters = {
    symbol: value
    for symbol, value in parameter_defaults.items()
    if (symbol.name.startswith("m_") and symbol not in masses)
    or symbol.name.startswith(R"\Gamma_")
    or symbol in prod_couplings
}
fixed_parameters = {
    symbol: value
    for symbol, value in parameter_defaults.items()
    if symbol not in free_parameters
}

intensity_func = create_parametrized_function(
    unfolded_intensity_expr.xreplace(fixed_parameters),
    parameters=free_parameters,
    backend="jax",
)

polarimetry_funcs = tuple(
    create_parametrized_function(
        expr.xreplace(fixed_parameters),
        parameters=free_parameters,
        backend="jax",
    )
    for expr in unfolded_polarimetry_exprs
)

Phase space#

The \(\Lambda_c^+ \to p K \pi\) kinematics is fully described by two dynamic variables, \(m_{K\pi}\) and \(m_{pK}\) (see Phase space for a three-body decay). The third Mandelstam variable can be computed from the other two and the masses of the initial and final state:

\[\begin{split}\displaystyle \begin{array}{rcl} \sigma_{3} & = & m_{K}^{2} + m_{p}^{2} + m_{\pi}^{2} + m_{\Lambda_c}^{2} - \sigma_{1} - \sigma_{2} \\ \end{array}\end{split}\]

Values for the angles will be computed form the Mandelstam values with a data transformer for the symbolic angle definitions:

kinematic_variables = {
    symbol: expression.doit().xreplace(masses).xreplace(fixed_parameters)
    for symbol, expression in angles.items()
}
kinematic_variables.update({s: s for s in [σ1, σ2, σ3]})  # include identity
transformer = SympyDataTransformer.from_sympy(kinematic_variables, backend="jax")

We now define phase space over a grid that contains the space in the Dalitz plane that is kinematically ‘available’ to the decay:

m0_val, m1_val, m2_val, m3_val = masses.values()
σ1_min = (m2_val + m3_val) ** 2
σ1_max = (m0_val - m1_val) ** 2
σ2_min = (m1_val + m3_val) ** 2
σ2_max = (m0_val - m2_val) ** 2


def generate_phsp_grid(resolution: int):
    x = np.linspace(σ1_min, σ1_max, num=resolution)
    y = np.linspace(σ2_min, σ2_max, num=resolution)
    X, Y = np.meshgrid(x, y)
    Z = compute_third_mandelstam.function(X, Y)
    phsp = {"sigma1": X, "sigma2": Y, "sigma3": Z}
    return X, Y, transformer(phsp)


X, Y, phsp = generate_phsp_grid(resolution=500)

Intensity distribution#

Finally, all intensities can be computed as follows:

Fit fractions#

The total decay rate for \(\Lambda_c^+ \to pK\pi\) can be broken into fractions that correspond to the different decay chains and interference terms. The total rate is computed as an integral of the intensity over decay kinematics:

\[ \begin{align} I_\text{tot}(\{\mathcal{H}\}) = \int d m_{pK}^2 d m_{K\pi}^2\, I_0(m_{pK}, m_{K\pi} | \{\mathcal{H}\}) \approx \frac{\Phi_0}{N_\text{MC}} \sum_{e=1}^{N_\text{MC}}\,\,I_0(m_{pK,e}, m_{K\pi,e} | \{\mathcal{H}\})\,, \end{align} \]

where \(\Phi_0\) is an (irrelevant) constant equal to the flat phase-space integral, \((m_{pK,e}, m_{K\pi,e})\) is a vector of the kinematic variables for the \(e\)-th point in the MC sample.

The conditional argument \(\{\mathcal{H}\}\) indicates dependence of the rate on the value of the couplings. The individual fractions are found by computing the total rate for a subset of couplings set to zero,

\[\begin{split} \begin{align} I_\text{tot}^{K} &= I_\text{tot}\left(\{\mathcal{H}^{\Lambda_c^+\to\Delta^{**} K}, \mathcal{H}^{\Lambda_c^+\to\Lambda^{**} \pi} = 0\}\right)\,,\\ I_\text{tot}^{\Delta} &= I_\text{tot}\left(\{\mathcal{H}^{\Lambda_c^+\to K^{**} p}, \mathcal{H}^{\Lambda_c^+\to\Lambda^{**} \pi} = 0\}\right)\,,\\ I_\text{tot}^{\Lambda} &= I_\text{tot}\left(\{\mathcal{H}^{\Lambda_c^+\to\Delta^{**} K}, \mathcal{H}^{\Lambda_c^+\to K^{**} p} = 0\}\right)\,,\\ I_\text{tot}^{K/\Lambda} &= I_\text{tot}\left(\{\mathcal{H}^{\Lambda_c^+\to\Delta^{**} K} = 0\}\right) - I_\text{tot}^{K} - I_\text{tot}^{\Lambda}\,,\\ & \dots\,, \end{align} \end{split}\]

where the terms with a single chain index are the rate of the decay chain. The sum of all fractions should give the total rate:

\[ \begin{align} I_\text{tot}\left(\{\mathcal{H}\}\right) = \sum_{R} I_\text{tot}^{R} + \sum_{R < R'} I_\text{tot}^{R/R'} \end{align} \]

\[\begin{split}\displaystyle \begin{array}{crr} & I_\mathrm{sub}\,/\,I \\ \hline K^{**} & 0.371 \\ \Lambda^{**} & 0.051 \\ \color{gray}{\Lambda(1520)} & \color{gray}{0.031} \\ \color{gray}{\Lambda(1670)} & \color{gray}{0.020} \\ \Delta^{**} & 0.582 \\ \Delta/\Lambda & 0.013 \\ K/\Delta & -0.018 \\ K/\Lambda & 0.001 \\ \hline \mathrm{total} & 1.000 \\ \end{array}\end{split}\]

Polarimetry distributions#

\[\begin{split}\displaystyle \begin{array}{cccc} & \bar{|\alpha|} & \bar\alpha_x & \bar\alpha_y & \bar\alpha_z \\ K^{**} & 0.873 \pm 0.040 & 0.000 \pm 0.554 & 0.000 \pm 0.396 & +0.100 \pm 0.538 \\ \Lambda^{**} & 0.906 \pm 0.101 & -0.529 \pm 0.214 & 0.000 \pm 0.190 & -0.338 \pm 0.596 \\ \Delta^{**} & 0.540 \pm 0.000 & +0.320 \pm 0.153 & -0.000 \pm 0.000 & -0.324 \pm 0.246 \\ \end{array}\end{split}\]

Define sliders for the widget Hide code cell content

# Slider construction
sliders = {}
for symbol, _value in free_parameters.items():
    if symbol.name.startswith(R"\mathcal{H}"):
        real_slider = create_slider(symbol)
        imag_slider = create_slider(symbol)
        sliders[f"{symbol.name}_real"] = real_slider
        sliders[f"{symbol.name}_imag"] = imag_slider
        real_slider.description = R"\(\mathrm{Re}\)"
        imag_slider.description = R"\(\mathrm{Im}\)"
    else:
        slider = create_slider(symbol)
        sliders[symbol.name] = slider

# Slider ranges
σ3_max = (m0_val - m3_val) ** 2
σ3_min = (m1_val + m2_val) ** 2

for name, slider in sliders.items():
    slider.continuous_update = True
    slider.step = 0.01
    if name.startswith("m_"):
        if "K" in name:
            slider.min = np.sqrt(σ1_min)
            slider.max = np.sqrt(σ1_max)
        elif R"\Lambda" in name:
            slider.min = np.sqrt(σ2_min)
            slider.max = np.sqrt(σ2_max)
        elif R"\Delta" in name:
            slider.min = np.sqrt(σ3_min)
            slider.max = np.sqrt(σ3_max)
    elif name.startswith(R"\Gamma_"):
        slider.min = 0
        slider.max = max(0.5, 2 * slider.value)
    elif name.startswith(R"\mathcal{H}"):
        slider.min = -15
        slider.max = +15


# Slider values
def reset_sliders(click_event):
    for symbol, value in free_parameters.items():
        if symbol.name.startswith(R"\mathcal{H}"):
            set_slider(sliders[symbol.name + "_real"], value)
            set_slider(sliders[symbol.name + "_imag"], value)
        else:
            set_slider(sliders[symbol.name], value)


def set_coupling_to_zero(filter_pattern):
    if isinstance(filter_pattern, Combobox):
        filter_pattern = filter_pattern.value
    for name, _slider in sliders.items():
        if not name.startswith(R"\mathcal{H}"):
            continue
        if filter_pattern not in name:
            continue
        set_slider(_slider, 0)


def set_slider(slider, value):
    if slider.description == R"\(\mathrm{Im}\)":
        value = complex(value).imag
    else:
        value = complex(value).real
    n_decimals = -round(np.log10(slider.step))
    if slider.value != round(value, n_decimals):  # widget performance
        slider.value = value


reset_sliders(click_event=None)
reset_button = Button(description="Reset slider values")
reset_button.on_click(reset_sliders)

all_resonances = [r.latex for r_list in resonance_particles.values() for r in r_list]
filter_button = Combobox(
    placeholder="Enter coupling filter pattern",
    options=all_resonances,
    description=R"$\mathcal{H}=0$",
)
filter_button.on_submit(set_coupling_to_zero)

# UI design
latex = {symbol.name: sp.latex(symbol) for symbol in free_parameters}
mass_sliders = [sliders[n] for n in sliders if n.startswith("m_")]
width_sliders = [sliders[n] for n in sliders if n.startswith(R"\Gamma_")]
coupling_sliders = {}
for res_list in resonance_particles.values():
    for res in res_list:
        coupling_sliders[res.name] = (
            [s for n, s in sliders.items() if n.endswith("_real") and res.latex in n],
            [s for n, s in sliders.items() if n.endswith("_imag") and res.latex in n],
            [
                HTMLMath(f"${latex[n[:-5]]}$")
                for n in sliders
                if n.endswith("_real") and res.latex in n
            ],
        )
slider_tabs = Tab(
    children=[
        Tab(
            children=[
                VBox([HBox(s) for s in zip(*pair)])
                for pair in coupling_sliders.values()
            ],
            titles=tuple(coupling_sliders),
        ),
        VBox([HBox([r, i]) for r, i in zip(mass_sliders, width_sliders)]),
    ],
    titles=("Couplings", "Masses and widths"),
)
ui = VBox([slider_tabs, HBox([reset_button, filter_button])])

Create interactive plot Hide code cell source

%config InlineBackend.figure_formats = ['png']

fig, axes = plt.subplots(
    figsize=(12, 6.2),
    ncols=2,
    sharey=True,
)
ax1, ax2 = axes
ax1.set_title("Intensity distribution")
ax2.set_title("Polarimeter vector field")
ax1.set_xlabel(Rf"${s1_label[1:-1]}, \alpha_x$")
ax2.set_xlabel(Rf"${s1_label[1:-1]}, \alpha_x$")
ax1.set_ylabel(Rf"${s2_label[1:-1]}, \alpha_z$")
for ax in axes:
    ax.set_box_aspect(1)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False

mesh = None
quiver = None
Xα, Yα, phsp_α = generate_phsp_grid(resolution=35)
XI, YI, phsp_I = generate_phsp_grid(resolution=200)


def plot3(**kwargs):
    global quiver, mesh
    kwargs = to_complex_kwargs(**kwargs)
    for func in [*list(polarimetry_funcs), intensity_func]:
        func.update_parameters(kwargs)
    intensity = intensity_func(phsp_I)
    αx, αy, αz = tuple(func(phsp_α).real for func in polarimetry_funcs)
    abs_α = jnp.sqrt(αx**2 + αy**2 + αz**2)
    if mesh is None:
        mesh = ax1.pcolormesh(XI, YI, intensity, cmap=plt.cm.Reds)
        c_bar = fig.colorbar(mesh, ax=ax1, pad=0.01, fraction=0.0473)
        c_bar.ax.set_yticks([])
    else:
        mesh.set_array(intensity)
    if quiver is None:
        quiver = ax2.quiver(Xα, Yα, αx, αz, abs_α, cmap=plt.cm.viridis_r, clim=(0, 1))
        c_bar = fig.colorbar(quiver, ax=ax2, pad=0.01, fraction=0.0473)
        c_bar.ax.set_ylabel(R"$\left|\vec\alpha\right|$")
    else:
        quiver.set_UVC(αx, αz, abs_α)
    fig.canvas.draw_idle()


def to_complex_kwargs(**kwargs):
    complex_valued_kwargs = {}
    for key, value in dict(kwargs).items():
        if key.endswith("real"):
            symbol_name = key[:-5]
            imag = kwargs[f"{symbol_name}_imag"]
            complex_valued_kwargs[symbol_name] = complex(value, imag)
        elif key.endswith("imag"):
            continue
        else:
            complex_valued_kwargs[key] = value
    return complex_valued_kwargs


output = interactive_output(plot3, controls=sliders)
fig.tight_layout()
display(ui, output)