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
import matplotlib.pyplot as plt
import numpy as np
import phasespace
import tensorflow as tf
from pylorentz import Momentum4
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)
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([[ -929.08221581, 219.7790733 , 2434.34426185, 2618.58901417],
[ 485.02418198, 2554.7278664 , -275.03756384, 2618.58901417],
[-1822.52343773, 1710.19628907, -768.87291651, 2618.58901417],
...,
[-1706.03641318, 917.30598422, -1756.5642824 , 2618.58901417],
[ 1536.69102564, -2097.02618315, 280.33991668, 2618.58901417],
[-2610.75054005, -136.17195971, 54.47389274, 2618.58901417]])>,
'p_1': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
array([[ 929.08221581, -219.7790733 , -2434.34426185, 2661.06098583],
[ -485.02418198, -2554.7278664 , 275.03756384, 2661.06098583],
[ 1822.52343773, -1710.19628907, 768.87291651, 2661.06098583],
...,
[ 1706.03641318, -917.30598422, 1756.5642824 , 2661.06098583],
[-1536.69102564, 2097.02618315, -280.33991668, 2661.06098583],
[ 2610.75054005, 136.17195971, -54.47389274, 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
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
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([[ 1672.3757309 , -259.34077472, 1925.96514886, 2715.77793272],
[ 2316.89032833, 1093.53593912, -98.17501568, 2715.77793272],
[-2180.9908383 , -240.25285876, 1326.2720326 , 2715.77793272],
...,
[-2136.77175711, -1402.44756565, -201.95856254, 2715.77793272],
[ -235.68083642, -102.18384633, 2550.97098813, 2715.77793272],
[ 182.94239133, 2496.74733071, -553.37584463, 2715.77793272]])>,
'gamma': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
array([[-1672.3757309 , 259.34077472, -1925.96514886, 2563.87206728],
[-2316.89032833, -1093.53593912, 98.17501568, 2563.87206728],
[ 2180.9908383 , 240.25285876, -1326.2720326 , 2563.87206728],
...,
[ 2136.77175711, 1402.44756565, 201.95856254, 2563.87206728],
[ 235.68083642, 102.18384633, -2550.97098813, 2563.87206728],
[ -182.94239133, -2496.74733071, 553.37584463, 2563.87206728]])>,
'K+': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
array([[ 1.43885543e+03, -3.56979640e+02, 1.32406299e+03,
2.04807206e+03],
[ 1.99226875e+03, 6.88082671e+02, -1.77873456e+02,
2.17208390e+03],
[-1.05779565e+03, 1.63192505e-02, 9.30828371e+02,
1.49301376e+03],
...,
[-1.97970190e+03, -1.07237045e+03, -2.62366923e+02,
2.31986020e+03],
[-3.08550502e+02, 6.19817898e+01, 2.31556348e+03,
2.38842969e+03],
[ 2.55478715e+02, 1.09131233e+03, -4.22385716e+02,
1.29551482e+03]])>,
'pi-': <tf.Tensor: shape=(100000, 4), dtype=float64, numpy=
array([[ 233.52030328, 97.63886521, 601.90216088, 667.70586822],
[ 324.62157833, 405.45326771, 79.69844041, 543.69402973],
[-1123.19518758, -240.26917801, 395.44366163, 1222.76417676],
...,
[ -157.06986182, -330.07711634, 60.40836056, 395.91773059],
[ 72.86966587, -164.16563612, 235.40750618, 327.34824483],
[ -72.53632379, 1405.43500357, -130.99012912, 1420.26311527]])>}
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
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
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()
