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, 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).
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\) with a fixed mass, we would write:
from phasespace import Particle
B0_MASS = 5279.58
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', 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 expressed as arrays of dimension (n_events, 4).
N_EVENTS = 1000
weights, particles = bz.generate(n_events=N_EVENTS)
The generate
method directly produces numpy arrays; for advanced usage, generate_tensor
returns the same objects with the
numpy arrays replaced by tf.Tensor
of the same shape.
So one can do, equivalent to the previous example:
import tensorflow as tf
with tf.Session() as sess:
weights, particles = sess.run(bz.generate_tensor(n_events=N_EVENTS))
In both cases, 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 both
generate
and generate_tensor
. 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).
Additionally, it is possible to obtain the unnormalized weights by using the generate_unnormalized
flag in
generate
and generate_tensor
. In this case, the method returns the unnormalized weights, the per-event maximum weight
and the particle dictionary.
>>> particles
{'K*': array([[ 1732.79325872, -1632.88873127, 950.85807735, 2715.78804872],
[-1633.95329448, 239.88921123, -1961.0402768 , 2715.78804872],
[ 407.15613764, -2236.6569286 , -1185.16616251, 2715.78804872],
...,
[ 1091.64603395, -1301.78721269, 1920.07503991, 2715.78804872],
[ -517.3125083 , 1901.39296899, 1640.15905194, 2715.78804872],
[ 656.56413668, -804.76922982, 2343.99214816, 2715.78804872]]),
'K+': array([[ 750.08077976, -547.22569019, 224.6920906 , 1075.30490935],
[-1499.90049089, 289.19714633, -1935.27960292, 2514.43047106],
[ 97.64746732, -1236.68112923, -381.09526192, 1388.47607911],
...,
[ 508.66157459, -917.93523639, 1474.7064148 , 1876.11771642],
[ -212.28646168, 540.26381432, 610.86656669, 976.63988936],
[ 177.16656666, -535.98777569, 946.12636904, 1207.28744488]]),
'gamma': array([[-1732.79325872, 1632.88873127, -950.85807735, 2563.79195128],
[ 1633.95329448, -239.88921123, 1961.0402768 , 2563.79195128],
[ -407.15613764, 2236.6569286 , 1185.16616251, 2563.79195128],
...,
[-1091.64603395, 1301.78721269, -1920.07503991, 2563.79195128],
[ 517.3125083 , -1901.39296899, -1640.15905194, 2563.79195128],
[ -656.56413668, 804.76922982, -2343.99214816, 2563.79195128]]),
'pi+': array([[ 982.71247896, -1085.66304109, 726.16598675, 1640.48313937],
[ -134.0528036 , -49.3079351 , -25.76067389, 201.35757766],
[ 309.50867032, -999.97579937, -804.0709006 , 1327.31196961],
...,
[ 582.98445936, -383.85197629, 445.36862511, 839.6703323 ],
[ -305.02604662, 1361.12915468, 1029.29248526, 1739.14815935],
[ 479.39757002, -268.78145413, 1397.86577911, 1508.50060384]])}
It is worth noting that the graph generation is cached even when using generate
, so iterative generation
can be performed using normal python loops without loss in performance:
for i in range(10):
weights, particles = bz.generate(n_events=1000)
...
(do something with weights and particles)
...
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 a tf.Tensor with the generated masses.
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
):
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):
ones = tf.ones((n_events, ), dtype=tf.float64)
kstar_mass = KSTARZ_MASS * ones
return tfp.distributions.TruncatedNormal(loc=kstar_mass,
scale=ones * KSTARZ_WIDTH,
low=min_mass,
high=max_mass).sample()
bz = Particle('B0', B0_MASS).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 takes
- The mass of the top particle.
- The mass of children particles as a list.
- The number of events to generate.
- The optional
boost_to
tensor.
For example, to generate \(B^0\to K\pi\), one would do:
import phasespace
N_EVENTS = 1000
B0_MASS = 5279.58
PION_MASS = 139.57018
KAON_MASS = 493.677
weights, particles = phasespace.generate(B0_MASS,
[PION_MASS, KAON_MASS],
n_events=N_EVENTS)
Internally, this function builds a decay chain using Particle
, and therefore the same considerations as before apply.
To avoid running the TensorFlow graph, one can set the as_numpy
flag to False
to get the graphs instead of the
numpy arrays.