Phase Space (hazma.phase_space
)¶
Overview¶
The phase_space
module contains the hazma.phase_space.Rambo
class for working with N-body phase-space and the specialized
hazma.phase_space.ThreeBody
class for three-body phase-space.
Rambo
contains methods for generating phase-space
points, integrating over N-body phase-space, computing cross-sections and
computing decay widths.
Rambo¶
The Rambo
class is instantiated by supplying the
center-of-mass energy, the masses of the final-state particles and optionally a
function to compute the squared matrix element. The function to compute the
squared matrix element must be a vectorized unary function that return the
squared matrix elements given four-momenta. The internal momenta generated by
Rambo
have a shape of (4, nfsp, n)
where nfsp in the number of final-state particles and n is the
number of points generated simultaneously. The leading dimension holds the
energy, x, y and z-components of the four-momenta. If no function to compute
the squared matrix element is supplied, it will be taken to be 1.
In the following snippet, we create a Rambo
object
for a process with 3 final-state particles, and a squared matrix element equal
to \(p_{1}\cdot p_{2}\). Note the Utilities (hazma.utils) module contains
a function hazma.utils.ldot()
to compute the Lorentzian scalar product
between two four-vectors.
from hazma import phase_space
from hazma import utils
cme = 10.0
masses = [1.0, 2.0, 3.0]
def msqrd(momenta):
p1 = momenta[:, 0] # Pick out the 4-momenta of particle 1
p2 = momenta[:, 1] # Pick out the 4-momenta of particle 2
return utils.ldot(p1, p2) # Compute p1.p2
phase_space = phase_space.Rambo(cme, masses, msqrd)
Integrating over phase-space¶
The Rambo.integrate()
method computes the following
integral:
This function takes in an integer specifying the number of points to use for the Monte-Carlo integration and computes the integral using the RAMBO algorithm. As a simple example, let’s compute the N-body integral for massless particles. The analytical result is:
We can verify this using the Rambo.integrate()
method:
import math
import numpy as np
from hazma import phase_space
def analytic(n):
fact = math.factorial(n - 2) * math.factorial(n - 1)
return (0.5 * math.pi) ** (n - 1) * (2 * np.pi) ** (4 - 3 * n) / fact
analytic_integrals = [analytic(n) for n in range(2, 10)]
integrals = []
for n in range(2, 10):
rambo = phase_space.Rambo(1.0, [0.0] * n)
integral, error = rambo.integrate(n=10)
integrals.append(integral)
np.max([abs(i1 - i2) / i2 for i1, i2 in zip(integrals, analytic_integrals)])
# Possible output: 2.2608966769988504e-15
Generating phase-space points¶
Sometimes it is useful to have access the momenta and weights of N-body
phase-space. The are two methods for generating momenta and weights:
Rambo.generate()
and func:Rambo.generator. The func:Rambo.generate
:py:method will return a specified number of momenta and weights. The
Rambo.generator()
method returns a python generator, allowing for
iterating over batches of momenta and weights.
As described above, the momenta will have a shape of (4, nfsp, n) where nfsp is the number of final-state particles and n is the number of requested points. The weights will have a shape of (n,). To see this explicitly, consider the following. Here we generate 10 phase-space points:
from hazma import phase_space
rambo = phase_space.Rambo(cme=10.0, masses=[1.0, 2.0, 3.0])
momenta, weights = rambo.generate(n=10)
print(momenta.shape)
print(weights.shape)
# (4, 3, 10)
# (10,)
In some cases, one may not want to generate all points at once (since it is
costly memory-wise or maybe one wants to monitor convergence.) We supply the
Rambo.generator()
method for generating batches of
phase-space points. For example, suppose we want to integrate over phase-space
ourselves using 1M points and batches of 50,000 points at a time. To do this,
we can use:
import numpy as np
from hazma import phase_space
rambo = phase_space.Rambo(cme=10.0, masses=[1.0, 2.0, 3.0])
n = int(1e6)
batch_size = 50_000
integrals = []
for momenta, weights in rambo.generator(n, batch_size, seed=1234):
integrals.append(np.nanmean(weights))
np.average(integrals)
# Output: 0.0036118278252665406
Note
The above code is equivalent to using
hazma.phase_space.Rambo.generate()
with n=int(1e6)
and
batch_size=50_000
. All methods accept a batch_size
argument used can
to split the computation into chunks in cases where the user has limited memory.
Computing decay widths and cross sections¶
The most common use of the Rambo()
is computing cross sections or decay
widths. The functions hazma.phase_space.Rambo.cross_section()
and
hazma.phase_space.Rambo.decay_width()
can be used for these purposes.
These methods just wrap hazma.phase_space.Rambo.integrate()
and append
the appropriate prefactors. As an example, let’s compute the muon decay width
for the process \(\mu\to e\nu_{e}\nu_{\mu}\). The squared matrix element
(ignoring the electron mass) is:
where \(G_{F}\) is the Fermi constant and \(t=(p_{e}+p_{\nu_{\mu}})^2\). The analytic result is
To compute this, we use the following:
from hazma import phase_space
from hazma import utils
from hazma.parameters import GF
from hazma.parameters import muon_mass as MMU
def msqrd(momenta):
p1 = momenta[:, 0]
p3 = momenta[:, 2]
t = utils.lnorm_sqr(p1 + p3)
return 16.0 * GF**2 * t * (MMU**2 - t)
rambo = phase_space.Rambo(MMU, [0.0, 0.0, 0.0], msqrd=msqrd)
width, error = rambo.decay_width(n=50_000, seed=1234)
analytic = GF**2 * MMU**5 / (192 * np.pi**3)
actual_error = abs(width - analytic)
print(f"width = {width:.2e} +- {error:.2e}")
print(f"actual error = {actual_error:.2e} = {actual_error / analytic * 100:.2f} %")
# Output:
# width = 3.02e-19 +- 5.99e-22
# actual error = 9.50e-22 = 0.32 %
Note
The Utilities (hazma.utils) module contains a couple of functions useful for dealing
with four-vectors, namely, hazma.utils.ldot()
for computing scalar
products and hazma.utils.lnorm_sqr()
for computing the squared
norm.
Energy and Invariant Mass Distributions¶
The Rambo
class can also compute energy distributions as well as the
invariant mass distributions of pairs of final state particles. To compute
energies distributions, use Rambo.energy_distributions
:
from hazma import phase_space
from hazma.parameters import standard_model_masses as sm_masses
import matplotlib.pyplot as plt
states = ["pi", "e", "mu", "k"]
masses = [sm_masses[s] for s in states]
cme = 3 * sum(masses)
rambo = phase_space.Rambo(cme, masses)
energy_dists = rambo.energy_distributions(n=1<<16, nbins=25)
plt.figure(dpi=150)
labels=[r"$\pi$", r"$e$", r"$\mu$", r"$K$"]
for i, dist in enumerate(energy_dists):
plt.plot(dist.bin_centers, dist.probabilities, label=labels[i])
plt.ylabel(r"$P(\epsilon) \ [\mathrm{MeV}^{-1}]$", fontdict=dict(size=16))
plt.xlabel(r"$\epsilon \ [\mathrm{MeV}]$", fontdict=dict(size=16))
plt.tight_layout()
plt.legend()
We note the ‘choppy’ behavior of the curves. Since we are performing naive
Monte-Carlo integration, we need a large number of points to properly sample
phase-space. The invariant mass distributions are generated in a similar
fashion. Note that the return value of the invariant mass distributions is a
dictionary with keys given by a pair of integers that specify the pair of
particles the distribution corresponds to. See
hazma.phase_space.PhaseSpaceDistribution1D
for more information on
the distribution objects.
Three Body Phase Space (hazma.phase_space.ThreeBody
)¶
Attention
This class is meant for cases where the squared matrix element only depends on the total momentum and the momenta of the final state particles. If this isn’t the case, then this class is not suitable.
As an example of a case where this class does apply, recall that the squared matrix element (summed over spins) for the decay of an unstable particle into a three-body final state only depends on the momenta of the final-state particles.
For three-body phase space where the squared matrix element only depends on the final-state momenta, the phase space integral simplifies. The result is:
where \(Q\) is the center-of-mass energy, \(s=(p_2+p_3)^2\) and \(t=(p_1+p_3)^2\). The limits of integration are:
Here, \(\lambda(a,b,c)=a^2+b^2+c^2-2ab-2ac-2bc\).
Integrating over phase-space¶
The interface for ThreeBody
is similar to Rambo
. To
integrate over phase space, use ThreeBody.integrate()
. The is an
important difference between ThreeBody
and Rambo
: the
matrix element must be a binary function taking the variables
\(s=(p_2+p_3)^2\) and \(t=(p_1+p_3)^2\).
As an example, let’s consider the muon decay again. To integrate over phase space, we use:
import numpy as np
from hazma import phase_space
from hazma.parameters import GF
from hazma.parameters import muon_mass as MMU
def msqrd(s, t):
return 16.0 * GF**2 * t * (MMU**2 - t)
integral, error = phase_space.ThreeBody(MMU, [0.0, 0.0, 0.0], msqrd=msqrd).integrate()
width = integral / (2.0 * MMU)
error = error / (2.0 * MMU)
analytic = GF**2 * MMU**5 / (192 * np.pi**3)
actual_error = abs(width - analytic)
print(f"analytic = {analytic:.8e}")
print(f"width = {width:.8e} +- {error:.2e}")
print(f"actual error = {actual_error:.2e} = {actual_error / analytic * 100:.2e} %")
# Output:
# analytic = 3.00917842e-16
# width = 3.00917842e-16 +- 6.68e-30
# actual error = 4.93e-32 = 1.64e-14 %
Energy and Invariant Mass Distributions¶
The interface for computing energy and invariant-mass distributions using
ThreeBody
is nearly identical to Rambo
. You can simply
construct your ThreeBody
instance and call either
ThreeBody.energy_distributions
or
ThreeBody.invariant_mass_distributions
. The main difference is the
signature of the msqrd
parameter. As before, it must be a binary function
accepting \(s=(p_2+p_3)^2\) and \(t=(p_1+p_3)^2\).
The hazma.form_factors
module uses ThreeBody
for the 3-body form
factors. As an example of how to generate distributions, let’s use the
hazma.form_factors.vector.VectorFormFactorPiKK0
to look at the energy
distributions of the final state mesons.
import numpy as np
import matplotlib.pyplot as plt
from hazma import phase_space
import hazma.form_factors.vector as ffv
def msqrd(s, t, q, form_factor: ffv.VectorFormFactorPiKK0, gvuu, gvdd, gvss):
ff = form_factor.form_factor(q, s, t, gvuu=gvuu, gvdd=gvdd, gvss=gvss)
lor = form_factor.squared_lorentz_structure(q, s, t)
return np.abs(ff) ** 2 * lor
# Specify parameters
gvuu, gvdd, gvss = 2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0
pk0_ff = ffv.VectorFormFactorPiKK0()
q = 1.5e3 # 1.5 GeV
# Construct `ThreeBody` object
tb = phase_space.ThreeBody(
q,
pk0_ff.fsp_masses,
msqrd=lambda s, t: msqrd(s, t, q, pk0_ff, gvuu, gvdd, gvss),
)
# Construct distributions
energy_dists = tb.energy_distributions(nbins=100)
inv_mass_dists = tb.invariant_mass_distributions(nbins=100)
# Make plot
plt.figure(dpi=150, figsize=(9, 3))
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)
labels = [r"$K^0$", r"$K$", r"$\pi$"]
lss = ["-", "--", "-."]
for i, dist in enumerate(energy_dists):
ax1.plot(dist.bin_centers, dist.probabilities, label=labels[i], ls=lss[i])
ax1.set_ylabel(r"$P(\epsilon)$", fontsize=16)
ax1.set_xlabel(r"$\epsilon$", fontsize=16)
labels = {
(0, 1): r"$(K^0,K)$",
(0, 2): r"$(K^0,\pi)$",
(1, 2): r"$(K,\pi)$",
}
lss = {(0, 1): "-", (0, 2): "--", (1, 2): "-."}
for key, dist in inv_mass_dists.items():
ax2.plot(dist.bin_centers, dist.probabilities, label=labels[key], ls=lss[key])
ax2.set_ylabel(r"$P(s_{ij})$", fontsize=16)
ax2.set_xlabel(r"$s_{ij}$", fontsize=16)
ax1.legend()
ax2.legend()
plt.tight_layout()
plt.show()