Step 3: Perform fit¶
Warning
The pycompwa is no longer maintained. Use the ComPWA packages QRules, AmpForm, and TensorWaves instead!
As explained in the previous step, an Intensity object behaves just like a mathematical function that takes a DataSet as an argument and returns a list of intensities (real numbers). We now 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 pycompwa.ui as pwa
2021-10-27 10:27:35,894 INFO [default] Logging to file disabled!
2021-10-27 10:27:35,894 [INFO] Log level: INFO
2021-10-27 10:27:35,894 [INFO] Current date and time: Wed Oct 27 10:27:35 2021
particle_list = pwa.read_particles("model.xml")
kinematics = pwa.create_helicity_kinematics("model.xml", particle_list)
kinematics.create_all_subsystems()
data_sample = pwa.read_root_data(input_file="generated_data.root")
phsp_sample = pwa.read_root_data(input_file="generated_phsp.root")
intensity_builder = pwa.IntensityBuilderXML(
"model.xml", particle_list, kinematics, phsp_sample
)
intensity = intensity_builder.create_intensity()
2021-10-27 10:27:36,292 [INFO] HelicityKinematics::HelicityKinematics() | Initialized kinematics for reaction ( J/psi )->( gamma[ID=2] pi0[ID=3] pi0[ID=4] )
Event position to final state ID mapping:
0: 2
1: 3
2: 4
2021-10-27 10:27:36,296 [INFO] creating all Subsystems!
2021-10-27 10:27:36,728 [INFO] Setting phase space sample weights...
2021-10-27 10:27:36,738 [INFO] Updating data container content...
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. To create an estimator we therefore generally need:
The data sample of momentum tuples can be converted to a DataSet as follows:
data_set = kinematics.convert(data_sample)
Now we can create an estimator. In this example, we use an unbinned log likelihood estimator. This estimator is used most commonly in high-energy physics.
(
estimator,
initial_parameters,
) = pwa.create_unbinned_log_likelihood_function_tree_estimator(
intensity, data_set
)
Notice that you not only receive a estimator object, but also a list of fit parameters. You use this list of fit parameters to initialize the optimizer later on.
The parameters contains the following info:
initial value
whether they are free or fixed
optional upper and lower limits
errors that give certain optimizers hints on the step size
initial_parameters
[strength_incoherent: 1 (fixed),
mass_f21270: 1.2751 +- 0.0012 (fixed),
Mass_gamma: 0 (fixed),
Mass_jpsi: 3.0969 (fixed),
Mass_neutralPion: 0.134977 +- 6e-06 (fixed),
width_f21270: 0.185 (fixed),
Radius_f21270: 1 (fixed),
Magnitude_J/psi_to_f2(1270)_0+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1270)_0+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 0,
mass_f21950: 1.944 +- 0.012 (fixed),
width_f21950: 0.472 (fixed),
Radius_f21950: 1 (fixed),
Magnitude_J/psi_to_f2(1950)_0+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1950)_0+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 0,
Magnitude_J/psi_to_f2(1270)_-1+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1270)_-1+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 0,
Magnitude_J/psi_to_f2(1950)_-1+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1950)_-1+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 0,
Magnitude_J/psi_to_f2(1270)_-2+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1270)_-2+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 0,
Magnitude_J/psi_to_f2(1950)_-2+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f2(1950)_-2+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 0,
mass_f0980: 0.994 +- 0.001 (fixed),
width_f0980: 0.07 (fixed),
Radius_f0980: 1 (fixed),
Magnitude_J/psi_to_f0(980)_0+gamma_-1;f0(980)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f0(980)_0+gamma_-1;f0(980)_to_pi0_0+pi0_0;: 0,
mass_f01500: 1.505 +- 0.006 (fixed),
width_f01500: 0.109 (fixed),
Radius_f01500: 1 (fixed),
Magnitude_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;: 1,
Phase_J/psi_to_f0(1500)_0+gamma_-1;f0(1500)_to_pi0_0+pi0_0;: 0]
3.2 Modify parameters¶
The fit parameters in this list have been initialized with the values stated in the XML file or with default values if unspecified, but they can easily be modified in Python as well:
print("\nthis parameter is initially not fixed:")
print(initial_parameters[8])
initial_parameters[7].is_fixed = True
initial_parameters[8].is_fixed = True
print("\nand now it is fixed:")
print(initial_parameters[8])
this parameter is initially not fixed:
Phase_J/psi_to_f2(1270)_0+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 0
and now it is fixed:
Phase_J/psi_to_f2(1270)_0+gamma_-1;f2(1270)_to_pi0_0+pi0_0;: 0 (fixed)
To make the fit process a bit more interesting, we modify one of the parameters to a different initial value than the true value:
original_value = initial_parameters[12].value
initial_parameters[12]
Magnitude_J/psi_to_f2(1950)_0+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 1
initial_parameters[12].value = 2.0
initial_parameters[12]
Magnitude_J/psi_to_f2(1950)_0+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 2
3.3 Optimize the intensity model¶
Now it’s time to perform the fit! This is quite simple: just create an optimizer instance of your choice, here Minuit2, and call its optimize() method to start the fitting process.
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.
print("data size", len(data_sample.events))
print("phsp size", len(phsp_sample.events))
print(
"free parameters:",
len([par for par in initial_parameters if not par.is_fixed]),
)
data size 10000
phsp size 100000
free parameters: 14
Note
With 14 free parameters, this fit takes about 15 seconds on an Intel(R) Core(TM) i7-6820HQ 2.70GHz CPU.
optimizer = pwa.MinuitIF()
fit_result = optimizer.optimize(estimator, initial_parameters)
2021-10-27 10:27:38,987 [INFO] MinuitIF::optimize() | Number of parameters (free): 32 (14)
2021-10-27 10:27:38,988 [INFO] Minuit2 strategy: medium
GradientNCycles: 3
GradientStepTolerance: 0.3
GradientTolerance: 0.05
HessianNCycles: 5
HessianGradientNCycles: 2
HessianStepTolerance: 0.3
HessianG2Tolerance: 0.05
2021-10-27 10:27:38,989 [INFO] MinuitIF::optimize() | Starting migrad: maxCalls=0 tolerance=0.1
2021-10-27 10:27:54,736 [INFO] MinuitIF::optimize() | Migrad finished! Minimum is valid = 1
2021-10-27 10:27:54,739 [INFO] MinuitIF::optimize() | Starting hesse
2021-10-27 10:27:57,930 [INFO] MinuitIF::optimize() | Hesse finished
2021-10-27 10:27:57,933 [INFO] MinuitIF::optimize() | Minimization finished! LH = -5496.086088
fit_result.fit_duration_in_seconds
18
Let’s check whether the fit parameters are close to the true values.
value = fit_result.final_parameters[12].value
error = fit_result.final_parameters[12].error
if abs(value - original_value) < 3 * error[0]:
print("This value is close to", original_value, "again:")
print(value, "+-", error[0]) # upper and lower error are equal in hesse
else:
raise Exception("Value does not lie within 3 sigma!")
This value is close to 1.0 again:
0.9747981228967589 +- 0.10721062179072079
Warning
The intensity instance still needs to be notified about this optimal set of parameters. They can be applied with updateParametersFrom().
intensity.updateParametersFrom(fit_result.final_parameters)
3.4 Exporting and importing¶
Optimizing an intensity model can take a long time, so it is important that you store the fit result in the end. We again write to an XML file (but note that this has a different format than that of the amplitude model).
fit_result.write("fit_result.xml")
imported_fit_result = pwa.load("fit_result.xml")
exported_value = fit_result.final_parameters[12].value
imported_value = imported_fit_result.final_parameters[12].value
if exported_value != imported_value:
raise Exception
else:
print(imported_fit_result.final_parameters[12])
Magnitude_J/psi_to_f2(1950)_0+gamma_-1;f2(1950)_to_pi0_0+pi0_0;: 0.974798 +- 0.107211