Step 3: Perform fit¶

As explained in the previous step, a Model object and its lambdified Function form behave 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’.

First, load the relevant data from the previous steps.

import pickle

import qrules

from tensorwaves.model import SympyModel

reaction = qrules.io.load("transitions.json")
with open("helicity_model.pickle", "rb") as stream:
    model = pickle.load(stream)
with open("data_set.pickle", "rb") as stream:
    data_set = pickle.load(stream)
with open("phsp_set.pickle", "rb") as stream:
    phsp_set = pickle.load(stream)
sympy_model = SympyModel(
    expression=model.expression.doit(),
    parameters=model.parameter_defaults,
    max_complexity=100,
)

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.

Generally, the intensity model is not normalized, but a log likelihood estimator requires a normalized function. This is where the phase space dataset comes into play again: the intensity is evaluated separately with the phase space dataset so that the output of the Function can be normalized with it. The phase space sample is therefore a required argument!

from tensorwaves.estimator import UnbinnedNLL

estimator = UnbinnedNLL(
    sympy_model,
    data_set,
    phsp_set,
    backend="jax",
)

Note that the UnbinnedNLL can be expressed with different backends (it creates a LambdifiedFunction internally from the template Model). Here, we use jax, which is turns out to be the fastest backend for this model.

3.2 Optimize intensity model¶

Specify fit parameters¶

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. The optimize() method requires a mapping of parameter names to their initial values. Only the parameters listed in the mapping are optimized.

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 fitting to a data sample that was generated with this very same amplitude model after all).

initial_parameters = {
    R"C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}": 1.0
    + 0.0j,
    "Gamma_f(0)(980)": 0.15,
    "Gamma_f(0)(1500)": 0.2,
    "m_f(0)(1710)": 1.78,
}

Recall that a Function object 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 1 and 2, 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
from matplotlib import cm

reaction_info = model.adapter.reaction_info
intermediate_states = sorted(
    (
        p
        for p in model.particles
        if p not in reaction_info.final_state.values()
        and p not in reaction_info.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]


def indicate_masses():
    plt.xlabel("$m$ [GeV]")
    for i, p in enumerate(intermediate_states):
        plt.gca().axvline(
            x=p.mass, linestyle="dotted", label=p.name, color=colors[i]
        )


def compare_model(
    variable_name,
    data_set,
    phsp_set,
    intensity_model,
    bins=100,
):
    data = data_set[variable_name]
    phsp = phsp_set[variable_name]
    intensities = intensity_model(phsp_set)
    fig, ax = plt.subplots(figsize=(9, 4))
    axis = plt.gca()
    axis.hist(
        data,
        bins=bins,
        alpha=0.5,
        label="data",
        density=True,
    )
    axis.hist(
        phsp,
        weights=intensities,
        bins=bins,
        histtype="step",
        color="red",
        label="initial fit model",
        density=True,
    )
    indicate_masses()
    plt.legend()

from tensorwaves.model import LambdifiedFunction

intensity = LambdifiedFunction(sympy_model, backend="numpy")

intensity.update_parameters(initial_parameters)
compare_model("m_12", data_set, phsp_set, intensity)

Custom callbacks¶

The Minuit2 class allows one to insert certain callbacks. Callbacks are used to insert behavior into the Optimizer.optimize() method. The callbacks module provides several handy callback classes, but it’s also possible to define custom callbacks into the Optimizer by defining a custom Callback class. Here’s one that live updates a plot of the latest fit model!

import os

from IPython.display import clear_output

from tensorwaves.optimizer.callbacks import Callback

variable_to_plot = "m_12"
is_doc = "READTHEDOCS" in os.environ or "EXECUTE_NB" in os.environ


class PyplotCallback(Callback):
    def __init__(self, step_size=10):
        self.__step_size = step_size
        self.__fig = None
        self.__ax = None
        self.__latest_parameters = {}

    def on_optimize_start(self, logs):
        self.__fig, self.__ax = plt.subplots(1, figsize=(8, 5))

    def on_optimize_end(self, logs):
        self.update_plot(nbins=200)
        self.__ax = None
        self.__fig = None

    def on_iteration_end(self, iteration, logs=None):
        pass

    def on_function_call_end(self, function_call, logs):
        self.__latest_parameters = logs["parameters"]
        if is_doc:
            return
        if function_call % self.__step_size != 0:
            return
        if function_call > 75:
            return
        self.update_plot(nbins=80)
        clear_output(wait=True)
        display(plt.gcf())

    def update_plot(self, nbins: int):
        if self.__fig is None or self.__ax is None:
            self.__fig, self.__ax = plt.subplots(1, figsize=(8, 5))
        data = data_set[variable_to_plot]
        phsp = phsp_set[variable_to_plot]
        intensity.update_parameters(self.__latest_parameters)
        intensities = intensity(phsp_set)
        self.__ax.cla()
        self.__ax.hist(data, bins=nbins, alpha=0.5, label="data", density=True)
        self.__ax.hist(
            phsp,
            weights=intensities,
            bins=nbins,
            histtype="step",
            color="red",
            label="fit model",
            density=True,
        )
        self.__ax.set_xlim((0.25, 2.5))
        self.__ax.set_ylim((0, 1.9))
        indicate_masses()
        plt.gcf().legend()

In the fit example below, we want to use several callbacks together. To do so, we ‘stack’ them together with CallbackList:

from tensorwaves.optimizer.callbacks import (
    CallbackList,
    CSVSummary,
    TFSummary,
    YAMLSummary,
)

callbacks = CallbackList(
    [
        CSVSummary("fit_traceback_minuit.csv"),
        PyplotCallback(),  # comment this line for raw fit performance
        TFSummary(),
        YAMLSummary("current_fit_result.yaml"),
    ]
)

Perform fit¶

Finally, we create an Optimizer to optimize() the model (which is embedded in the Estimator). Here, we choose the Minuit2 optimizer, which is the most common optimizer in high-energy physics. Since we constructed the UnbinnedNLL with jax, can an optimize the model with analytic gradient over the LambdifiedFunction:

from tensorwaves.optimizer.minuit import Minuit2

minuit2 = Minuit2(
    callback=callbacks,
    use_analytic_gradient=False,
)
fit_result = minuit2.optimize(estimator, initial_parameters)
fit_result

FitResult(
 minimum_valid=True,
 execution_time=11.673398494720459,
 function_calls=165,
 estimator_value=-4888.608860592183,
 parameter_values={
  'Gamma_f(0)(980)': 0.06365616742326215,
  'Gamma_f(0)(1500)': 0.11174387905327125,
  'm_f(0)(1710)': 1.7095291603108285,
  'C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}': (1.0430793357711892-0.01756807653557114j),
 },
 parameter_errors={
  'Gamma_f(0)(980)': 0.003210271911638299,
  'Gamma_f(0)(1500)': 0.005785942332913687,
  'm_f(0)(1710)': 0.002462836497190077,
  'C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}': (0.03169614381083045+0.038182833991177334j),
 },
)

As can be seen, the values of the optimized parameters in the result are again comparable to the original values that are contained in the model (SympyModel.parameters):

Covariance matrix¶

Each of the optimizers offer more specific information about the fit result. This information can be accessed with FitResult.specifics. A common example would be to get the ~iminuit.Minuit.covariance matrix:

covariance_matrix = fit_result.specifics.covariance
covariance_matrix

	real_C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}	imag_C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}	Gamma_f(0)(980)	Gamma_f(0)(1500)	m_f(0)(1710)
real_C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}	0.001	-6e-05 (-0.050)	2.22e-05 (0.218)	-5.64e-05 (-0.308)	6.52e-08
imag_C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}	-6e-05 (-0.050)	0.00146	-4.64e-06 (-0.038)	-2.48e-05 (-0.112)	-3.67e-05 (-0.391)
Gamma_f(0)(980)	2.22e-05 (0.218)	-4.64e-06 (-0.038)	1.03e-05	4.04e-06 (0.217)	8.91e-08 (0.011)
Gamma_f(0)(1500)	-5.64e-05 (-0.308)	-2.48e-05 (-0.112)	4.04e-06 (0.217)	3.35e-05	3.47e-06 (0.244)
m_f(0)(1710)	6.52e-08	-3.67e-05 (-0.391)	8.91e-08 (0.011)	3.47e-06 (0.244)	6.07e-06

AIC and BIC¶

As an example, here is how to compute the AIC and BIC from this FitResult:

n_real_par = fit_result.count_number_of_parameters(complex_twice=True)
n_events = len(list(data_set.values())[0])
log_likelihood = -fit_result.estimator_value

bic = n_real_par * np.log(n_events) - 2 * log_likelihood
aic = 2 * n_real_par - 2 * log_likelihood

AIC: -9767.217721184366
BIC: -9731.166019324486

optimized_parameters = fit_result.parameter_values
for p in optimized_parameters:
    print(p)
    print(f"  initial:   {initial_parameters[p]:.3}")
    print(f"  optimized: {optimized_parameters[p]:.3}")
    print(f"  original:  {sympy_model.parameters[p]:.3}")

Gamma_f(0)(980)
  initial:   0.15
  optimized: 0.0637
  original:  0.06
Gamma_f(0)(1500)
  initial:   0.2
  optimized: 0.112
  original:  0.112
m_f(0)(1710)
  initial:   1.78
  optimized: 1.71
  original:  1.7
C_{J/\psi(1S) \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500) \to \pi^{0}_{0} \pi^{0}_{0}}
  initial:   (1+0j)
  optimized: (1.04-0.0176j)
  original:  (1+0j)

3.3 Export and import¶

In 3.2 Optimize intensity model, we initialized Minuit2 with some callbacks that are Loadable. Such callback classes offer the possibility to Loadable.load_latest_parameters(), so you can pick up the optimize process in case it crashes or if you pause it. Loading the latest parameters goes as follows:

latest_parameters = YAMLSummary.load_latest_parameters(
    "current_fit_result.yaml"
)
latest_parameters

{'Gamma_f(0)(980)': 0.06365616742326215,
 'Gamma_f(0)(1500)': 0.11174387905327125,
 'm_f(0)(1710)': 1.7095291603108285,
 'C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}': (1.0430793357711892-0.01756807653557114j)}

To restart the fit with the latest parameters, simply rerun as before.

minuit2 = Minuit2()
fit_result = minuit2.optimize(estimator, latest_parameters)
fit_result

FitResult(
 minimum_valid=True,
 execution_time=3.6539642810821533,
 function_calls=72,
 estimator_value=-4888.608862759369,
 parameter_values={
  'Gamma_f(0)(980)': 0.06365683600198997,
  'Gamma_f(0)(1500)': 0.11174546145693998,
  'm_f(0)(1710)': 1.7095299785532052,
  'C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}': (1.0431010041357933-0.01750920916912178j),
 },
 parameter_errors={
  'Gamma_f(0)(980)': 0.0032103302622486806,
  'Gamma_f(0)(1500)': 0.005786187605015158,
  'm_f(0)(1710)': 0.0024628323456297786,
  'C_{J/\\psi(1S) \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500) \\to \\pi^{0}_{0} \\pi^{0}_{0}}': (0.03169688150079737+0.03818331883381536j),
 },
)

Lo and behold: the parameters were already optimized, so the fit converged faster!

3.4 Visualize¶

Plot optimized model¶

Using the same method as above, we renew the parameters of the Function and plot it again over the data distribution.

intensity.update_parameters(latest_parameters)
compare_model("m_12", data_set, phsp_set, intensity)

Intensity components¶

Notice that the HelicityModel contains components attribute. This is dict of component names to sympy.Exprs helps us to identify amplitudes and intensities in the total amplitude.

sorted(model.components)

['A_{J/\\psi(1S)_{+1} \\to f_{0}(1370)_{0} \\gamma_{+1}; f_{0}(1370)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(1370)_{0} \\gamma_{-1}; f_{0}(1370)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(1500)_{0} \\gamma_{-1}; f_{0}(1500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(1710)_{0} \\gamma_{+1}; f_{0}(1710)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(1710)_{0} \\gamma_{-1}; f_{0}(1710)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(500)_{0} \\gamma_{+1}; f_{0}(500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(500)_{0} \\gamma_{-1}; f_{0}(500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(980)_{0} \\gamma_{+1}; f_{0}(980)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{+1} \\to f_{0}(980)_{0} \\gamma_{-1}; f_{0}(980)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1370)_{0} \\gamma_{+1}; f_{0}(1370)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1370)_{0} \\gamma_{-1}; f_{0}(1370)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1500)_{0} \\gamma_{+1}; f_{0}(1500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1500)_{0} \\gamma_{-1}; f_{0}(1500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1710)_{0} \\gamma_{+1}; f_{0}(1710)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(1710)_{0} \\gamma_{-1}; f_{0}(1710)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(500)_{0} \\gamma_{+1}; f_{0}(500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(500)_{0} \\gamma_{-1}; f_{0}(500)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(980)_{0} \\gamma_{+1}; f_{0}(980)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'A_{J/\\psi(1S)_{-1} \\to f_{0}(980)_{0} \\gamma_{-1}; f_{0}(980)_{0} \\to \\pi^{0}_{0} \\pi^{0}_{0}}',
 'I_{J/\\psi(1S)_{+1} \\to \\gamma_{+1} \\pi^{0}_{0} \\pi^{0}_{0}}',
 'I_{J/\\psi(1S)_{+1} \\to \\gamma_{-1} \\pi^{0}_{0} \\pi^{0}_{0}}',
 'I_{J/\\psi(1S)_{-1} \\to \\gamma_{+1} \\pi^{0}_{0} \\pi^{0}_{0}}',
 'I_{J/\\psi(1S)_{-1} \\to \\gamma_{-1} \\pi^{0}_{0} \\pi^{0}_{0}}']

model.components[
    R"A_{J/\psi(1S)_{+1} \to f_{0}(1370)_{0} \gamma_{+1}; f_{0}(1370)_{0} \to"
    R" \pi^{0}_{0} \pi^{0}_{0}}"
].subs(model.parameter_defaults).doit()

\[\displaystyle \frac{0.23625 \sqrt{2} \sqrt{\frac{\left(m_{012}^{2} - \left(- m_{0} + m_{12}\right)^{2}\right) \left(m_{012}^{2} - \left(m_{0} + m_{12}\right)^{2}\right)}{m_{012}^{2} \left(1 + \frac{\left(m_{012}^{2} - \left(- m_{0} + m_{12}\right)^{2}\right) \left(m_{012}^{2} - \left(m_{0} + m_{12}\right)^{2}\right)}{4 m_{012}^{2}}\right)}} \left(\frac{1}{2} - \frac{\cos{\left(\theta_{1+2} \right)}}{2}\right) e^{i \phi_{1+2}}}{- m_{12}^{2} + 1.8225 - \frac{0.86113125 i \sqrt{\frac{\left(m_{12}^{2} - \left(m_{1} - m_{2}\right)^{2}\right) \left(m_{12}^{2} - \left(m_{1} + m_{2}\right)^{2}\right)}{m_{12}^{2}}}}{\sqrt{\left(1.8225 - \left(m_{1} - m_{2}\right)^{2}\right) \left(1.8225 - \left(m_{1} + m_{2}\right)^{2}\right)} \left|{m_{12}}\right|}}\]

Just like in 2.2 Generate intensity-based sample, these intensity components can each be expressed in a computational backend. This can be done with the method sum_components() from the original model and creating a new LambdifiedFunction from that expression. After that we, update the parameters with the optimized parameter values we found in Perform fit. The two components in the example below should be the same:

added_components = model.sum_components(
    components=[
        R"A_{J/\psi(1S)_{+1} \to f_{0}(500)_{0} \gamma_{+1}; f_{0}(500)_{0}"
        R" \to \pi^{0}_{0} \pi^{0}_{0}}",
        R"A_{J/\psi(1S)_{+1} \to f_{0}(980)_{0} \gamma_{+1}; f_{0}(980)_{0}"
        R" \to \pi^{0}_{0} \pi^{0}_{0}}",
        R"A_{J/\psi(1S)_{+1} \to f_{0}(1370)_{0} \gamma_{+1}; f_{0}(1370)_{0}"
        R" \to \pi^{0}_{0} \pi^{0}_{0}}",
        R"A_{J/\psi(1S)_{+1} \to f_{0}(1500)_{0} \gamma_{+1}; f_{0}(1500)_{0}"
        R" \to \pi^{0}_{0} \pi^{0}_{0}}",
        R"A_{J/\psi(1S)_{+1} \to f_{0}(1710)_{0} \gamma_{+1}; f_{0}(1710)_{0}"
        R" \to \pi^{0}_{0} \pi^{0}_{0}}",
    ]
)
from_amplitudes = LambdifiedFunction(
    model=SympyModel(
        expression=added_components.doit(),
        parameters=model.parameter_defaults,
    ),
    backend="numpy",
)
from_amplitudes.update_parameters(latest_parameters)

added_components = model.sum_components(
    components=[R"I_{J/\psi(1S)_{+1} \to \gamma_{+1} \pi^{0}_{0} \pi^{0}_{0}}"]
)
from_intensity = LambdifiedFunction(
    model=SympyModel(
        expression=added_components.doit(),
        parameters=model.parameter_defaults,
    ),
    backend="numpy",
)
from_intensity.update_parameters(latest_parameters)

import numpy as np

difference = np.average(from_amplitudes(phsp_set) - from_intensity(phsp_set))
assert np.round(difference, decimals=15) == 0.0

The result is a LambdifiedFunction that can be plotted just like in Plot optimized model:

fig, ax = plt.subplots(1, figsize=(8, 5))
bins = 150
ax.hist(
    phsp_set["m_12"],
    weights=intensity(phsp_set),
    bins=bins,
    alpha=0.2,
    label="full intensity",
)
ax.hist(
    phsp_set["m_12"],
    weights=from_intensity(phsp_set),
    bins=bins,
    histtype="step",
    label=R"$J/\psi(1S)_{-1} \to \gamma_{-1} \pi^0 \pi^0$",
)
ax.set_xlim(0.25, 2.5)
ax.set_xlabel(R"$m_{\pi^0\pi^0}$ [GeV]")
ax.set_yticks([])
plt.legend();

Or generically, so that we can stack all the sub-intensities:

import logging

from tensorwaves.data import generate_data

logging.basicConfig()
logging.getLogger().setLevel(logging.ERROR)
intensity_components = [
    LambdifiedFunction(
        SympyModel(
            expression=model.sum_components([c]).doit(),
            parameters=model.parameter_defaults,
        ),
        backend="numpy",
    )
    for c in model.components
    if c.startswith("I")
]
reaction_info = model.adapter.reaction_info
initial_state_mass = reaction_info.initial_state[-1].mass
final_state_masses = {i: p.mass for i, p in reaction_info.final_state.items()}
sub_intensities = [
    generate_data(
        size=5_000,
        initial_state_mass=initial_state_mass,
        final_state_masses=final_state_masses,
        data_transformer=model.adapter,
        intensity=component,
    )
    for component in intensity_components
]

fig, ax = plt.subplots(1, figsize=(8, 5))
plt.hist(
    [model.adapter.transform(i)["m_12"] for i in sub_intensities],
    bins=100,
    stacked=True,
    alpha=0.6,
)
ax.set_xlim(0.25, 2.5)
ax.set_xlabel(R"$m_{\pi^0\pi^0}$ [GeV]")
ax.set_yticks([])
ax.legend(
    labels=[f"${c[3:-1]}$" for c in model.components if c.startswith("I")]
);

Analyze optimization process¶

Note that in Step 3.2, we initialized Minuit2 with a TFSummary callback as well. Its output files provide a nice, interactive representation of the fit process and can be viewed with TensorBoard as follows:

Terminal

tensorboard --logdir logs

Python

import tensorboard as tb

tb.notebook.list()  # View open TensorBoard instances
tb.notebook.start(args_string="--logdir logs")

See more info here

Jupyter notebook

%load_ext tensorboard
%tensorboard --logdir logs