Lecture 12 – Phase space simulation

How to perform a phase space simulation with Python?

Prepare the notebook by importing numpy, matplotlib.pyplot, and pylorentz and download the data file (CSV format) from Google Drive.

Import Python libraries

Hide code cell content

import matplotlib.pyplot as plt
import numpy as np
import phasespace
import tensorflow as tf
from pylorentz import Momentum4

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
E0000 00:00:1754936581.462080    2872 cuda_dnn.cc:8310] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
E0000 00:00:1754936581.466200    2872 cuda_blas.cc:1418] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered

Two-body decay

Let’s start with a simple two body decay at rest: \(B^0\rightarrow K^+\pi^-\).

B0_MASS = 5279.65
PION_MASS = 139.57018
KAON_MASS = 493.677
n_events = 100_000

decay = phasespace.nbody_decay(B0_MASS, [PION_MASS, KAON_MASS])
weights, four_momenta = decay.generate(n_events=n_events)

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1754936587.878944    2893 service.cc:148] XLA service 0x7fe81c005bd0 initialized for platform Host (this does not guarantee that XLA will be used). Devices:
I0000 00:00:1754936587.879010    2893 service.cc:156]   StreamExecutor device (0): Host, Default Version
I0000 00:00:1754936587.901172    2894 device_compiler.h:188] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

The simulation produces a dictionary (four_momenta) of tf.Tensor objects. Each object can be addressed with particles['p_i'], where i is the number of the \(i\)-th generated particle.

four_momenta

{'p_0': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[ 2150.03754434,   649.22069653, -1339.17124936,  2618.58901417],
        [ -806.95001273,   880.12921747,  2326.31314047,  2618.58901417],
        [  335.5975645 ,  2553.19756917,  -453.96589899,  2618.58901417],
        ...,
        [ 1496.44531762,  1552.22180007, -1479.45513103,  2618.58901417],
        [ 1532.19215393,  1070.62181806, -1828.57450384,  2618.58901417],
        [  225.95680655,  2564.10723349,   460.24581119,  2618.58901417]])>,
 'p_1': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[-2150.03754434,  -649.22069653,  1339.17124936,  2661.06098583],
        [  806.95001273,  -880.12921747, -2326.31314047,  2661.06098583],
        [ -335.5975645 , -2553.19756917,   453.96589899,  2661.06098583],
        ...,
        [-1496.44531762, -1552.22180007,  1479.45513103,  2661.06098583],
        [-1532.19215393, -1070.62181806,  1828.57450384,  2661.06098583],
        [ -225.95680655, -2564.10723349,  -460.24581119,  2661.06098583]])>}

Each tf.Tensor can be converted to a NumPy array, which can then be converted to a pylorentz.

def to_lorentz(p: tf.Tensor) -> Momentum4:
    p = p.numpy().T
    return Momentum4(p[3], *p[:3])

pion = to_lorentz(four_momenta["p_0"])
kaon = to_lorentz(four_momenta["p_1"])

These objects can be used to do kinematic computations. Let’s first verify that the invariant mass of the kaon+pion system corresponds to the mass of the mother \(B^0\):

B0 = pion + kaon
np.testing.assert_almost_equal(B0.m.mean(), B0_MASS)

Let’s also plot the momentum components of the two daugther particles.

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 4, tight_layout=True, figsize=(11, 3.5))
ax[0].hist(B0.m, bins=100, color="CornFlowerBlue", range=(5279, 5280))
ax[0].set_xlabel(R"i.m.($\pi$K) [MeV/$c^2$]")
ax[0].set_title("(pion-kaon) i.m. \n")

ax[1].hist(kaon.p_x, bins=70, color="lightcoral", hatch="//")
ax[1].hist(
    pion.p_x,
    bins=70,
    color="springgreen",
    histtype="barstacked",
    hatch="\\",
    alpha=0.5,
)
ax[1].set_xlabel("$p_x$ [MeV/$c$]")
ax[1].set_title("x mom. component \n")

ax[2].hist(kaon.p_y, bins=70, color="lightcoral", hatch="//")
ax[2].hist(
    pion.p_y,
    bins=70,
    color="springgreen",
    histtype="barstacked",
    hatch="\\",
    alpha=0.5,
)
ax[2].set_xlabel("$p_y$ [MeV/$c$]")
ax[2].set_title("y mom. component \n")

ax[3].hist(kaon.p_z, bins=70, color="lightcoral", hatch="//")
ax[3].hist(
    pion.p_z,
    bins=70,
    color="springgreen",
    histtype="barstacked",
    hatch="\\",
    alpha=0.5,
)
ax[3].set_xlabel("$p_z$ [MeV/$c$]")
ax[3].set_title("z mom. component \n")
plt.show()

But it’s monochromatic!! of course it is… it’s a decay at rest. The momentum components are uniformly distributed in the available phase space.

Three-body decay

Let’s consider now a three body decay like \(B^0\rightarrow K^+\pi^-\pi^0\) and repeat the plot of the relevant kinematic variables. We can also make Dalitz plots this time.

n_events = 50_000
PION0_MASS = 134.9766
decay = phasespace.nbody_decay(B0_MASS, [PION_MASS, PION0_MASS, KAON_MASS])
weights, four_momenta = decay.generate(n_events=n_events)

pim = to_lorentz(four_momenta["p_0"])
pi0 = to_lorentz(four_momenta["p_1"])
kaon = to_lorentz(four_momenta["p_2"])
s1 = (kaon + pim).m2
s2 = (kaon + pi0).m2
s3 = (pim + pi0).m2

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3.5))
f0 = ax[0].hist2d(s1, s3, bins=70, cmap="turbo", cmin=1)
fig.colorbar(f0[3], ax=ax[0])
ax[0].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[0].set_ylabel(R"i.m.$^2(\pi^-\pi^0)$ [(MeV/$c^2)^2]$")

f1 = ax[1].hist2d(s2, s3, bins=70, cmap="turbo", cmin=1)
fig.colorbar(f1[3], ax=ax[1])
ax[1].set_xlabel(R"i.m.$^2(\pi^0K^+)$ [(MeV/$c^2)^2]$")
ax[1].set_ylabel(R"i.m.$^2(\pi^-\pi^0)$ [(MeV/$c^2)^2]$")

f2 = ax[2].hist2d(s1, s2, bins=70, cmap="turbo", cmin=1)
fig.colorbar(f2[3], ax=ax[2])
ax[2].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[2].set_ylabel(R"i.m.$^2(\pi^0K^+)$ [(MeV/$c^2)^2]$")
plt.show()

Decay chain

The phasespace package allows to treat also multiple decays. Let’s consider the \(B^0\rightarrow K^{\ast 0}\gamma\) decay, followed by \(K^{\ast 0}\rightarrow \pi^-K^+\). It can be simulated using the following procedure:

from phasespace import GenParticle

B0_MASS = 5279.65
K0STAR_MASS = 895.55
PION_MASS = 139.57018
KAON_MASS = 493.677
GAMMA_MASS = 0.0

Kp = GenParticle("K+", KAON_MASS)
pim = GenParticle("pi-", PION_MASS)
Kstar = GenParticle("KStar", K0STAR_MASS).set_children(Kp, pim)
gamma = GenParticle("gamma", GAMMA_MASS)
B0 = GenParticle("B0", B0_MASS).set_children(Kstar, gamma)

weights, four_momenta = B0.generate(n_events=100_000)

four_momenta

{'KStar': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[ 1203.27865264,  1712.88496772, -1480.40046907,  2715.77793272],
        [ -517.7365893 ,  2506.83302917,   145.52307494,  2715.77793272],
        [ -950.14307847,  2183.30204899,   950.71566237,  2715.77793272],
        ...,
        [ -318.02301541, -2498.86075021,  -477.48957068,  2715.77793272],
        [ 2559.24620387,   153.67397429,    -9.11893428,  2715.77793272],
        [ 1601.93602851, -1420.02330457,  1410.94817495,  2715.77793272]])>,
 'gamma': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[-1203.27865264, -1712.88496772,  1480.40046907,  2563.87206728],
        [  517.7365893 , -2506.83302917,  -145.52307494,  2563.87206728],
        [  950.14307847, -2183.30204899,  -950.71566237,  2563.87206728],
        ...,
        [  318.02301541,  2498.86075021,   477.48957068,  2563.87206728],
        [-2559.24620387,  -153.67397429,     9.11893428,  2563.87206728],
        [-1601.93602851,  1420.02330457, -1410.94817495,  2563.87206728]])>,
 'K+': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[  654.66001818,   565.25505916,  -692.7926642 ,  1213.16597266],
        [ -212.61215218,  2067.97425862,   257.3916297 ,  2152.13589091],
        [-1002.5676664 ,  1876.93434178,   930.80844311,  2374.47800309],
        ...,
        [ -101.38414792, -2210.27642836,  -521.98857086,  2326.32536066],
        [ 1889.05013489,   379.96441923,   -84.00239424,  1990.89345636],
        [  916.37171002,  -481.43276486,   629.20975627,  1308.10416853]])>,
 'pi-': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
 array([[ 548.61863445, 1147.62990856, -787.60780487, 1502.61196006],
        [-305.12443711,  438.85877055, -111.86855476,  563.64204181],
        [  52.42458794,  306.36770721,   19.90721926,  341.29992963],
        ...,
        [-216.63886749, -288.58432185,   44.49900019,  389.45257206],
        [ 670.19606898, -226.29044495,   74.88345995,  724.88447636],
        [ 685.5643185 , -938.5905397 ,  781.73841868, 1407.6737642 ]])>}

gamma = to_lorentz(four_momenta["gamma"])
pion = to_lorentz(four_momenta["pi-"])
kaon = to_lorentz(four_momenta["K+"])
Kstar = to_lorentz(four_momenta["KStar"])

Let’s build the Dalitz plots matching particle pairs. The particles measured in the final state are \(K^-,\; \pi^-\) and \(\gamma\).

s1 = (pion + kaon).m2
s2 = (gamma + kaon).m2
s3 = (gamma + pion).m2

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3.5))

f0 = ax[0].hist2d(s1, s2, bins=70, cmap="turbo")
fig.colorbar(f0[3], ax=ax[0])
ax[0].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[0].set_ylabel(R"i.m.$^2(K^+\gamma)$ [(MeV/$c^2)^2]$")

f1 = ax[1].hist2d(s2, s3, bins=70, cmap="turbo")
fig.colorbar(f1[3], ax=ax[1])
ax[1].set_xlabel(R"i.m.$^2(K^+\gamma)$ [(MeV/$c^2)^2]$")
ax[1].set_ylabel(R"i.m.$^2(\pi^-\gamma)$ [(MeV/$c^2)^2]$")

f2 = ax[2].hist2d(s1, s3, bins=70, cmap="turbo")
fig.colorbar(f2[3], ax=ax[2])
ax[2].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[2].set_ylabel(R"i.m.$^2(\pi^-\gamma)$ [(MeV/$c^2)^2]$")
plt.show()

Width distribution

These distributions aren’t so interesting, because the masses of each particle are one fixed value. So let’s simulate a more realistic \(K^\ast\) particle; not monochromatic, but with a width of 47 MeV.[1] The mass is extracted from a Gaussian distribution centered at the B0_MASS value and with \(\sigma = 47/2.36 \sim 20\) MeV. See more info on how to do this with the phasespace package here.

import tensorflow as tf
import tensorflow_probability as tfp

K0STAR_WIDTH = 47 / 2.36


def kstar_mass(min_mass, max_mass, n_events):
    min_mass = tf.cast(min_mass, tf.float64)
    max_mass = tf.cast(max_mass, tf.float64)
    kstar_mass_cast = tf.cast(K0STAR_MASS, dtype=tf.float64)
    tf.cast(K0STAR_WIDTH, tf.float64)
    tf.broadcast_to(kstar_mass_cast, shape=(n_events,))
    return tfp.distributions.TruncatedNormal(
        loc=K0STAR_MASS,
        scale=K0STAR_WIDTH,
        low=min_mass,
        high=max_mass,
    ).sample()

K = GenParticle("K+", KAON_MASS)
pion = GenParticle("pi-", PION_MASS)
Kstar = GenParticle("KStar", kstar_mass).set_children(K, pion)
gamma = GenParticle("gamma", GAMMA_MASS)
B0 = GenParticle("B0", B0_MASS).set_children(Kstar, gamma)
weights, four_momenta = B0.generate(n_events=100_000)

gamma = to_lorentz(four_momenta["gamma"])
pion = to_lorentz(four_momenta["pi-"])
kaon = to_lorentz(four_momenta["K+"])
Kstar = to_lorentz(four_momenta["KStar"])

Now you have all the 4-vectors to plot the invariant mass distributions for the different steps of the decay chains.

s1 = (pion + kaon).m2
s2 = (gamma + kaon).m2
s3 = (gamma + pion).m2

Show code cell source

Hide code cell source

fig, ax = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3.5))
f0 = ax[0].hist2d(s1, s2, bins=70, cmap="rainbow")
fig.colorbar(f0[3], ax=ax[0])
ax[0].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[0].set_ylabel(R"i.m.$^2(K^+\gamma)$ [(MeV/$c^2)^2]$")

f1 = ax[1].hist2d(s2, s3, bins=70, cmap="rainbow")
fig.colorbar(f1[3], ax=ax[1])
ax[1].set_xlabel(R"i.m.$^2(K^+\gamma)$ [(MeV/$c^2)^2]$")
ax[1].set_ylabel(R"i.m.$^2(\pi^-\gamma)$ [(MeV/$c^2)^2]$")

f2 = ax[2].hist2d(s1, s3, bins=70, cmap="rainbow")
fig.colorbar(f2[3], ax=ax[2])
ax[2].set_xlabel(R"i.m.$^2(\pi^-K^+)$ [(MeV/$c^2)^2]$")
ax[2].set_ylabel(R"i.m.$^2(\pi^-\gamma)$ [(MeV/$c^2)^2]$")
plt.show()

---

[1]

Do you remember the difference between a BW width and a Lorentzian FWHM? The widths you get out of the PDG booklet are FWHM’s.