Step 2: Generate data samples¶

In this section, we will use the AmplitudeModel that we created with the expert system in the previous step to generate a data sample via hit & miss Monte Carlo. We do this with the tensorwaves.data.generate module.

First, we load a AmplitudeModel that was created in the previous step:

from expertsystem import io

model = io.load_amplitude_model("amplitude_model_helicity.yml")

2.1 Define kinematics of the particle reaction¶

An AmplitudeModel defines the kinematics, the particles involved in the reaction, the dynamics used for the model on which to perform the eventual optimization, etc.

from tensorwaves.physics.helicity_formalism.kinematics import HelicityKinematics

kin = HelicityKinematics.from_model(model)
print("Initial state mass:", kin.reaction_kinematics_info.initial_state_masses)
print("Final state masses:", kin.reaction_kinematics_info.final_state_masses)
print("Involved particles:", model.particles.names)

Initial state mass: [3.0969]
Final state masses: [0.0, 0.1349768, 0.1349768]
Involved particles: {'f(0)(980)', 'f(0)(1500)', 'f(0)(1710)', 'gamma', 'J/psi(1S)', 'f(0)(500)', 'pi0', 'f(0)(1370)'}

A Kinematics object defines the kinematics of the particle reaction we are studying. From the final state masses listed here, we can see we are dealing with the reaction \(J/\psi \to \gamma\pi^0\pi^0\).

2.2 Generate phase space sample¶

Kinematics define the constraints of the phase space. As such, we have enough information to generate a phase-space sample for this particle reaction. We do this with the generate_phsp function. By default, this function uses TFPhaseSpaceGenerator as a, well… phase-space generator (using tensorflow in the back-end) and generates random numbers with TFUniformRealNumberGenerator. You can specify this with the arguments of generate_phsp function.

from tensorwaves.data.generate import generate_phsp

phsp_sample = generate_phsp(300000, kin)
phsp_sample.shape

(3, 300000, 4)

As you can see, the phase-space sample is a three-dimensional array: 300.000 events of four-momentum tuples for three particles.

2.3 Generate intensity-based sample¶

‘Data samples’ are more complicated than phase space samples in that they represent the intensity profile resulting from a reaction. You therefore need an IntensityTF object (or, more generally, a Function instance) and a phase space over which to generate that intensity distribution. We call such a data sample an intensity-based sample.

An intensity-based sample is generated with the function generate_data. Its usage is similar to generate_phsp, but now you have to give an IntensityTF in addition to the Kinematics object. An IntensityTF object can be created with the IntensityTF class:

from tensorwaves.physics.helicity_formalism.amplitude import IntensityBuilder

builder = IntensityBuilder(model.particles, kin, phsp_sample)
intensity = builder.create_intensity(model)

That’s it, now we have enough info to create an intensity-based data sample. Notice how the structure is the sample as the phase-space sample we saw before:

from tensorwaves.data.generate import generate_data

data_sample = generate_data(30000, kin, intensity)
data_sample.shape

(3, 30000, 4)

2.4 Visualize kinematic variables¶

We now have a phase space sample and an intensity-based sample. Their data structure isn’t the most informative though: it’s just a collection of four-momentum tuples. However, we can use the Kinematics class to convert those four-momentum tuples to a data set of kinematic variables.

Now we can use the Kinematics.convert method to convert the phase space and data samples of four-momentum tuples to kinematic variables.

phsp_set = kin.convert(phsp_sample)
data_set = kin.convert(data_sample)
list(data_set)

['mSq_2_3_4',
 'mSq_3_4',
 'mSq_2',
 'mSq_3',
 'mSq_4',
 'theta+3_4+2',
 'phi+3_4+2',
 'theta+3+4_vs_2',
 'phi+3+4_vs_2']

The data set is just a dict of kinematic variables (keys are the names, values is a list of computed values for each event). The numbers you see here are final state IDs as defined in the AmplitudeModel member of the AmplitudeModel:

for state_id, particle in model.kinematics.final_state.items():
    print(f"ID {state_id}:", particle.name)

ID 2: gamma
ID 3: pi0
ID 4: pi0

Available kinematic variables

By default, tensorwaves only generates invariant masses of the SubSystems that are of relevance to the decay problem. In this case, we only have resonances \(f_0 \to \pi^0\pi^0\). If you are interested in more invariant mass combinations, you can do so with the method HelicityKinematics.register_invariant_mass(), e.g.:

kin.register_invariant_mass([2, 4])

The data set is dict, which allows us to easily convert it to a pandas.DataFrame and plot its content in the form of a histogram:

import pandas as pd

data_frame = pd.DataFrame(data_set)
phsp_frame = pd.DataFrame(phsp_set)
data_frame

	mSq_2_3_4	mSq_3_4	mSq_2	mSq_3	mSq_4	theta+3_4+2	phi+3_4+2	theta+3+4_vs_2	phi+3+4_vs_2
0	9.59079	0.456118	0.0	0.018219	0.018219	1.509823	2.029175	1.164715	-2.472122
1	9.59079	2.944665	0.0	0.018219	0.018219	0.630947	-1.732462	1.332378	-2.157104
2	9.59079	0.956728	0.0	0.018219	0.018219	1.235279	1.745086	1.236766	2.902517
3	9.59079	1.833047	0.0	0.018219	0.018219	2.482181	-0.669168	1.123778	-2.019197
4	9.59079	2.148748	0.0	0.018219	0.018219	2.188719	0.051362	1.546433	2.460967
...	...	...	...	...	...	...	...	...	...
29995	9.59079	2.049257	0.0	0.018219	0.018219	3.034526	0.578213	2.304229	-0.539174
29996	9.59079	2.358109	0.0	0.018219	0.018219	1.228739	1.090566	1.721073	-0.204690
29997	9.59079	1.742805	0.0	0.018219	0.018219	1.114133	-0.424786	0.785586	-0.548252
29998	9.59079	2.389674	0.0	0.018219	0.018219	0.630924	1.279562	0.737599	0.805103
29999	9.59079	0.451225	0.0	0.018219	0.018219	0.835934	-0.064061	2.603734	0.181360

30000 rows × 9 columns

This also means that we can use all kinds of fancy plotting functionality of for instance matplotlib.pyplot to see what’s going on. Here’s an example:

import numpy as np
from matplotlib import cm

intermediate_states = sorted(
    (
        p
        for p in model.particles
        if p not in model.kinematics.final_state.values()
        and p not in model.kinematics.initial_state.values()
    ),
    key=lambda p: p.mass,
)

evenly_spaced_interval = np.linspace(0, 1, len(intermediate_states))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]

import matplotlib.pyplot as plt

np.sqrt(data_frame["mSq_3_4"]).hist(bins=100, alpha=0.5, density=True)
plt.xlabel("$m$ [GeV]")
for i, p in enumerate(intermediate_states):
    plt.axvline(x=p.mass, linestyle="dotted", label=p.name, color=colors[i])
plt.legend();

2.5 Export data sets¶

tensorwaves currently has no export functionality for data samples. However, as we work with numpy.ndarray, it’s easy to just ‘pickle’ these data samples with numpy.save:

import numpy as np

np.save("data_sample", data_sample)
np.save("phsp_sample", phsp_sample)

In the next step, we will illustrate how to ‘perform a fit’ with tensorwaves by optimizing the intensity model to these data samples.

TensorWaves