Usage¶
The base of phasespace is the Particle 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
a number or tf.constant, 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 a TensorFlow Tensor:
- The minimum mass allowed by the decay chain.
- The maximum mass available.
- The number of events to generate, i.e., the size of the output tensor, which must be of shape (1,
n_events).
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).
With these considerations in mind, one can build a decay chain by using the set_children method of the Particle
class (which returns the class itself). As an example, to build the \(B^{0}\to K^{*}\gamma\) decay in which
\(K^*\to K\pi\), we would write:
from phasespace import Particle
KSTARZ_MASS = 895.81
PION_MASS = 139.57018
KAON_MASS = 493.677
pion = Particle('pi+', PION_MASS)
kaon = Particle('K+', KAON_MASS)
kstar = Particle('K*', KSTARZ_MASS).set_children(pion, kaon)
gamma = Particle('gamma', 0)
bz = Particle('B0').set_children(kstar, gamma)
The graph can then be generated by using the generate method, which, similarly to TGenPhaseSpace, gets an input
4-vector (shape (4,)) of the top particle, and, optionally, the number of events to generate.
If the input 4-vector is of shape (4, n), then n events are generated.
Since generate just creates a TensorFlow graph, we need to use tf.Session.run() to actually run the generation:
N_EVENTS = 1000
B0_MASS = 5279.58
B0_AT_REST = [0.0, 0.0, 0.0, B0_MASS]
with tf.Session() as sess:
weights, particles = sess.run(bz.generate(B0_AT_REST, N_EVENTS))
As discussed, this would be equivalent to:
B0_MASS = 5279.58
B0s_AT_REST = [[0.0, 0.0, 0.0, B0_MASS] for _ in range(N_EVENTS)]
with tf.Session() as sess:
weights, particles = sess.run(bz.generate(B0s_AT_REST))
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 objects returned by Particle.generate are the normalized weights vector (the weight of each event divided by
the maximum size of the phase space given the decay chain), of shape (N_EVENTS, ), and a dictionary containing
the generated 4-momenta for each particle, of shape (4, N_EVENTS), labelled by their name:
>>> particles
{'K*': array([[-2259.88717495, 742.20158838, -1419.57804967, ...,
385.51632682, 890.89417859, -1938.80489221],
[ -491.3119786 , -2348.67021741, -2049.19459865, ...,
-932.58261761, -1054.16217965, -1669.40481126],
[-1106.5946257 , 711.27644522, -598.85626591, ...,
-2356.84025605, -2160.57372728, -164.77965753],
[ 2715.78804872, 2715.78804872, 2715.78804872, ...,
2715.78804872, 2715.78804872, 2715.78804872]]),
'K+': array([[-1918.74294565, 363.10302225, -830.13803095, ...,
9.28960349, 850.87382095, -895.29815921],
[ -566.15415012, -956.94044749, -1217.14751182, ...,
-243.52446264, -1095.04308712, -1078.03237584],
[-1108.26109897, 534.79579335, -652.41135612, ...,
-901.56453631, -2069.39723754, -244.1159568 ],
[ 2339.67191226, 1255.90698132, 1685.21060224, ...,
1056.37401241, 2539.53293518, 1505.66336806]]),
'gamma': array([[2259.88717495, -742.20158838, 1419.57804967, ..., -385.51632682,
-890.89417859, 1938.80489221],
[ 491.3119786 , 2348.67021741, 2049.19459865, ..., 932.58261761,
1054.16217965, 1669.40481126],
[1106.5946257 , -711.27644522, 598.85626591, ..., 2356.84025605,
2160.57372728, 164.77965753],
[2563.79195128, 2563.79195128, 2563.79195128, ..., 2563.79195128,
2563.79195128, 2563.79195128]]),
'pi+': array([[ -341.14422931, 379.09856613, -589.44001872, ...,
376.22672333, 40.02035764, -1043.506733 ],
[ 74.84217153, -1391.72976992, -832.04708683, ...,
-689.05815497, 40.88090746, -591.37243542],
[ 1.66647327, 176.48065186, 53.55509021, ...,
-1455.27571974, -91.17648974, 79.33629927],
[ 376.11613646, 1459.8810674 , 1030.57744648, ...,
1659.41403631, 176.25511354, 1210.12468065]])}
If interested in the unnormalized event weights, one can use the generate_unnormalized method, which
returns the raw weights, the per-event maximum weight and the particle dictionary as before.
To generate the mass of a resonance, we need to give a function as its mass. 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:
import tensorflow as tf
import tensorflow_probability as tfp
from phasespace import Particle
KSTARZ_MASS = 895.81
KSTARZ_WIDTH = 47.4
def kstar_mass(min_mass, max_mass, n_events):
kstar_mass = KSTARZ_MASS * ones
min_mass = tf.broadcast_to(min_mass, (1, n_events))
max_mass = tf.broadcast_to(max_mass, (1, n_events))
kstar_mass = tfp.distributions.TruncatedNormal(loc=KSTARZ_MASS,
scale=KSTARZ_WIDTH,
low=min_mass,
high=max_mass).sample()
bz = Particle('B0').set_children(Particle('K*0', mass=kstar_mass)
.set_children(Particle('K+', mass=KAON_MASS),
Particle('pi-', mass=PION_MASS)),
Particle('gamma', mass=0.0))
Shortcut for simple decays¶
The generation of simple n-body decays can be done using the generate function of phasespace, which has a
very similar interface to TGenPhaseSpace. For example, to generate \(B^0\to K\pi\), we would do:
import phasespace
import tensorflow as tf
N_EVENTS = 1000
B0_MASS = 5279.58
B0_AT_REST = [0.0, 0.0, 0.0, B0_MASS]
PION_MASS = 139.57018
KAON_MASS = 493.677
with tf.Session() as sess:
weights, particles = sess.run(phasespace.generate(B0_AT_REST,
[PION_MASS, KAON_MASS],
N_EVENTS))
In this case, since particles are unnamed, the particles object contains a list of (4, N_EVENTS) tensors in the order
of the particles specified in the generate call.
Internally, this function builds a decay chain using Particle, and therefore the same considerations as before apply;
for example, it is possible to not specify the number of events and give a list of input momenta.