Usage¶
The base of phasespace
is the GenParticle
object.
This object, which represents a particle, either stable or decaying, has only one mandatory argument, its name.
In most cases (except for the top particle of a decay), one wants to also specify its mass, which can be either be array-like, or a function. Functions are used to specify the mass of particles such as resonances, which are not fixed but vary according to a broad distribution. These mass functions get three arguments, and must return an array-like object of the same shape as the input arguments.
The minimum mass allowed by the decay chain, which will be of shape (n_events,).
The maximum mass available, which will be of shape (n_events,).
The number of events to generate.
This function signature allows to handle threshold effects cleanly, giving enough information to produce kinematically
allowed decays (NB: phasespace
will throw an error if a kinematically forbidden decay is requested).
A simple example¶
With these considerations in mind, one can build a decay chain by using the set_children
method of the :py:class:~`phasespace.GenParticle`
class. As an example, to build the \(B^{0}\to K^{*}\gamma\) decay in which
\(K^*\to K\pi\) with a fixed mass, one would write:
from phasespace import GenParticle
import numpy as np
from particle import literals as lp
import vector
B0_MASS = lp.B_0.mass
KSTARZ_MASS = lp.Kst_892_0.mass
PION_MASS = lp.pi_plus.mass
KAON_MASS = lp.K_plus.mass
pion = GenParticle('pi-', PION_MASS)
kaon = GenParticle('K+', KAON_MASS)
kstar = GenParticle('K*', KSTARZ_MASS).set_children(pion, kaon)
gamma = GenParticle('gamma', 0)
bz = GenParticle('B0', B0_MASS).set_children(kstar, gamma)
Phasespace events can be generated using the generate
method, which gets the number of events to generate as input.
The method returns:
The normalized weights of each event, as an array of dimension (n_events,).
The 4-momenta of the generated particles as values of a dictionary with the particle name as key. These momenta are either expressed as arrays of dimension (n_events, 4) or :py:class:~`vector.Vector4D` objects, depending on the
as_vectors
flag given togenerate
.
N_EVENTS = 1000
weights, particles = bz.generate(n_events=N_EVENTS, as_vectors=True)
# or
weights, particles = bz.generate(n_events=N_EVENTS)
(The array-like objects can always directly be converted to a numpy array (if really needed) through np.asarray(obj)
.)
Boosting the particles¶
The particles are generated in the rest frame of the top particle.
To produce them at a given momentum of the top particle, one can pass these momenta with the boost_to
argument in
generate
. This latter approach can be useful if the momentum of the top particle
is generated according to some distribution, for example the kinematics of the LHC (see test_kstargamma_kstarnonresonant_lhc
and test_k1gamma_kstarnonresonant_lhc
in tests/test_physics.py
to see how this could be done).
The boost_to
argument can be a 4-momentum array of shape (n_events, 4) with (px, py, yz, energy) or a :py:class:~`vector.Vector4D` (both a momentum and a Lorentz vector).
N_EVENTS = 1000
# Generate the top particle with a momentum of 100 GeV
top_momentum = np.array([0, 0, 100, np.sqrt(100**2 + B0_MASS**2)])
# or
top_momentum = vector.array({'px': [0], 'py': [0], 'pz': [100], 'E': [np.sqrt(100**2 + B0_MASS**2)]})
weights, particles = bz.generate(n_events=N_EVENTS, boost_to=top_momentum)
Weights¶
Additionally, it is possible to obtain the unnormalized weights by using the generate_unnormalized
flag in
generate
. In this case, the method returns the unnormalized weights, the per-event maximum weight
and the particle dictionary.
print(particles)
{'K*': <tf.Tensor: shape=(1000, 4), dtype=float64, numpy=
array([[ 692.24719713, -1312.41100297, 2142.70756462, 2755.87269891],
[ 1275.42679112, -2192.33502973, 426.23233161, 2723.36743483],
[ -22.29634617, -1958.69487113, 1705.96732696, 2747.60205423],
...,
[ -42.27486169, -1435.14629279, 2175.97453851, 2756.50268282],
[-2552.74098348, 234.5511906 , 95.78122635, 2717.10961052],
[ -846.55711145, 1931.681667 , -1406.70263309, 2688.65674752]])>, 'gamma': <tf.Tensor: shape=(1000, 4), dtype=float64, numpy=
array([[ -692.24719713, 1312.41100297, -2042.70756462, 2524.73424685],
[-1275.42679112, 2192.33502973, -326.23233161, 2557.23951093],
[ 22.29634617, 1958.69487113, -1605.96732696, 2533.00489153],
...,
[ 42.27486169, 1435.14629279, -2075.97453851, 2524.10426294],
[ 2552.74098348, -234.5511906 , 4.21877365, 2563.49733524],
[ 846.55711145, -1931.681667 , 1506.70263309, 2591.95019823]])>, 'pi-': <tf.Tensor: shape=(1000, 4), dtype=float64, numpy=
array([[ 275.7024959 , -149.37122026, 776.30675493, 848.79661828],
[ 762.33680325, -781.51782636, 66.6050566 , 1102.65300155],
[ 248.66784768, -576.76596603, 683.23186583, 938.49899003],
...,
[ -80.94259439, -844.986719 , 775.90579069, 1158.47483731],
[-1382.72652864, 349.29357275, 150.52671913, 1440.8597572 ],
[ -406.15527628, 328.11724388, -203.32162091, 577.44489735]])>, 'K+': <tf.Tensor: shape=(1000, 4), dtype=float64, numpy=
array([[ 416.54470123, -1163.03978271, 1366.40080969, 1907.07608063],
[ 513.08998788, -1410.81720337, 359.627275 , 1620.71443328],
[ -270.96419385, -1381.9289051 , 1022.73546113, 1809.1030642 ],
...,
[ 38.6677327 , -590.15957379, 1400.06874782, 1598.02784551],
[-1170.01445484, -114.74238215, -54.74549278, 1276.24985332],
[ -440.40183517, 1603.56442312, -1203.38101217, 2111.21185017]])>}
Iterative generation can be performed using normal python loops without loss in performance:
for i in range(5):
weights, particles = bz.generate(n_events=100)
# ...
# (do something with weights and particles)
# ...
Resonances with variable mass¶
To generate the mass of a resonance, we need to give a function as its mass instead of a floating number. This function should take as input the per-event lower mass allowed, per-event upper mass allowed and the number of events, and should return an array-like object with the generated masses and shape (nevents,). Well suited for this task are the TensorFlow Probability distributions or, for more customized mass shapes, the zfit pdfs (currently an experimental feature is needed, contact the zfit developers to learn more).
Following with the same example as above, and approximating the resonance shape by a gaussian, we could
write the \(B^{0}\to K^{*}\gamma\) decay chain as (more details can be found in tests/helpers/decays.py
):
from phasespace import numpy as tnp
import tensorflow_probability as tfp
from phasespace import GenParticle
KSTARZ_MASS = 895.81
KSTARZ_WIDTH = 47.4
def kstar_mass(min_mass, max_mass, n_events):
min_mass = tnp.asarray(min_mass, tnp.float64)
max_mass = tnp.asarray(max_mass, tnp.float64)
kstar_width_cast = tnp.asarray(KSTARZ_WIDTH, tnp.float64)
kstar_mass_cast = tnp.asarray(KSTARZ_MASS, tnp.float64)
kstar_mass = tnp.broadcast_to(kstar_mass_cast, shape=(n_events,))
if KSTARZ_WIDTH > 0:
kstar_mass = tfp.distributions.TruncatedNormal(loc=kstar_mass,
scale=kstar_width_cast,
low=min_mass,
high=max_mass).sample()
return kstar_mass
bz = GenParticle('B0', B0_MASS).set_children(GenParticle('K*0', mass=kstar_mass)
.set_children(GenParticle('K+', mass=KAON_MASS),
GenParticle('pi-', mass=PION_MASS)),
GenParticle('gamma', mass=0.0))
bz.generate(n_events=500)
Shortcut for simple decays¶
The generation of simple n-body decay chains can be done using the nbody_decay
function of phasespace
, which takes
The mass of the top particle.
The mass of children particles as a list.
The name of the top particle (optional).
The names of the children particles (optional).
If the names are not given, top and p_{i} are assigned. For example, to generate \(B^0\to K\pi\), one would do:
import phasespace
from particle import literals as lp
N_EVENTS = 1000
B0_MASS = lp.B_0.mass
PION_MASS = lp.pi_plus.mass
KAON_MASS = lp.K_plus.mass
decay = phasespace.nbody_decay(B0_MASS, [PION_MASS, KAON_MASS],
top_name="B0", names=["pi", "K"])
weights, particles = decay.generate(n_events=N_EVENTS)
In this example, decay
is simply a GenParticle
with the corresponding children.
Eager execution¶
By default, phasespace uses JIT (just-in-time) compilation of TensorFlow to greatly speed up the generation of events. Simplified, this means that the first time a decay is generated, a symbolic array without a concrete value is used and the computation is remembered. As a user calling the function, you will not notice this, the output will be the same as if the function was executed eagerly. The consequence is two-fold: on one hand the initial overhead is higher with a significant speedup for subsequent generations, on the other hand, the values of the generated particles inside the function are not available in pure Python (e.g. for debugging basically).
If you need to debug the internals, using tf.config.run_functions_eagerly(True)
(or the environment variable "PHASESPACE_EAGER=1"
) will make everything run numpy-like.
Random numbers¶
The random number generation inside phasespace is transparent in order to allow for deterministic behavior if desired. A function that uses random number generation inside always takes a seed (or rng) argument. The behavior is as follows
if no seed is given, the global random number generator of TensorFlow will be used. Setting this instance explicitly or by setting the seed via tf.random.set_seed allows for a deterministic execution of a whole _script_.
if the seed is a number it will be used to create a random number generator from this. Using the same seed again will result in the same output.
if the seed is an instance of
tf.random.Generator
, this instance will directly be used and advances an undefined number of steps.