Download this notebook here.
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.
[1]:
import numpy as np
import yaml
from tensorwaves.physics.helicity_formalism.kinematics import HelicityKinematics
from tensorwaves.physics.particle import extract_particles
from tensorwaves.physics.helicity_formalism.amplitude import IntensityBuilder
[2]:
with open("amplitude_model_helicity.yml") as input_file:
recipe = yaml.load(input_file.read(), Loader=yaml.SafeLoader)
kinematics = HelicityKinematics.from_recipe(recipe)
particles = extract_particles(recipe)
data_sample = np.load("data_sample.npy")
phsp_sample = np.load("phsp_sample.npy")
intensity_builder = IntensityBuilder(particles, kinematics, phsp_sample)
intensity = intensity_builder.create_intensity(recipe)
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:
[3]:
from tensorwaves.estimator import UnbinnedNLL
estimator = UnbinnedNLL(intensity, data_set)
The parameters that the estimator.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:
[4]:
estimator.parameters
[4]:
{'strength_incoherent': 1.0,
'Magnitude_J/psi_to_f0(980)_0+gamma_-1;f0(980)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi_to_f0(980)_0+gamma_-1;f0(980)_to_pi0_0+pi0_0;': 0.0,
'Magnitude_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;': 1.0,
'Phase_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;': 0.0,
'Mass_J/psi': 3.0969,
'Width_J/psi': 0.0,
'MesonRadius_J/psi': 1.0,
'Mass_f0(980)': 0.994,
'Width_f0(980)': 0.07,
'MesonRadius_f0(980)': 1.0,
'Mass_f0(1500)': 1.505,
'Width_f0(1500)': 0.109,
'MesonRadius_f0(1500)': 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:
[5]:
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. Parameters not listed in that `dict`, are not optimized.
So let’s select a few of the parameters that we saw in 3.1 Define estimator 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).
[6]:
initial_parameters = {
"Phase_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;": 0.0,
"Mass_f0(980)": 1.1,
"Width_f0(980)": 0.1,
"Mass_f0(1500)": 1.4,
"Width_f0(1500)": 0.2,
}
[7]:
from tensorwaves.optimizer.minuit import Minuit2
minuit2 = Minuit2()
result = minuit2.optimize(estimator, initial_parameters)
result
[7]:
{'params': {'Phase_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;': (0.01732103038947869,
0.029060391801374985),
'Mass_f0(980)': (0.9941593955051807, 0.0004910702276528474),
'Width_f0(980)': (0.06961706764720732, 0.0007732360735722078),
'Mass_f0(1500)': (1.5045574190841913, 0.0005611945200997931),
'Width_f0(1500)': (0.11017162926673399, 0.001061707129755457)},
'log_lh': -34640.28968481776,
'iterations': 214,
'func_calls': 214,
'time': 66.16642260551453}
Note
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.