Step 3: Perform fit¶
As explained in the previous step, an IntensityTF object behaves just like a mathematical function that takes a set of data points as an argument and returns a list of intensities (real numbers). At this stage, we want to optimize the parameters of the intensity model, so that it matches the distribution of our data sample. This is what we call ‘performing a fit’.
Note
We first load the relevant data from the previous steps.
import numpy as np
from expertsystem import io
from tensorwaves.physics.helicity_formalism.amplitude import IntensityBuilder
from tensorwaves.physics.helicity_formalism.kinematics import HelicityKinematics
model = io.load_amplitude_model("amplitude_model_helicity.yml")
kinematics = HelicityKinematics.from_model(model)
data_sample = np.load("data_sample.npy")
phsp_sample = np.load("phsp_sample.npy")
intensity_builder = IntensityBuilder(model.particles, kinematics, phsp_sample)
intensity = intensity_builder.create_intensity(model)
phsp_set = kinematics.convert(phsp_sample)
data_set = kinematics.convert(data_sample)
3.1 Define estimator¶
To perform a fit, you need to define an estimator. An estimator is a measure for the discrepancy between the intensity model and the data distribution to which you fit. In PWA, we usually use an unbinned negative log likelihood estimator, which can be created as follows with the tensorwaves.estimator module:
from tensorwaves.estimator import UnbinnedNLL
estimator = UnbinnedNLL(intensity, data_set)
The parameters that the UnbinnedNLL carries, have been taken from the IntensityTF object and you’ll recognize them from the YAML recipe file that we created in Step 1:
estimator.parameters
{'MesonRadius_J/psi(1S)': 1.0,
'Magnitude_J/psi(1S)_to_f(0)(1710)_0+gamma_1;f(0)(1710)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi(1S)_to_f(0)(1710)_0+gamma_1;f(0)(1710)_to_pi0_0+pi0_0;': 0.0,
'Magnitude_J/psi(1S)_to_f(0)(980)_0+gamma_1;f(0)(980)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi(1S)_to_f(0)(980)_0+gamma_1;f(0)(980)_to_pi0_0+pi0_0;': 0.0,
'Magnitude_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;': 0.0,
'Magnitude_J/psi(1S)_to_f(0)(1370)_0+gamma_1;f(0)(1370)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi(1S)_to_f(0)(1370)_0+gamma_1;f(0)(1370)_to_pi0_0+pi0_0;': 0.0,
'Magnitude_J/psi(1S)_to_f(0)(500)_0+gamma_1;f(0)(500)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi(1S)_to_f(0)(500)_0+gamma_1;f(0)(500)_to_pi0_0+pi0_0;': 0.0,
'strength_incoherent': 1.0,
'Mass_J/psi(1S)': 3.0969,
'Width_J/psi(1S)': 9.29e-05,
'Mass_f(0)(1710)': 1.704,
'Width_f(0)(1710)': 0.123,
'MesonRadius_f(0)(1710)': 1.0,
'Mass_f(0)(980)': 0.99,
'Width_f(0)(980)': 0.06,
'MesonRadius_f(0)(980)': 1.0,
'Mass_f(0)(1500)': 1.506,
'Width_f(0)(1500)': 0.112,
'MesonRadius_f(0)(1500)': 1.0,
'Mass_f(0)(1370)': 1.35,
'Width_f(0)(1370)': 0.35,
'MesonRadius_f(0)(1370)': 1.0,
'Mass_f(0)(500)': 0.475,
'Width_f(0)(500)': 0.55,
'MesonRadius_f(0)(500)': 1.0}
3.3 Optimize the intensity model¶
Now it’s time to perform the fit! This is the part where the TensorFlow backend comes into play. For this reason, we suppress a few warnings:
import tensorflow as tf
tf.get_logger().setLevel("ERROR")
Starting the fit itself is quite simple: just create an optimizer instance of your choice, here Minuit2, and call its optimize method to start the fitting process. Notice that the optimize method requires a second argument. This is a mapping of parameter names that you want to fit to their initial values.
Let’s first select a few of the parameters that we saw in Step 3.1 and feed them to the optimizer to run the fit. Notice that we modify the parameters slightly to make the fit more interesting (we are running using a data sample that was generated with this very amplitude model after all).
initial_parameters = {
"Phase_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;": 0.0,
"Width_f(0)(500)": 0.1,
"Mass_f(0)(980)": 0.9,
"Mass_f(0)(1500)": 1.55,
"Mass_f(0)(1710)": 1.8,
"Width_f(0)(1710)": 0.3,
}
Recall that an IntensityTF object is really just a Function that computes the intensity for a certain dataset. This can be seen now nicely when we use these intensities as weights on the phase space sample and plot it together with the original dataset. Here, we look at the invariant mass distribution projection of the final states 3 and 4, which, as we saw before, are the final state particles \(\pi^0,\pi^0\).
Don’t forget to use update_parameters first!
import matplotlib.pyplot as plt
import numpy as np
def compare_model(
variable_name, data_set, phsp_set, intensity_model, transform=lambda p: p, bins=100
):
data = transform(data_set[variable_name])
phsp = transform(phsp_set[variable_name])
intensities = intensity_model(phsp_set)
plt.hist(data, bins=bins, alpha=0.5, label="data", density=True)
plt.hist(
phsp,
weights=intensities,
bins=bins,
histtype="step",
color="red",
label="initial fit model",
density=True,
)
plt.legend()
intensity.update_parameters(initial_parameters)
compare_model("mSq_3_4", data_set, phsp_set, intensity, np.sqrt)
Finally, we create an Optimizer to optimize the model (which is embedded in the Estimator. Only the parameters that are given to the optimize method are optimized.
from tensorwaves.optimizer.minuit import Minuit2
minuit2 = Minuit2()
result = minuit2.optimize(estimator, initial_parameters)
result
{'params': {'Phase_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;': (-0.023798691615680538,
0.03836138218557237),
'Width_f(0)(500)': (0.5597695025691618, 0.015149731647513634),
'Mass_f(0)(980)': (0.9882650309317992, 0.000916535928361252),
'Mass_f(0)(1500)': (1.5048838543453684, 0.0018579090480874351),
'Mass_f(0)(1710)': (1.7031612363982251, 0.0012262007503469697),
'Width_f(0)(1710)': (0.12067176356289812, 0.003731431127665419)},
'log_lh': -13356.093975609852,
'iterations': 242,
'func_calls': 242,
'time': 62.56543302536011}
Computation time ― around a minute in this case
The computation time depends on the complexity of the model, number of data events, the size of the phase space sample, and the number of free parameters. This model is rather small and has but a few free parameters, so the optimization shouldn’t take more than a few minutes.
Using the same method as above, we renew the parameters of the IntensityTF and plot it again over the data distribution.
intensity.update_parameters(initial_parameters)
compare_model("mSq_3_4", data_set, phsp_set, intensity, np.sqrt)
3.4 Export fit result¶
Not implemented yet
This functionality has not yet been implemented. See issue 99.
Optimizing an intensity model can take a long time, so it is important that you store the fit result in the end. We can simply dump the file to json:
import json
with open("tensorwaves_fit_result.json", "w") as stream:
json.dump(result, stream, indent=2)
…and load it back:
with open("tensorwaves_fit_result.json") as stream:
result = json.load(stream)
result
{'params': {'Phase_J/psi(1S)_to_f(0)(1500)_0+gamma_1;f(0)(1500)_to_pi0_0+pi0_0;': [-0.023798691615680538,
0.03836138218557237],
'Width_f(0)(500)': [0.5597695025691618, 0.015149731647513634],
'Mass_f(0)(980)': [0.9882650309317992, 0.000916535928361252],
'Mass_f(0)(1500)': [1.5048838543453684, 0.0018579090480874351],
'Mass_f(0)(1710)': [1.7031612363982251, 0.0012262007503469697],
'Width_f(0)(1710)': [0.12067176356289812, 0.003731431127665419]},
'log_lh': -13356.093975609852,
'iterations': 242,
'func_calls': 242,
'time': 62.56543302536011}