Amplitude Analysis 101 (PWA 101)

Amplitude Analysis 101 (PWA 101)#

Abstract

This document introduces Amplitude Analysis / Partial Wave Analysis (PWA) by demonstrating its application to a specific reaction channel and amplitude model. It aims to equip readers with a basic understanding of the full workflow and methodologies of PWA in hadron physics through a practical, hands-on example. Only basic Python programming and libraries (e.g. numpy, scipy, etc.) are used to illustrate the more fundamental steps in a PWA. Calculations with 4-vectors in this report are performed with the vector package.

Amplitude model#

Theory#

This amplitude model is adapted from the Lecture 11 in STRONG2020 HaSP School by Vincent Mathieu.

The reaction \( \gamma p \to \eta \pi^0 p\) is one of the reaction channels that are studied in photo-production experiments like GlueX. For simplicity, we assume that the decay proceeds through three possible resonances—\(a_2\), \(\Delta^+\), and \(N^*\)—with each in a different subsystem.

https://github.com/ComPWA/compwa-org/assets/17490173/ec6bf191-bd5f-43b0-a6cb-da470b071630

In the following, we denote \(1 \equiv \eta, 2 \equiv \pi^0, 3 \equiv p\). Given these three subsystems in this particle transition, we can construct three amplitudes,

(1)#\[\begin{split} \begin{eqnarray} A^{12} &=& \frac{\sum a_m Y_2^m (\Omega_1)}{s_{12}-m^2_{a_2}+im_{a_2} \Gamma_{a_2}} \times s_{12}^{0.5+0.9u_3} \,, \nonumber \\ A^{23} &=& \frac{\sum b_m Y_1^m (\Omega_2)}{s_{23}-m^2_{\Delta}+im_{\Delta} \Gamma_{\Delta}} \times s_{23}^{0.5+0.9t_1} \,, \nonumber \\ A^{31} &=& \frac{c_0}{s_{31}-m^2_{N^*}+im_{N^*} \Gamma_{N^*}} \times s_{31}^{1.08+0.2t_2} \,, \end{eqnarray} \end{split}\]

where \(s, t, u\) are the Mandelstam variables (\(s_{ij}=(p_i+p_j)^2\), \(t_i=(p_a-p_i)^2\), and \(u_i=(p_b-p_i)^2\)), \(m\) is the resonance mass, \(\Gamma\) is the resonance width, \(Y^m_l\) are spherical harmonics functions, \(\Omega_i\) are pairs of Euler angles (polar angle \(\theta\) and azimuthal angle \(\phi\)) that represent the helicity decay angles, and \(a_i\), \(b_i\), and \(c_i\) are complex-valued coefficients. Note that the Mandelstam variables and angles come from measured events, while the other symbols are parameters that need to be modified in order to have the amplitude model match the data.

Note

All that an amplitude model requires as data input, are four-momenta. For this amplitude model, there are just three required four-momenta per event, one for each final state in \(p\gamma \to \eta\pi^0 p\).

The original full amplitude model from the Lecture 11 in STRONG2020 HaSP School is shown in Equation (1). In this report, only the Breit–Wigner and Spherical harmonics terms are kept, while the exponential factor is abandoned, i.e.

(2)#\[\begin{split} \begin{eqnarray} 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{eqnarray} \end{split}\]

The intensity \(I\) that describes our measured distributions is then expressed as a coherent sum of the amplitudes \(A^{ij}\),

(3)#\[ \begin{eqnarray} I &=& |A^{12} + A^{23} + A^{31}|^2 \,. \end{eqnarray} \]

Ultimately, the amplitude model in Equations (2) and (3) in this tutorial consists of three resonances, and each of them are formed by two components: a Breit-Wigner times some spherical harmonics (\(l = 2, 1, 0\)).

Assumption: spinless final state!

While the spin of the \(\eta\) meson and the \(\pi^0\) meson are both \(0\), the spin of the proton is spin-\(\frac{1}{2}\). However, in order to simplify the amplitude model, we treat the proton as a spin-0 particle. Overall, the \(\eta\), \(\pi^0\) and \(p\) are therefore all treated as spin-0 particles.

The primary motivation for assuming proton to be spin-\(0\) is to avoid the necessity of aligning the amplitudes (see e.g. ComPWA/ampform#6). This assumption enables the intensity \(I\) to be written as a coherent sum of the amplitudes of the subsystems without the need for additional Wigner rotations. In addition, the amplitude for each decay chain contains only one spherical harmonic, namely that of the resonance.

The spherical harmonics in Equation (2) are therefore relevant only to the resonances. Here, \(l\) represents the spin of the resonances and \(m\) represents its spin projection in the decay chain. The total angular momentum and coupled spin (for the two-body state of the two decay products of the resonance) are not considered. According to [Richman, 1984] and other classical references on helicity, this is known as the helicity basis. In contrast, the canonical basis does not sum over \(L\) and \(S\), resulting in a more complex coherent sum of amplitudes. The transformation between these bases is also discussed here on the AmpForm website.

In our case:

\(A^{12}\) represents a d-wave interaction, because we assume there is a \(a_2\) resonance (spin 2) in this subsystem. The possible \(m\) values are \(−2,−1,0,1,2\). Each of these values corresponds to different orientations of the d-wave.
\(A^{23}\) represents a p-wave interaction, because we assume this subsystem has a (spin-1) \(\Delta^+\) resonance. The possible \(m\) values are \(−1,0,1\).
\(A^{31}\) represents an s-wave interaction, because we assume there is one spin-0 \(N^*\) resonance in this subsystem. The only possible \(m\) value is 0 and, since \(Y_0^0=0\), the amplitude only consists of a Breit–Wigner.

Implementation#

In the following, we define a few functions that implement parts of the amplitude model of Equation (2). Later on in Visualization, we can use these functions to visualize each components of the model individually.

Breit-Wigner (only) Model#

The following functions define the Breit–Wigner function, as well as the following intensity function that only contains the Breit–Wigner (dynamics) component of each amplitude.

(4)#\[ \begin{array}{rcl} I &= & \left\lvert \frac{1}{s-m^2_{a_2}+im_{a_2} \Gamma_{a_2}} + \frac{1}{s-m^2_{\Delta}+im_{\Delta} \Gamma_{\Delta}} + \frac{1}{s-m^2_{N^*}+im_{N^*} \Gamma_{N^*}} \right\rvert ^2 \end{array} \]

def BW_model(s12, s23, s31, *, M12, Gamma12, M23, Gamma23, M31, Gamma31, **kwargs):
    A12 = BW(s12, M12, Gamma12)
    A23 = BW(s23, M23, Gamma23)
    A31 = BW(s31, M31, Gamma31)
    return np.abs(A12 + A23 + A31) ** 2


def BW(s, m, Gamma):
    return 1 / (s - m**2 + m * Gamma * 1j)

Spherical Harmonics (only) Model#

The calculation of \(Y_l^m(\phi, \theta)\) is done via scipy.special.sph_harm(). However, we use a different definition of \(\phi\) and \(\theta\), following a notation that is more common in hadron physics.

\(\phi\) is the azimuthal angle and ranges from -\(\pi\) to \(\pi\). SciPy represents this as \(\theta\), ranging from \(0\) to \(2\pi\).
\(\theta\) is the polar angle and ranges from \(0\) to \(\pi\). SciPy represents this as \(\phi\) with the same range.

This leads to

(5)#\[ Y_\ell^m(\phi, \theta) = \sqrt{\frac{2n+1}{4\pi}\frac{(n-m)!}{(n+m)!}}e^{im\phi}P_\ell^m\left(\cos\theta\right) \,. \]

Tip

Alternatively, we can formulate the spherical harmonics in terms of a Wigner-\(d\) or Wigner-\(D\) function, as

(6)#\[ \begin{eqnarray} Y_\ell^m(\theta,\phi) = \sqrt{\frac{2\ell+1}{4\pi}}D^\ell_{m0}(-\phi, \theta, 0) = \sqrt{\frac{2\ell+1}{4\pi}}e^{im\phi} d^\ell_{m0}(\theta) \,. \end{eqnarray} \]

In the following, we define functions to compute the spherical harmonics of the three subsystem amplitudes in Equation (2). Note how the function signature consists of two input data columns, theta and phi, and how the rest of the arguments are parameters. The final kwargs (key-word arguments) is there so that we can compose larger functions from these function definitions.

def Ylm12(
    theta: np.ndarray,
    phi: np.ndarray,
    *,
    a_minus2,
    a_minus1,
    a_0,
    a_plus1,
    a_plus2,
    **kwargs,
) -> np.ndarray:
    return (
        a_plus2 * sp.special.sph_harm(2, 2, phi, theta)
        + a_plus1 * sp.special.sph_harm(1, 2, phi, theta)
        + a_0 * sp.special.sph_harm(0, 2, phi, theta)
        + a_minus1 * sp.special.sph_harm(-1, 2, phi, theta)
        + a_minus2 * sp.special.sph_harm(-2, 2, phi, theta)
    )

def Ylm23(
    theta: np.ndarray, phi: np.ndarray, *, b_minus1, b_0, b_plus1, **kwargs
) -> np.ndarray:
    return (
        b_plus1 * sp.special.sph_harm(1, 1, phi, theta)
        + b_0 * sp.special.sph_harm(0, 1, phi, theta)
        + b_minus1 * sp.special.sph_harm(-1, 1, phi, theta)
    )

def Ylm31(theta: np.ndarray, phi: np.ndarray, c_0: complex, **kwargs) -> np.ndarray:
    return c_0 * sp.special.sph_harm(0, 0, phi, theta)

We now have the ingredients to define an intensity function that only contains spherical harmonics, that is, the angular part of the amplitude model:

(7)#\[ \begin{array}{rcl} I &=& \left\lvert \sum a_m Y_2^m (\Omega_1) + \sum b_m Y_1^m (\Omega_2) + c_0 \right\rvert ^2 \,. \end{array} \]

def SH_model(phi1, theta1, phi2, theta2, *, c_0, **pars):
    return np.abs(Ylm12(phi1, theta1, **pars) + Ylm23(phi2, theta2, **pars) + c_0) ** 2

Breit-Wigner \(\times\) Spherical Harmonics Model#

Finally, we combine the Breit-Wigner (only) Model and Spherical Harmonics (only) Model to get an implementation for Equation (2),

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

def BW_SH_model(
    s12,
    s23,
    s31,
    phi1,
    theta1,
    phi2,
    theta2,
    *,
    M12,
    Gamma12,
    M23,
    Gamma23,
    M31,
    Gamma31,
    c_0,
    **parameters,
):
    A12 = BW(s12, M12, Gamma12) * Ylm12(phi1, theta1, **parameters)
    A23 = BW(s23, M23, Gamma23) * Ylm23(phi2, theta2, **parameters)
    A31 = BW(s31, M31, Gamma31) * c_0
    return np.abs(A12 + A23 + A31) ** 2

Visualization#

As discussed in Theory, the amplitude model contains variables that are provided by experiment (“events” as data input) as well as parameters that are to be optimized. The data input is usually in the form of four-momenta, while the model is formulated in the form of Mandelstam variables and helicity angles. We therefore have to compute these variables from the four-momenta.

Phase space generation#

First, however, we need to generate a phase space sample of four-momenta in order to plot the amplitude model as a distribution over each of the variables. In this section, we use the phasespace package for generating four-momenta for the reaction \(p\gamma \to p\eta\pi^0\). The phase space sample will also be used later on to normalize the model when calculating the likelihood over the data sample (see Fitting).

Lab frame and CM frame#

In Center-of-Mass (CM) frame, the 4-momentum of the total system can be acquired by 4-momentum conservation:

(8)#\[\begin{split}\begin{pmatrix} E_0 \\ 0 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} E_{\gamma} \\ 0 \\ 0 \\ p_z \end{pmatrix} + \begin{pmatrix} E_p \\ 0 \\ 0 \\ -p_z \end{pmatrix} = \begin{pmatrix} E_{\eta} \\ p_{\eta,x} \\ p_{\eta,y} \\ p_{\eta,z} \end{pmatrix} + \begin{pmatrix} E_{\pi} \\ p_{\pi,x} \\ p_{\pi,y} \\ p_{\pi,z} \end{pmatrix} + \begin{pmatrix} E_p \\ p_{p,x} \\ p_{p,y} \\ p_{p,z} \end{pmatrix} \end{split}\]

Caution

The calculation here involved using values from the CM frame. While this frame is commonly used for theoretical calculations, experimental data is often analyzed in the lab frame. However, it’s worth noting that oin some collider experiments that do not have a fixed target, the CM frame can coincide with the lab frame.

The GlueX experiment at Jefferson Lab uses a fixed proton target with a linearly polarized photon beam, and the beam energy range in the lab frame is typically from 8 to 9 GeV.

We use the \(\gamma\) beam energy in the lab frame as input, e.g. \(E_{\gamma, \text{lab}} = 8.5 \; \text{GeV}\), and want to know the collision energy in the CM frame.

We can calculate the energy of the photon in the CM frame as follows. The four-momentum of photon (beam) in the lab frame is

(9)#\[ p_{\gamma,\text{lab}} = \left(E_{\gamma, \text{lab}}, \vec{p}_{\gamma,\text{lab}}\right) \,. \]

Since the photon is massless, \(m_\gamma=0\), and

(10)#\[ m^2 = E^2 - \left|\vec{p}\right|^2 \,, \]

we get

(11)#\[ E_{\gamma} = |\vec{p_{\gamma}}| . \]

The proton (target) is at rest in the lab frame, so we have

(12)#\[ p_{p,\text{lab}} = (m_p, \vec{0})\,, \]

where \(m_p\) is the proton mass. We have a total four-momentum in the lab frame of

(13)#\[ p_{\text{tot},\text{lab}} = p_{\gamma,\text{lab}} + p_{p,\text{lab}} \]

and we know

(14)#\[ E_p = \sqrt{p_p^2+ m_p^2} \,. \]

The CM total energy \(E_0\) expressed in terms of quantities from the lab frame is thus

(15)#\[ m_0 \equiv E_0 = \sqrt{2 E_{\gamma,\text{lab}} m_p + m_p^2} \,. \]

Equivalently, from the CM frame perspective, since \(\vec{p}_{\gamma} = -\vec{p}_{p}\) and \(|\vec{p}_{\gamma}| = |\vec{p}_{p}|= p_{z}\), we find

(16)#\[ E_0 = m_0 = |p_{z}|+\sqrt{p_{z}^2 + m_p^2} \,. \]

Our implementation based on Equation (15) thus becomes:

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

Thus, we then have the mass of the system \(m_0\) (or the mass of a ‘virtual’ particle \(p\gamma\)) in CM frame of

\[ E_{0} = m_0\approx 4.102 \;\; \text{GeV} \,. \]

Final state four-momenta#

The phasespace library is a Python package designed to simulate particle decays according to the principles of relativistic kinematics and phase space distributions. We first use the phasespace to generate decay particles.

phsp_events = 500_000
weights, particles = phasespace.nbody_decay(m_0, [m_eta, m_pi, m_proton]).generate(
    n_events=phsp_events
)

Note that each event comes with a statistical weights for each generated event. These weights represent how likely each particular event configuration (set of momenta and energies) is, based on phase space considerations. In order to generate a flat distribution, we will have to use a hit-and-miss method over these weights.

def generate_phsp_decay(
    size: int, seed: int | None = None, bunch_size: int = 500_000
) -> tuple[MomentumNumpy4D, MomentumNumpy4D, MomentumNumpy4D]:
    rng = np.random.default_rng(seed=seed)
    generated_seed = rng.integers(1_000_000)
    phsp_sample = generate_phsp_bunch(bunch_size, generated_seed)
    while get_size(phsp_sample) < size:
        bunch = generate_phsp_bunch(bunch_size, generated_seed)
        phsp_sample = concatenate(phsp_sample, bunch)
    phsp_sample = remove_overflow(phsp_sample, size)
    return tuple(to_vector(tensor) for tensor in phsp_sample)


def generate_phsp_bunch(
    size: int, seed: int | None = None
) -> tuple[tf.Tensor, tf.Tensor, tf.Tensor]:
    rng = np.random.default_rng(seed=seed)
    final_state = [m_eta, m_pi, m_proton]
    weights, particles = phasespace.nbody_decay(m_0, final_state).generate(
        n_events=size, seed=seed
    )
    random_weights = rng.uniform(0, weights.numpy().max(), size=weights.shape)
    selector = weights > random_weights
    return tuple(particles[f"p_{i}"][selector] for i in range(len(final_state)))


def get_size(phsp: tuple[tf.Tensor, tf.Tensor, tf.Tensor]):
    return len(phsp[0])


def concatenate(
    phsp1: tuple[tf.Tensor, tf.Tensor, tf.Tensor],
    phsp2: tuple[tf.Tensor, tf.Tensor, tf.Tensor],
) -> tuple[tf.Tensor, tf.Tensor, tf.Tensor]:
    return tuple(tf.concat([phsp1[i], phsp2[i]], axis=0) for i in range(3))


def remove_overflow(phsp: tuple[tf.Tensor, tf.Tensor, tf.Tensor], size: int):
    return tuple(tensor[:size] for tensor in phsp)


def to_vector(tensor: tf.Tensor) -> MomentumNumpy4D:
    return vector.array({
        key: tensor.numpy().T[j] for j, key in enumerate(["px", "py", "pz", "E"])
    })

Hit-and-miss sampling

Step 1: Generate a small phase space sample bunch with generate_phsp_bunch()
- A NumPy random number generator (rng) is initialized with a fixed seed for reproducibility. This ensures that every time the simulation is run, it produces the same random numbers.
- The phasespace package is used to simulate a decay process where a parent particle decays into three daughter particles. This is set up by specifying the masses of the parent particle (m_0) and the daughter particles (m_eta, m_pi, m_proton). The function nbody_decay followed by generate is used to create the phase space and return both the weights of each decay event and the momenta of the decay products.
- As mentioned, the phasespace package generates four-momenta that are not evenly distributed .The distributions appear even only when multiplied by the generated weights. To get a sample bunch that is uniformly distributed, rejection sampling (“hit-and-miss”) is applied over the weights. To achieve this, a random set of weights (random_weights) is generated, and only those events where the original weight exceeds the random weight are selected.
Step 2, Ensure sufficient sample size with generate_phsp_decay()
- The function starts by generating an initial phase space sample. If this initial sample does not meet the required size, more samples are generated and concatenated until the requested size is reached.
- Once the target sample size is reached or exceeded, remove_overflow() trims the sample to precisely match the requested size.
Step 3, Convert tensors to four-momentum vectors objects The phasespace package generates four-momenta in the form of \(4\times N\) TensorFlow tensors, with \(N\) the number of events. For each tensor in the phase space tuple (representing the momentum components like \(E, p_x, p_y, p_z\)), the to_vector function converts them into a structured MomentumNumpy4D array. This structured array includes keys for each momentum component and possibly the energy component and provides convenient methods for computing boosts etc.

Initial state four-momenta#

Because of the simple relation in Equation (8), the four-momenta for the initial state, \(\gamma\) (a) and \(p\) (b), do not have to be generated by a phase space generator.

To find \(p_a\) and \(p_b\) by 4-momentum conservation, we can use

\[\begin{split} \begin{array}{rcl} m_a &=& m_{\gamma} =0 \\ m_b &=& m_p \\ p_{a,x} &=& p_{a,y} = 0 = p_{b,x} = p_{b,y} \\ p_{a,z} &=& - p_{b,z} \,. \end{array} \end{split}\]

Due to energy conservation, we have

\[\begin{split} \begin{array}{rcl} E_a + E_b &=& E_1 + E_2 + E_3 = E_0 \\ E_a + E_b &=& \sqrt{m_a^2 + p_a^2} + \sqrt{m_b^2 + p_b^2} = E_0 \,. \end{array} \end{split}\]

Since \(m_a=0\) and \(p_a = -p_b\), we have

\[ p_a + \sqrt{m_b^2 + (-p_a)^2} = E_0 \,. \]

Reorganizing, we get,

\[ p_{a,z} + \sqrt{m_p^2 + p_a^2} - E_0 = 0 \]

and

\[ p_{a,z} = \frac{E_0^2 - m_p^2}{2E_0} \]

Finally, this gives us

(17)#\[\begin{split} \begin{array}{rcl} p_{b,z} &=& -p_{a,z} \\ E_a &=& p_a \\ E_{b} &=& \sqrt{m_b^2+p_b^2} \,. \end{array} \end{split}\]

With this, our numerical implementation to compute four-momentum of \(p_a\) and \(p_b\) from Equation (17) becomes:

def compute_pa_pb(
    p1: MomentumNumpy4D, p2: MomentumNumpy4D, p3: MomentumNumpy4D
) -> tuple[MomentumNumpy4D, MomentumNumpy4D]:
    shape = p1.shape
    E0 = (p1 + p2 + p3).e
    px = np.zeros(shape)
    py = np.zeros(shape)
    pz = np.ones(shape) * (E0**2 - p3.m.mean() ** 2) / (2 * E0)
    E = np.ones(shape) * np.sqrt(p3.m.mean() ** 2 + pz.mean() ** 2)
    pa = MomentumNumpy4D({"E": pz, "px": px, "py": py, "pz": pz})
    pb = MomentumNumpy4D({"E": E, "px": px, "py": py, "pz": -pz})
    return pa, pb

Four-momenta of all particles#

We now create a function generate_phsp_all() to create and compute all particles in phase space all-at-once, combining the two functions from the previous sections into one function.

def generate_phsp_all(
    size: int, seed: int | None = None, bunch_size: int = 500_000
) -> tuple[
    MomentumNumpy4D, MomentumNumpy4D, MomentumNumpy4D, MomentumNumpy4D, MomentumNumpy4D
]:
    p1, p2, p3 = generate_phsp_decay(size, seed, bunch_size)
    pa, pb = compute_pa_pb(p1, p2, p3)
    return p1, p2, p3, pa, pb

%%time
p1_phsp, p2_phsp, p3_phsp, pa_phsp, pb_phsp = generate_phsp_all(
    size=phsp_events, seed=42, bunch_size=1_000_000
)

Kinematic variable calculation#

Spherical coordinate system#

Before introducing CM and helicity angles, we first introduce polar angles and azimuthal angles in spherical coordinates.

In the spherical coordinates, the polar angle \(\theta\) and azimuthal angle \(\phi\) are defined as

(18)#\[\begin{split} \begin{array}{rcl} \theta &=& \arccos \frac{p_z}{|p|} \nonumber \\ \phi &=& \arctan2(p_y , p_x) \,. \end{array} \end{split}\]

https://github.com/ComPWA/compwa.github.io/assets/29308176/89e1ad6e-3a7b-4895-b747-8bd644652504

In Equation (18), \(p_z\) is equivalent to \(z\), and \(|p|\) is \(r\) in figure above, while \(p_y\) equivalent to \(y\), and \(p_x\) is \(x\).

The sample is plotted to check whether the distribution looks uniformly distributed in the Dalitz plane. The Mandelstam variable \(s\) of each of the three subsystems can be easily computed with from the four-momentum objects as follows. Here, we use the methods and attributes provided by the vector package.

p12_phsp = p1_phsp + p2_phsp
p23_phsp = p2_phsp + p3_phsp
p31_phsp = p3_phsp + p1_phsp
s12_phsp = p12_phsp.m2
s23_phsp = p23_phsp.m2
s31_phsp = p31_phsp.m2

Tip

There are different ways to represent the Dalitz plot, each with its advantages.

Scatter Plot: This method plots individual events as points, offering a clear view of the density and distribution of events within the phase space. It is particularly useful for visualizing smaller datasets or when high resolution is needed to identify specific features or clusters.
2D Histogram: This approach divides the phase space into bins and counts the number of events within each bin, representing the density of events using a color scale. It is effective for large datasets, providing a smooth and continuous representation of the phase space that highlights overall trends and structures.

In conclusion, to check if the phase space sample is evenly distributed, a scatter plot is typically more straightforward and visually clear.

CM Angles#

Angles in the CM frame are the polar and azimuthal angles of the spatial components of the four-momenta in the CM frame (i.e. the frame that satisfies the relations in Equation (8)). They are different than the helicity angles in each subsystem (which is after rotation and boost into the subsystem).

The values for phase space can be obtained directly in the CM frame, without boosting into a different frame after the generation. We denote these angles as

\[ \theta_1^\text{CM}, \theta_2^\text{CM}, \theta_3^\text{CM}, \phi_1^\text{CM}, \phi_2^\text{CM}, \phi_3^\text{CM} \,. \]

theta1_CM_phsp = p1_phsp.theta
theta2_CM_phsp = p2_phsp.theta
theta3_CM_phsp = p3_phsp.theta
phi1_CM_phsp = p1_phsp.phi
phi2_CM_phsp = p2_phsp.phi
phi3_CM_phsp = p3_phsp.phi

Helicity angles#

Helicity

Helicity \(\lambda\) is the projection of a particle’s spin \(S\) onto its direction of motion \(\vec{p}\). Mathematically, it’s given by

\[\lambda = \frac{S \cdot \vec{p}}{|\vec{p}|} \,.\]

For massless particles (photons in particular), helicity is a Lorentz-invariant quantity, meaning it remains unchanged under boosts along the direction of motion.

Helicity state

Helicity states are eigenstates of the helicity operator. For a particle moving in the \(z\) direction, the helicity operator is \(S \cdot \hat{p}\).

Wigner-\(D\) matrix

Wigner-\(D\) matrices are used in the helicity formalism to describe the rotation of spin states. These matrices depend on the Euler angles of the rotation and are denoted by

\[D^j_{m'm} (\alpha,\beta,\gamma) = e^{-im'\alpha} d^j_{m'm}(\beta) e^{-im\gamma} \,,\]

where \(j\) is the spin, \(m'\) and \(m\) are the magnetic quantum numbers, \(\alpha, \beta, \gamma\) are the Euler angles of the rotation, and \(d^j_{m'm}(\beta)\) is an element of the orthogonal Wigner’s (small) \(d\)-matrix.

As can be seen in the amplitude model of Equation (2), that the rotations (represented by \(Y^m_l\)) contain solid angles \(\Omega=\phi,\theta\) (see spherical coordinates). They have to be computed in the helicity frame of the resonance, which means we have to boost into each of the three subsystems. For instance, for the \(a_2\) resonance (\(A^{12}\)), this would be a boost into subsystem \(p_1+p_2\), plus a rotation such that the \(z\) axis points in the direction of \(p_1+p_2\). The helicity angles \(\phi\) and \(\theta\) can then easily be computed from the spherical coordinates of the (boosted) \(p_1'\) or \(p_2'\).

Reference frames after boost

CM frame	Helicity frame
Before boosting	After boosting into the \(p_1+p_2\) subsystem

The production plane (cyan) is defined by the momenta of the incoming particles that participate in the production of a particular state or particle. In the figure, we have \(a+b \to (R\to 1+2)+3\) (subsystem \(p_1+p_2\)), meaning that the production plane is spanned by the momenta of \(a\) (or \(b\)) and \(1\).

The decay plane (magenta) is defined by the momenta of the decay products of a particle. In the figure, the production plane is spanned by the momenta of \(1\) and \(2\).

Numerical angle computation#

To calculate the helicity angles \(\theta\) and \(\phi\), we define functions for boosting a combination of boost and rotation (around the \(y\) and \(z\) axis).

def theta_helicity(p_i: MomentumNumpy4D, p_ij: MomentumNumpy4D) -> np.ndarray:
    return _boost_to_helicity_frame(p_i, p_ij).theta


def phi_helicity(p_i: MomentumNumpy4D, p_ij: MomentumNumpy4D) -> np.ndarray:
    return _boost_to_helicity_frame(p_i, p_ij).phi


def _boost_to_helicity_frame(
    p_i: MomentumNumpy4D, p_ij: MomentumNumpy4D
) -> MomentumNumpy4D:
    return p_i.rotateZ(-p_ij.phi).rotateY(-p_ij.theta).boostZ(-p_ij.beta)

theta1_phsp = theta_helicity(p1_phsp, p12_phsp)
theta2_phsp = theta_helicity(p2_phsp, p23_phsp)
theta3_phsp = theta_helicity(p3_phsp, p31_phsp)
phi1_phsp = phi_helicity(p1_phsp, p12_phsp)
phi2_phsp = phi_helicity(p2_phsp, p23_phsp)
phi3_phsp = phi_helicity(p3_phsp, p31_phsp)

The distribution of the phase space sample over the kinematic variables is shown later in this section, together the generated data and weighted phase space (model).

Toy model parameter values#

Now that we have a numerical implementation for the amplitude model of Equation (2) and a means to compute the kinematic variables appearing in that expression, we are ready to visualize the model. First, however, we have to decide on some toy values for the parameters in the model. The toy model parameter values can be obtained from the data file from the Lecture 11 in STRONG2020 HaSP School. In this tutorial, the values are modified to make the structures in the Dalitz plot more visible.

toy_parameters = dict(
    a_minus2=0,
    a_minus1=0.5,
    a_0=3.5,
    a_plus1=4,
    a_plus2=2.5,
    b_minus1=-1.5,
    b_0=4,
    b_plus1=0.5,
    c_0=2.5,
    M12=1.32,
    Gamma12=0.1,
    M23=1.24 + 0.3,
    Gamma23=0.1,
    M31=1.57 + 0.3,
    Gamma31=0.1,
)

Note that the masses \((M)\) and widths \((\Gamma)\) are properties of the three resonances.

(19)#\[\begin{split} \begin{array}{rcl} M_{12} &=& M_{\eta\pi^0} = M_{a_2} \\ M_{23} &=& M_{\pi^0p} = M_{\Delta^+} \\ M_{31} &=& M_{\eta p} = M_{N^*} \\ \Gamma_{12} &=& \Gamma_{\eta\pi^0} = \Gamma_{a_2} \\ \Gamma_{23} &=& \Gamma_{\pi^0p} = \Gamma_{\Delta^+} \\ \Gamma_{31} &=& \Gamma_{\eta p} = \Gamma_{N^*} \\ \end{array} \end{split}\]

Spherical harmonics#

We now have a look at the real part and imaginary part of \(\sum a_m Y_2^m (\Omega_1)\) as well as \(\sum b_m Y_1^m (\Omega_2)\) in Equation (2). For this, we define a grid of values for \(\phi\) and \(\theta\) over which to visualize the amplitudes.

PHI, THETA = np.meshgrid(
    np.linspace(-np.pi, +np.pi, num=1_000),
    np.linspace(0, np.pi, num=1_000),
)
Z12 = Ylm12(PHI, THETA, **toy_parameters)
Z23 = Ylm23(PHI, THETA, **toy_parameters)

Dalitz Plots of (sub)models#

Dalitz plots are ideal for visualizing the model over the Mandelstam variables \(s_{12}\), \(s_{23}\), and \(s_{31}\). In the following, we plot each of the ‘sub-models’ separately: Breit-Wigner (only), Spherical Harmonics (only), and Breit-Wigner \(\times\) Spherical Harmonics.

To integrate out the model dependence on the helicity angles, we plot the Dalitz plots as a histogram over the phase space sample, taking the computed intensities for each event in the sample as weight.

fig

Data Generation#

In an actual amplitude analysis, we are now ready to ‘fit’ the model we formulated to a measured data distribution. In this tutorial, however, we generate the data ourselves and then perform a test fit using a starting model with some ‘guessed’ parameter values.

Hit and miss intensity sample#

The data can be generated with a similar hit-and-miss strategy as we saw in Phase space generation. In this case, the hit-and-miss is performed over the intensities computed from the model.

The following function generates a data sample based on a model by using hit-and-miss filter applied to the phase space sample. The output is a tuple of five four-momenta \((p_1, p_2, p_3, p_a, p_b)\).

def generate_data(
    model: callable, size: int, bunch_size: int = 500_000, seed: int | None = None
) -> tuple[MomentumNumpy4D, ...]:
    rng = np.random.default_rng(seed=seed)
    generated_seed = rng.integers(low=0, high=1_000_000)
    progress_bar = tqdm(total=size)
    phase_space = generate_phsp_all(bunch_size, generated_seed)
    data_sample = _hit_and_miss_filter(model, phase_space, generated_seed)
    progress_bar.n = min(get_size(data_sample), size)
    progress_bar.update(n=0)
    while get_size_data(data_sample) < size:
        phase_space = generate_phsp_all(bunch_size, generated_seed)
        bunch = _hit_and_miss_filter(model, phase_space)
        data_sample = concatenate_data(data_sample, bunch)
        progress_bar.n = min(get_size(data_sample), size)
        progress_bar.update(n=0)
    progress_bar.close()
    return remove_overflow_data(data_sample, size)


def _hit_and_miss_filter(
    model: callable, phase_space: tuple[MomentumNumpy4D, ...], seed: int | None = None
) -> tuple[MomentumNumpy4D, ...]:
    p1, p2, p3, pa, pb = phase_space
    p12 = p1 + p2
    p23 = p2 + p3
    p31 = p3 + p1

    intensities: np.ndarray = model(
        s12=p12.m2,
        s23=p23.m2,
        s31=p31.m2,
        phi1=phi_helicity(p1, p12),
        theta1=theta_helicity(p1, p12),
        phi2=phi_helicity(p2, p23),
        theta2=theta_helicity(p2, p23),
        **toy_parameters,
    )
    rng = np.random.default_rng(seed=seed)  # FIX seed is used here for reproducibility
    random_intensities = rng.uniform(0, intensities.max(), size=intensities.shape)
    selector = intensities > random_intensities
    return (
        p1[selector],
        p2[selector],
        p3[selector],
        pa[selector],
        pb[selector],
    )


def get_size_data(
    data: tuple[
        MomentumNumpy4D,
        MomentumNumpy4D,
        MomentumNumpy4D,
        MomentumNumpy4D,
        MomentumNumpy4D,
    ],
):
    return len(data[0])


def concatenate_data(
    data1: tuple[MomentumNumpy4D, ...],
    data2: tuple[MomentumNumpy4D, ...],
) -> tuple[MomentumNumpy4D, ...]:
    return tuple(
        concatenate_vectors((pi1, pj2)) for pi1, pj2 in zip(data1, data2, strict=True)
    )


def concatenate_vectors(vectors: tuple[MomentumNumpy4D]) -> MomentumNumpy4D:
    return vector.array({
        "px": np.concatenate([p.px for p in vectors]),
        "py": np.concatenate([p.py for p in vectors]),
        "pz": np.concatenate([p.pz for p in vectors]),
        "E": np.concatenate([p.e for p in vectors]),
    })


def remove_overflow_data(
    data: tuple[MomentumNumpy4D, ...], size: int
) -> tuple[MomentumNumpy4D, ...]:
    return tuple(momentum[:size] for momentum in data)

Hit-and-miss data generation

A few functions are defined above for generating data sample from filtering phase space samples of particle decays and converting them into four-momentum vectors. The main ideas of them are:

generate_data(): The main function used to generate the data
- It simulates events until it gathers a specified number of valid data points that meet the model’s conditions.
- Parameters:
  - model: A callable that computes the intensity of events based on given kinematic variables.
  - size: The desired number of valid events to generate.
  - bunch_size: The number of events to generate per iteration if the initial batch doesn’t meet the required size. Defaults to size if not provided.
  - seed: Seed for random number generation to ensure reproducibility.
- Process:
  - Initialize a random number generator.
  - Generate an initial set of phase space data (phase_space) and filter it using the model to obtain valid events (data_sample).
  - If the filtered data (data_sample) is less than the desired size, additional data is generated in batches (bunch_size) and filtered until the size requirement is met.
  - The function ensures the returned dataset size matches exactly the requested size.
_hit_and_miss_filter: The other important function that filters the generated phase space data
- using the hit-and-miss Monte Carlo method based on the model intensities.
- Parameters:
  - model: The model function to calculate intensities.
  - phase_space: Tuple containing momentum vectors of particles involved in the event.
  - seed: Seed for random number generation.
- Process:
  - Calculate invariant masses and other kinematic variables required by the model.
  - Call the model with these variables to compute intensities for each event.
  - Generate random numbers and select events where the model intensity is greater than a random threshold, simulating a probability proportional to the intensity.

There are other assisted functions that are used in the functions above:

get_size_data:
- Simply returns the number of events in the data, used to check if more data needs to be generated.
concatenate_data: - Merges two datasets into one, maintaining the structure of tuples of momentum vectors.
concatenate_vectors: - Helper function to concatenate individual momentum vectors within the data tuples.
remove_overflow_data: - Trims the dataset to the desired size if it exceeds the specified number due to batch generation.

Now we can use the function generate_data() to generate our data sample:

data_events = 100_000
data = generate_data(BW_SH_model, data_events, seed=0, bunch_size=1_000_000)
p1_data, p2_data, p3_data, pa_data, pb_data = data

Dalitz plot of data distribution#

1D projections#

The 1D projection distribution of CM angles, the helicity angles and invariant mass of the model, phase space, and data are shown in this section.

Fitting#

We now pretend that we do not know the parameter values for the data distribution. In this particular situation, the parameters (to be fitted) are mass \(M\) (\(M_{12}\), \(M_{23}\), and \(M_{31}\)) and decay width \(\Gamma\) (\(\Gamma_{12}\), \(\Gamma_{23}\), and \(\Gamma_{31}\)). We also guess the coefficient values, but note that we keep one coefficient ‘fixed’. The reason is that the model is normalized through the likelihood estimator, so what matters for the fits are the ratios between the parameters.

Guessed parameter values#

guessed_parameters = dict(
    M12=1.4,
    Gamma12=0.15,
    M23=1.59,
    Gamma23=0.05,
    M31=1.82,
    Gamma31=0.2,
    # a_minus2 is set to fixed and not changed during fitting
    a_minus1=0.51,
    a_0=3.49,
    a_plus1=4.08,
    a_plus2=2.6,
    b_minus1=-1.2,
    b_0=4.04,
    b_plus1=0.41,
    c_0=2.66,
)

Function for plotting model Hide code cell source

def compare_model_to_data(parameters: dict):
    fig, axes = plt.subplots(figsize=(13, 12), ncols=3, nrows=5)
    fig.suptitle(
        R"CM angles, Helicity angles and invariant mass (Before fitting: starting"
        R" parameters for BW $\times$ SH model)",
        y=1,
    )

    thetaCM_ax, phiCM_ax, theta_ax, phi_ax, mass_ax = axes
    for i, ax in enumerate(thetaCM_ax, 1):
        ax.set_title(thetaCM_subtitles[i - 1])
        ax.set_xticks([-1, 0, 1])
    for i, ax in enumerate(phiCM_ax, 1):
        ax.set_title(phiCM_subtitles[i - 1])
        ax.set_xticks([-np.pi, 0, np.pi])
        ax.set_xticklabels([R"-$\pi$", 0, R"$\pi$"])
    for i, ax in enumerate(theta_ax, 1):
        ax.set_title(theta_subtitles[i - 1])
        ax.set_xticks([-1, 0, 1])
    for i, ax in enumerate(phi_ax, 1):
        ax.set_title(phi_subtitles[i - 1])
        ax.set_xticks([-np.pi, 0, np.pi])
        ax.set_xticklabels([R"-$\pi$", 0, R"$\pi$"])
    for i, ax in enumerate(mass_ax, 1):
        ax.set_title(mass_subtitles[i - 1])

    data_kwargs = dict(bins=100, label="data", density=True)
    data_phi_kwargs = dict(bins=50, label="data", density=True)

    thetaCM_ax[0].hist(np.cos(theta1_CM_data), **data_kwargs, alpha=0.5)
    thetaCM_ax[1].hist(np.cos(theta2_CM_data), **data_kwargs, alpha=0.5)
    thetaCM_ax[2].hist(np.cos(theta3_CM_data), **data_kwargs, alpha=0.5)
    theta_ax[0].hist(np.cos(theta1_data), **data_kwargs)
    theta_ax[1].hist(np.cos(theta2_data), **data_kwargs)
    theta_ax[2].hist(np.cos(theta3_data), **data_kwargs)
    phiCM_ax[0].hist(phi1_CM_data, **data_phi_kwargs, alpha=0.5)
    phiCM_ax[1].hist(phi2_CM_data, **data_phi_kwargs, alpha=0.5)
    phiCM_ax[2].hist(phi3_CM_data, **data_phi_kwargs, alpha=0.5)
    phi_ax[0].hist(phi1_data, **data_phi_kwargs)
    phi_ax[1].hist(phi2_data, **data_phi_kwargs)
    phi_ax[2].hist(phi3_data, **data_phi_kwargs)
    mass_ax[0].hist(p12_data.m, **data_kwargs)
    mass_ax[1].hist(p23_data.m, **data_kwargs)
    mass_ax[2].hist(p31_data.m, **data_kwargs)

    intensities = BW_SH_model(
        s12_phsp,
        s23_phsp,
        s31_phsp,
        phi1_phsp,
        theta1_phsp,
        phi2_phsp,
        theta2_phsp,
        **{**toy_parameters, **parameters},
    )
    model_kwargs = dict(
        bins=100,
        weights=intensities,
        color="red",
        histtype="step",
        label="BW&SH",
        density=True,
    )
    model_phi_kwargs = dict(
        bins=50,
        weights=intensities,
        color="red",
        histtype="step",
        label="BW&SH",
        density=True,
    )
    thetaCM_ax[0].hist(np.cos(theta1_CM_phsp), **model_kwargs)
    thetaCM_ax[1].hist(np.cos(theta2_CM_phsp), **model_kwargs)
    thetaCM_ax[2].hist(np.cos(theta3_CM_phsp), **model_kwargs)
    phiCM_ax[0].hist(phi1_CM_phsp, **model_phi_kwargs)
    phiCM_ax[1].hist(phi2_CM_phsp, **model_phi_kwargs)
    phiCM_ax[2].hist(phi3_CM_phsp, **model_phi_kwargs)
    theta_ax[0].hist(np.cos(theta1_phsp), **model_kwargs)
    theta_ax[1].hist(np.cos(theta2_phsp), **model_kwargs)
    theta_ax[2].hist(np.cos(theta3_phsp), **model_kwargs)
    phi_ax[0].hist(phi1_phsp, **model_phi_kwargs)
    phi_ax[1].hist(phi2_phsp, **model_phi_kwargs)
    phi_ax[2].hist(phi3_phsp, **model_phi_kwargs)
    mass_ax[0].hist(p12_phsp.m, **model_kwargs)
    mass_ax[1].hist(p23_phsp.m, **model_kwargs)
    mass_ax[2].hist(p31_phsp.m, **model_kwargs)

    for ax in axes.flatten():
        ax.legend()
    fig.tight_layout()
    plt.show()

Estimator: likelihood function#

To optimize the model parameters so that the model ‘fits’ the data distribution, we need to define a measure for the ‘distance’ between the model to the data. The common choice for complicated, multidimensional models like an amplitude model is the unbinned negative log likelihood function (NLL).

def unbinned_nll(**parameters: float) -> float:
    phsp = (
        s12_phsp,
        s23_phsp,
        s31_phsp,
        phi1_phsp,
        theta1_phsp,
        phi2_phsp,
        theta2_phsp,
    )
    data = (
        s12_data,
        s23_data,
        s31_data,
        phi1_data,
        theta1_data,
        phi2_data,
        theta2_data,
    )
    model_integral = BW_SH_model(*phsp, **parameters).mean()
    data_intensities = BW_SH_model(*data, **parameters)
    likelihoods = data_intensities / model_integral
    log_likelihood = np.log(likelihoods).sum()
    return -log_likelihood

Tip

The phase space (phsp) is used for the normalization calculation for the unbinned log likelihood function

Optimizer: Minuit#

As a final step, we find the parameter values that minimize the unbinned Log likelihood function, i.e., we optimize the model with respect to the data distribution. In Python, this can be performed by e.g. iminuit package.

arg_names = tuple(guessed_parameters)
progress_bar = None


def wrapped_estimator(*args: float) -> float:
    # minuit migrad expects positional-argument functions
    global progress_bar
    if progress_bar is None:
        progress_bar = tqdm(desc="Running fit")
    new_kwargs = {
        **toy_parameters,
        **dict(zip(arg_names, args, strict=True)),
    }
    estimator_value = unbinned_nll(**new_kwargs)
    progress_bar.set_postfix({"estimator": estimator_value})
    progress_bar.update()
    return estimator_value

optimizer = Minuit(
    wrapped_estimator,
    **guessed_parameters,
    name=guessed_parameters,
)
optimizer.errordef = Minuit.LIKELIHOOD

optimizer.migrad()
progress_bar.close()
progress_bar = None  # reset progress bar

Final fit result#

	Parameter	Toy value	Starting value	Optimized value	Difference
0	a_minus2	0.00	NaN	NaN	0.0%
1	a_minus1	0.50	0.51	0.518316	3.7%
2	a_0	3.50	3.49	3.643171	4.1%
3	a_plus1	4.00	4.08	4.172394	4.3%
4	a_plus2	2.50	2.60	2.671472	6.9%
5	b_minus1	-1.50	-1.20	-1.600444	-6.7%
6	b_0	4.00	4.04	4.120737	3.0%
7	b_plus1	0.50	0.41	0.578660	15.7%
8	c_0	2.50	2.66	2.592561	3.7%
9	M12	1.32	1.40	1.320359	0.0%
10	Gamma12	0.10	0.15	0.100516	0.5%
11	M23	1.54	1.59	1.539956	0.0%
12	Gamma23	0.10	0.05	0.098987	1.0%
13	M31	1.87	1.82	1.869946	0.0%
14	Gamma31	0.10	0.20	0.100307	0.3%

Now we see the fit results of CM angles, Helicity angles and invariant mass as a comparison of data in 1D projection plots

Dalitz plots of phsp, data, and fit result#

To conclude our analysis in this document, we examine the Dalitz plots for the following cases: phase space, generated data, default parameters of the model, and fitted parameters of the model.

Phase Space Dalitz Plot:
- The plot displays a flat distribution, representing the kinematically allowed region for the reaction. This flat distribution serves as a baseline, showing no dynamical effects or resonant structures.
Generated Data Dalitz Plot:
- This plot is derived from synthetic data that simulates the actual experimental conditions. It includes potential resonances (and background processes). The generated data reflects the physical phenomena expected in the real experimental scenario.
Model with Default Parameters Dalitz Plot:
- Here, the Dalitz plot illustrates the model’s predictions using the initial set of parameters. This provides a theoretical reference for comparison with actual data and the fitted model.
Fitted Parameters Dalitz Plot:
- After performing the fitting procedure, this plot shows the model’s predictions with optimized parameters. The resonances and their characteristics are clearly visible, indicating the physical information extracted from the data.

By comparing these Dalitz plots, we can observe how the model evolves from theoretical predictions to a fitted solution that aligns with the experimental or generated data. The physical information about possible resonances becomes more apparent, allowing us to extract meaningful insights into the underlying processes of the reaction.

This comprehensive analysis showcases the interplay between kinematics and dynamics in understanding particle interactions, highlighting the power of Partial Wave Analysis in uncovering the nuances of hadronic reactions.

Summary#

In this document, we have covered the general workflow of PWA, which consists of three major parts:

As a summary, what have we done in this document consists of three major parts:

Amplitude Model formulation
- Introduced the theoretical model for the amplitude of the reaction.
- Defined the mathematical expressions representing the partial waves and their respective contributions.
- Incorporated the relevant physical parameters, such as resonance properties.
Phase space and data sample generation (or data analysis)
- Generated phase space distributions for the three-body decay processes.
- Created synthetic data samples or analyzed experimental data.
- Ensured that the data accurately reflects the kinematic constraints and dynamics of the reaction.
Perform fitting, analyze and evaluate the results
- Performed a fit of the amplitude model to the generated or experimental data.
- Used statistical techniques to optimize the model parameters.
- Analyzed the fit results to extract physical parameters.
- Evaluated the quality of the fit and the consistency of the model with the data.
- Interpreted the physical significance of the extracted parameters and identified potential resonances and their properties.

This structured approach provides a comprehensive understanding of both the reaction dynamics and kinematics, and helps in extracting meaningful physical insights from the data.

Tip

Extra have a look for another example of PWA for \(J/\psi \to \gamma \pi^0 \pi^0\) by using ComPWA at here .