Introduction to zfit

zfit is a model building and fitting library, mostly for likelihood fits and with advanced features for HEP analyses.

It aims to define an API and workflow to be the stable core and allow other libraries to be built on top of it. Furthermore, it uses JIT compilation and autograd (current backend TensorFlow, future JAX (?)) to significantly speed up the computation.

Workflow

There are five mostly independent parts in zfit

import matplotlib.pyplot as plt
import mplhep
import numpy as np
import zfit
# znp is a subset of numpy functions with a numpy interface but using actually the zfit backend (currently TF)
import zfit.z.numpy as znp
import zfit_physics as zphysics  # extension (for user contribution) containing physics PDFs
from numpy import floating

Data

This component in general plays a minor role in zfit: it is mostly to provide a unified interface for data.

Preprocessing is therefore not part of zfit and should be done beforehand. Python offers many great possibilities to do so (e.g. Pandas).

zfit knows unbinned and binned datasets.

The unbinned Data can load data from various sources, such as Numpy, Pandas DataFrame, TensorFlow Tensor and ROOT (using uproot). It is also possible, for convenience, to convert it directly to_pandas. The constructors are named from_numpy, from_root etc.

A Data needs not only the data itself but also the observables: the human-readable string identifiers of the axes (corresponding to “columns” of a Pandas DataFrame). It is convenient to define the Space not only with the observable but also with a limit: this can directly be re-used as the normalization range in the PDF.

First, let’s define our observables

obs = zfit.Space('obs1', -5, 10)
# or more explicitly
obs = zfit.Space(obs='obs1', lower=-5, upper=10)

This Space has limits and offers the following functionality:

• volume: calculate the volume (e.g. distance between upper and lower)

• inside(): return a boolean Tensor corresponding to whether the value is inside the Space

• filter(): filter the input values to only return the one inside

size_normal = 10_000
data_normal_np = np.random.normal(size=size_normal, scale=2)

data_normal = zfit.Data(data_normal_np, obs=obs)

The main functionality is

• nevents: attribute that returns the number of events in the object

• space: a Space that defines the limits of the data; if outside, the data will be cut (!)

• n_obs: the number of dimensions in the dataset

• with_obs: returns a dataset, possibly a subset of the dataset with only the given obs

• with_weights: returns a dataset with the given weights

• weights: event weights

Using indexing data[obs] will return a Data object with only the given observables and follows a similar behavior as Pandas DataFrames. A string, i.e. a single observable, returns an array, using alist of strings returns a Data object with the given observables.

Furthermore, value returns a Tensor with shape (nevents, n_obs).

print(f"We have {data_normal.nevents} events in our dataset with the minimum of {np.min(data_normal.value())}")  # remember! The obs cut out some of the data

data_normal.n_obs

data_normal["obs1"]  # indexing Pandas DataFrame style

Model

Building models is by far the largest part of zfit. We will therefore cover an essential part, the possibility to build custom models, in an extra chapter. Let’s start out with the idea that you define your parameters and your observable space; the latter is the expected input data.

PDFs are normalized (over a specified range); this is the main model and is what we are going to use throughout the tutorials.

A PDF is defined by

(1)\[\begin{align} \mathrm{PDF}_{f(x)}(x; \theta) = \frac{f(x; \theta)}{\int_{a}^{b} f(x; \theta)} \end{align}\]

where a and b define the normalization range (norm), over which (by inserting into the above definition) the integral of the PDF is unity.

Let’s create a simple Gaussian PDF. We already defined the Space for the data before, now we only need the parameters. These are a different objects than a Space.

Parameter

A Parameter (there are different kinds actually, more on that later) takes the following arguments as input: Parameter(name, initial value[, lower limit, upper limit]) where the limits are highly recommended but not mandatory. Furthermore, step_size can be given (which is useful to be around the given uncertainty, e.g. for large yields or small values it can help a lot to set this). Also, a floating argument is supported, indicating whether the parameter is allowed to float in the fit or not (just omitting the limits does not make a parameter constant).

Parameters have a unique name. This is served as the identifier for e.g. fit results.

mu = zfit.Parameter('mu', 1, -3, 3, step_size=0.2)
sigma_num = zfit.Parameter('sigma42', 1, 0.1, 10, floating=False)

Now we have everything to create a Gaussian PDF:

gauss = zfit.pdf.Gauss(obs=obs, mu=mu, sigma=sigma_num)

Since this holds all the parameters and the observables are well-defined, we can retrieve them

gauss.n_obs  # dimensions

gauss.obs

gauss.space

gauss.norm

As we’ve seen, the obs we defined is the space of Gauss: this acts as the default limits whenever needed (e.g. for sampling). gauss also has a norm, which equals by default as well to the obs given, however, we can explicitly change that.

We can also access the parameters of the PDF in two ways, depending on our intention: either by name (the parameterization name, e.g. mu and sigma, as defined in the Gauss), which is useful if we are interested in the parameter that describes the shape

gauss.params

or to retrieve all the parameters that the PDF depends on. As this now may sounds trivial, we will see later that models can depend on other models (e.g. sums) and parameters on other parameters. There is one function that automatically retrieves all dependencies, get_params. It takes three arguments to filter:

• floating: whether to filter only floating parameters, only non-floating or don’t discriminate

• is_yield: if it is a yield, or not a yield, or both

• extract_independent: whether to recursively collect all parameters. This, and the explanation for why independent, can be found later on in the Simultaneous tutorial.

Usually, the default is exactly what we want if we look for all free parameters that this PDF depends on.

gauss.get_params()

The difference will also be clear if we e.g. use the same parameter twice:

gauss_only_mu = zfit.pdf.Gauss(obs=obs, mu=mu, sigma=mu)
print(f"params={gauss_only_mu.params}")
print(f"get_params={gauss_only_mu.get_params()}")

Functionality

PDFs provide a few useful methods. The main features of a zfit PDF are:

• pdf: the normalized value of the PDF. It takes an argument norm that can be set to False, in which case we retrieve the unnormalized value

• integrate: given a certain range, the PDF is integrated. As pdf, it takes a norm argument that integrates over the unnormalized pdf if set to False

• sample: samples from the pdf and returns a Data object

integral = gauss.integrate(limits=(-1, 3))  # corresponds to 2 sigma integral
integral

Tensors

As we see, many zfit functions return Tensors. Usually, these are array-Like objects that resemble numpy arrays close enough that they can be directly used in-place of actual numpy arrays.

If explicitly needed, we can always safely convert them to a numpy array by calling np.asarray(tensor).

More information about the backend can be found in Extended tutorial on TensorFlow

np.sqrt(integral)  # works out-of-the-box

They also have shapes, dtypes, can be slices etc. So do not convert them except you need it. More on this can be seen in the talk later on about zfit and TensorFlow 2.0.

sample = gauss.sample(n=1000)  # default space taken as limits
sample

sample['obs1'][:10]
# sample.value()[:10, 0]  # "value" returns a numpy array-like object

sample.n_obs

sample.obs

We see that sample returns also a zfit Data object with the same space as it was sampled in. This can directly be used e.g.

probs = gauss.pdf(sample)
probs[:10]

NOTE: In case you want to do this repeatedly (e.g. for toy studies), there is a more efficient way (see later on)

Plotting

So far, we have a dataset and a PDF. Before we go for fitting, we can make a plot. This functionality is not directly provided in zfit (but can be added to zfit-physics). It is however simple enough to do it:

def plot_model(model, data, scale=1, plot_data=True):  # we will use scale later on

    nbins = 50

    lower, upper = data.data_range.limit1d
    x = znp.linspace(lower, upper, num=1000)  # np.linspace also works
    y = model.pdf(x) * size_normal / nbins * data.data_range.volume
    y *= scale
    plt.plot(x, y, label='model')
    if plot_data:
        # legacy way, also works
        # data_plot = data.value()[:, 0]  # we could also use the `to_pandas` method
        # plt.hist(data_plot, bins=nbins)
        # modern way
        data_binned = data.to_binned(nbins)
        mplhep.histplot(data_binned, label='data')
    plt.legend()

plot_model(gauss, data_normal)

We can of course do better (and will see that later on, continuously improve the plots), but this is quite simple and gives us the full power of matplotlib/mplhep.

Different models

zfit offers a selection of predefined models (and extends with models from zfit-physics that contain physics specific models such as ARGUS shaped models).

print(zfit.pdf.__all__)

print(zphysics.pdf.__all__)

To create a more realistic model, we can build some components for a mass fit with a

• signal component: CrystalBall

• combinatorial background: Exponential

• partial reconstructed background on the left: Kernel Density Estimation

Extended tutorial on different KDEs in zfit

Shape fit

mass_obs = zfit.Space('mass', 0, 1000)

# Signal component

mu_sig = zfit.Parameter('mu_sig', 400, 100, 600)
sigma_sig = zfit.Parameter('sigma_sig', 50, 1, 100)
alpha_sig = zfit.Parameter('alpha_sig', 200, 100, 400, floating=False)  # won't be used in the fit
npar_sig = zfit.Parameter('n sig', 4, 0.1, 30, floating=False)
signal = zfit.pdf.CrystalBall(obs=mass_obs, mu=mu_sig, sigma=sigma_sig, alpha=alpha_sig, n=npar_sig)

# combinatorial background

lam = zfit.Parameter('lambda', -0.01, -0.05, -0.0001)
comb_bkg = zfit.pdf.Exponential(lam, obs=mass_obs)

part_reco_data = np.random.normal(loc=150, scale=160, size=70_000)
part_reco = zfit.pdf.KDE1DimGrid(obs=mass_obs, data=part_reco_data)

Composing models

We can also compose multiple models together. Here we’ll stick to one dimensional models, the extension to multiple dimensions is explained in the “custom models tutorial”.

Here we will use a SumPDF. This takes pdfs and fractions. If we provide n pdfs and:

• n - 1 fracs: the nth fraction will be 1 - sum(fracs)

• n fracs: no normalization attempt is done by SumPDf. If the fracs are not implicitly normalized, this can lead to bad fitting behavior if there is a degree of freedom too much

Extended tutorial on composed models, including multidimensional products

sig_frac = zfit.Parameter('sig_frac', 0.3, 0, 1)
comb_bkg_frac = zfit.Parameter('comb_bkg_frac', 0.25, 0, 1)
model = zfit.pdf.SumPDF([signal, comb_bkg, part_reco], [sig_frac, comb_bkg_frac])

In order to have a corresponding data sample, we can just create one. Since we want to fit to this dataset later on, we will create it with slightly different values. Therefore, we can use the ability of a parameter to be set temporarily to a certain value with

print(f"before: {sig_frac}")
with sig_frac.set_value(0.25):
    print(f"new value: {sig_frac}")
print(f"after 'with': {sig_frac}")

While this is useful, it does not fully scale up. We can use the zfit.param.set_values helper therefore. (Sidenote: instead of a list of values, we can also use a FitResult, the given parameters then take the value from the result)

with zfit.param.set_values([mu_sig, sigma_sig, sig_frac, comb_bkg_frac, lam], [370, 34, 0.08, 0.31, -0.0025]):
    data = model.sample(n=10_000)

plot_model(model, data);

Plotting the components is not difficult now: we can either just plot the pdfs separately (as we still can access them) or in a generalized manner by accessing the pdfs attribute:

def plot_comp_model(model, data):
    for mod, frac in zip(model.pdfs, model.params.values()):
        plot_model(mod, data, scale=frac, plot_data=False)
    plot_model(model, data)

plot_comp_model(model, data)

Now we can add legends etc. Btw, did you notice that actually, the frac params are zfit Parameters? But we just used them as if they were Python scalars and it works.

print(model.params)

Extended PDFs

So far, we have only looked at normalized PDFs that do contain information about the shape but not about the absolute scale. We can make a PDF extended by adding a yield to it.

The behavior of the new, extended PDF does NOT change, any methods we called before will act the same. Only exception, some may require an argument less now. All the methods we used so far will return the same values. What changes is that the flag model.is_extended now returns True. Furthermore, we have now a few more methods that we can use which would have raised an error before:

• get_yield: return the yield parameter (notice that the yield is not added to the shape parameters params)

• ext_{pdf,integrate}: these methods return the same as the versions used before, however, multiplied by the yield

• sample is still the same, but does not require the argument n anymore. By default, this will now equal to a poissonian sampled n around the yield.

The SumPDF now does not strictly need fracs anymore: if all input PDFs are extended, the sum will be as well and use the (normalized) yields as fracs

The preferred way to create an extended PDf is to use extended=yield in the PDF constructor.

Alternatively, one can use extended_pdf = pdf.create_extended(yield). (Since this relies on copying the PDF, which may does not work for different reasons, there is also a pdf.set_yield(yield) method that sets the yield in-place. This won’t lead to ambiguities, as everything is supposed to work the same.)

yield_model = zfit.Parameter('yield_model', 10000, 0, 20000, step_size=10)
# model_ext = model.create_extended(yield_model)  # using an existing model

# when creating a new model 
model_ext = zfit.pdf.SumPDF([signal, comb_bkg, part_reco], [sig_frac, comb_bkg_frac], extended=yield_model)

alternatively, we can create the models as extended and sum them up

sig_yield = zfit.Parameter('sig_yield', 2000, 0, 10000, step_size=1)
sig_ext = signal.create_extended(sig_yield)

comb_bkg_yield = zfit.Parameter('comb_bkg_yield', 6000, 0, 10000, step_size=1)
comb_bkg_ext = comb_bkg.create_extended(comb_bkg_yield)

part_reco_yield = zfit.Parameter('part_reco_yield', 2000, 0, 10000, step_size=1)
part_reco.set_yield(
    part_reco_yield)  # unfortunately, `create_extended` does not work here. But no problem, it won't change anything.
part_reco_ext = part_reco

model_ext_sum = zfit.pdf.SumPDF([sig_ext, comb_bkg_ext, part_reco_ext])

Loss

A loss combines the model and the data, for example to build a likelihood. Furthermore, it can contain constraints, additions to the likelihood. Currently, if the Data has weights, these are automatically taken into account.

nll_gauss = zfit.loss.UnbinnedNLL(gauss, data_normal)

The loss has several attributes to be transparent to higher level libraries. We can calculate the value of it using value.

nll_gauss.value()

Notice that due to graph building, this will take significantly longer on the first run. Rerun the cell above and it will be way faster.

Furthermore, the loss also provides a possibility to calculate the gradients or, often used, the value and the gradients.

nll_gauss.gradient()

We can access the data and models (and possible constraints)

nll_gauss.model

nll_gauss.data

nll_gauss.constraints

Similar to the models, we can also get the parameters via get_params.

nll_gauss.get_params()

Extended loss

More interestingly, we can now build a loss for our composite sum model using the sampled data. Since we created an extended model, we can now also create an extended likelihood, taking into account a Poisson term to match the yield to the number of events.

nll = zfit.loss.ExtendedUnbinnedNLL(model_ext_sum, data)

nll.get_params()

Minimization

While a likelihood is interesting, we usually want to minimize it. Therefore, we can use the minimizers in zfit, most notably Minuit, a wrapper around the iminuit minimizer.

The philosophy is to create a minimizer instance that is mostly stateless, e.g. does not remember the position.

Given that iminuit provides us with a very reliable and stable minimizer, it is usually recommended to use this. Others are implemented as well and could easily be wrapped, however, the convergence is usually not as stable.

Minuit has a few options:

• tol: the Estimated Distance to Minimum (EDM) criteria for convergence (default 1e-3)

• verbosity: between 0 and 10, 5 is normal, 7 is verbose, 10 is maximum

• gradient: if True, uses the Minuit numerical gradient instead of the TensorFlow gradient. This is usually more stable for smaller fits; furthermore, the TensorFlow gradient can (experience based) sometimes be wrong. Using "zfit" uses the TensorFlow gradient.

minimizer = zfit.minimize.Minuit()
# minimizer = zfit.minimize.NLoptLBFGSV1()

Other minimizers

There are a couple of other minimizers, such as SciPy, NLopt that are wrapped and can be used.

They all provide a highly standardize interface, good initial values and non-ambiguous docstring, i.e. every method has its own docstring.

When to stop: a crucial problem is the question of convergence. iminuit has the EDM, a comparably reliable metric. Other minimizers in zfit use the same stopping criteria, to be able to compare results, but that is often very inefficienc: many minimizer do not use the Hessian, required for the EDM, internally, or do not provide it. This means, the Hessian has to be calculated separately (expensive!). Better stopping criteria ideas are welcome!

help(zfit.minimize.ScipyNewtonCGV1.__init__)

print(zfit.minimize.__all__)

For the minimization, we can call minimize, which takes a

• loss as we created above

• optionally: the parameters to minimize

By default, minimize uses all the free floating parameters (obtained with get_params). We can also explicitly specify which ones to use by giving them (or better, objects that depend on them) to minimize; note however that non-floating parameters, even if given explicitly to minimize won ‘t be minimized.

iminuit compatibility

As a sidenote, a zfit loss provides the same API as required by iminuit. Therefore, we can also directly invoke iminuit and use it.

The advantage can be that iminuit offers a few things out-of-the-box, like contour drawings, which zfit does not. Disadvantage is that the zfit worflow is exited.

import iminuit

iminuit_min = iminuit.Minuit(nll, nll.get_params())
# iminuit_min.migrad()  # this fails sometimes, as iminuit doesn't use the limits and is therefore unstable

result = minimizer.minimize(nll)

plot_comp_model(model_ext_sum, data)

Fit result

The result of every minimization is stored in a FitResult. This is the last stage of the zfit workflow and serves as the interface to other libraries. Its main purpose is to store the values of the fit, to reference to the objects that have been used and to perform (simple) uncertainty estimation.

Extended tutorial on the FitResult

print(result)

This gives an overview over the whole result. Often we’re mostly interested in the parameters and their values, which we can access with a params attribute.

print(result.params)

This is a dict which stores any knowledge about the parameters and can be accessed by the parameter (object) itself:

result.params[mu_sig]

‘value’ is the value at the minimum. To obtain other information about the minimization process, result contains more attributes:

• valid: whether the result is actually valid or not (most important flag to check)

• fmin: the function minimum

• edm: estimated distance to minimum

• info: contains a lot of information, especially the original information returned by a specific minimizer

• converged: if the fit converged

result.fmin

Info provides additional information on the minimization, which can be minimizer dependent. For example, if iminuit was used, the instance can be retrieved.

print(result.info['minuit'])  # the actual instance used in minimization

Estimating uncertainties

The FitResult has mainly two methods to estimate the uncertainty:

• a profile likelihood method (like MINOS)

• Hessian approximation of the likelihood (like HESSE)

When using Minuit, this uses (currently) its own implementation. However, zfit has its own implementation, which can be invoked by changing the method name.

Hesse also supports the corrections for weights.

We can explicitly specify which parameters to calculate, by default it does for all.

result.hesse(name='hesse')

# result.hesse(method='hesse_np', name='numpy_hesse')
# result.hesse(method='minuit_hesse', name='iminuit_hesse')

We get the result directly returned. This is also added to result.params for each parameter and is nicely displayed with an added column

print(result.params)

errors, new_result = result.errors(params=[sig_yield, part_reco_yield, mu_sig],
                                   name='errors')  # just using three, but we could also omit argument -> all used

errorszfit, new_result_zfit = result.errors(params=[sig_yield, part_reco_yield, mu_sig], method='zfit_errors', name='errorszfit')

print(errors)

print(result)

result.update_params()
result.hesse()
result.errors()

loss = result.loss

poi = sig_yield
poimin = result.params[poi]['value']
poierr = result.params[poi]['hesse']['error']
loss.value(params={poi: poimin + 1 * poierr}) - loss.value(params={poi: poimin})

from typing import Dict, List, Optional, Union

import matplotlib.pyplot as plt
import numpy as np
import zfit


def create_impact_plot(
        result: zfit.result.FitResult,
        parameters_of_interest: Optional[list[zfit.Parameter]] = None,
        sort_by_impact: bool = True,
        figsize: tuple = (10, 6),
        title: str = "Parameter Impact Plot",
        color: str = "royalblue",
        error_bar_color: str = "black",
        save_path: Optional[str] = None
) -> plt.Figure:
    """
    Create an impact plot from a zfit FitResult to visualize how each parameter
    affects the fit value and uncertainty.

    Parameters:
    -----------
    result : zfit.minimize.FitResult
        The result object from a zfit minimization
    parameters_of_interest : List[zfit.Parameter], optional
        List of parameters to include in the plot. If None, all parameters are used.
    sort_by_impact : bool, default=True
        Whether to sort parameters by their impact magnitude
    figsize : tuple, default=(10, 6)
        Size of the figure to create
    title : str, default="Parameter Impact Plot"
        Title of the plot
    color : str, default="royalblue"
        Color for the parameter impact bars
    error_bar_color : str, default="black"
        Color for the error bars
    save_path : str, optional
        If provided, save the figure to this path

    Returns:
    --------
    plt.Figure
        The created figure
    """
    # Get parameters to plot
    if parameters_of_interest is None:
        parameters = result.params
    else:
        parameters = parameters_of_interest

    # Extract parameter values and errors
    param_names = [param.name for param in parameters]
    param_values = [result.params[param]["value"] for param in param_names]
    param_errors = [result.params[param]["hesse"]['error'] for param in param_names]

    # Calculate normalized impact (parameter value / error)
    impacts = np.array([val / err if err != 0 else 0 for val, err in zip(param_values, param_errors)])

    # Sort by impact magnitude if requested
    if sort_by_impact:
        # Get sorting indices by absolute impact
        sorted_indices = np.argsort(np.abs(impacts))[::-1]
        param_names = [param_names[i] for i in sorted_indices]
        param_values = [param_values[i] for i in sorted_indices]
        param_errors = [param_errors[i] for i in sorted_indices]
        impacts = impacts[sorted_indices]

    # Create the figure
    fig, ax = plt.subplots(figsize=figsize)

    # Plot parameter impacts
    y_pos = np.arange(len(param_names))
    bars = ax.barh(y_pos, impacts, align='center', color=color, alpha=0.7)

    # Add error bars (showing ±1)
    ax.axvline(x=0, color='grey', linestyle='-', alpha=0.3)
    ax.axvline(x=1, color='grey', linestyle='--', alpha=0.2)
    ax.axvline(x=-1, color='grey', linestyle='--', alpha=0.2)

    # Add labels and formatting
    ax.set_yticks(y_pos)
    ax.set_yticklabels(param_names)
    ax.invert_yaxis()  # Labels read top-to-bottom
    ax.set_xlabel('Impact (parameter value / uncertainty)')
    ax.set_title(title)

    # Add values as text
    for i, (impact, value, error) in enumerate(zip(impacts, param_values, param_errors)):
        ax.text(impact + np.sign(impact) * 0.1, i,
                f"{value:.3f} ± {error:.3f}",
                va='center', fontsize=9)

    # Set reasonable x limits
    max_impact = max(abs(np.max(impacts)), abs(np.min(impacts)))
    ax.set_xlim(-max(max_impact + 0.5, 2), max(max_impact + 0.5, 2))

    # Add grid lines
    ax.grid(axis='x', linestyle='--', alpha=0.3)

    plt.tight_layout()

    # Save if requested
    if save_path:
        plt.savefig(save_path, dpi=300, bbox_inches='tight')

    return fig


# Example usage:
"""
# After performing a fit with zfit:
import zfit

# Define model and perform fit
# ...
# result = minimizer.minimize(loss)

# Create impact plot
fig = create_impact_plot(result)
plt.show()
"""

pois = [sig_yield, part_reco_yield, mu_sig]
pois = nll.get_params()

result.hesse()
# result.errors(params=pois)

poi.set_value(poimin + poierr)
poi.floating = False
print(f"{poi} set to poimin={poimin} + poierr={poierr}")

# for p in nll.get_params(floating=None):
#     p.floating = True

from typing import Dict, List, Optional, Tuple, Union

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import zfit


def create_cms_impact_plot(
        result: zfit.result.FitResult,
        param_of_interest: zfit.Parameter,
        nuisance_params: list[zfit.Parameter],
        pulls: Optional[dict[str, float]] = None,
        impacts_plus: Optional[dict[str, float]] = None,
        impacts_minus: Optional[dict[str, float]] = None,
        poi_best_fit: Optional[float] = None,
        poi_uncertainty: Optional[tuple[float, float]] = None,
        figsize: tuple = (12, 10),
        title: str = "CMS Internal",
        sort_by_impact: bool = True,
        save_path: Optional[str] = None
) -> plt.Figure:
    """
    Create a CMS-style impact plot showing nuisance parameter pulls and their impact on a parameter of interest.

    Parameters:
    -----------
    result : zfit.minimize.FitResult
        The result object from a zfit minimization
    param_of_interest : zfit.Parameter
        The parameter of interest (POI) for which to show impacts
    nuisance_params : List[zfit.Parameter]
        List of nuisance parameters to include in the plot
    pulls : Dict[str, float], optional
        Dictionary of pre-calculated pulls for each nuisance parameter.
        If None, they will be calculated from the result.
    impacts_plus : Dict[str, float], optional
        Dictionary of pre-calculated +1σ impacts on the POI for each nuisance parameter
    impacts_minus : Dict[str, float], optional
        Dictionary of pre-calculated -1σ impacts on the POI for each nuisance parameter
    poi_best_fit : float, optional
        Best fit value of the parameter of interest. If None, extracted from result.
    poi_uncertainty : Tuple[float, float], optional
        Asymmetric uncertainty on the POI (upper, lower). If None, extracted from result.
    figsize : tuple, default=(12, 10)
        Size of the figure to create
    title : str, default="CMS Internal"
        Title of the plot
    sort_by_impact : bool, default=True
        Whether to sort parameters by their impact magnitude
    save_path : str, optional
        If provided, save the figure to this path

    Returns:
    --------
    plt.Figure
        The created figure
    """
    # Get parameter names
    param_names = [param.name for param in nuisance_params]
    n_params = len(param_names)

    # If pulls are not provided, extract them from the result
    # In a real implementation, you would need to get these from the fit result
    if pulls is None:
        pulls = {}
        for name in param_names:
            # This is a simplification - in real HEP analyses pulls would be calculated differently
            param_value = result.params[name]["value"]
            param_error = result.params[name]["hesse"]['error']
            pulls[name] = param_value / param_error if param_error != 0 else 0

    # If impacts are not provided, you'd need to calculate them
    # This would typically involve refitting with fixed parameters at ±1σ
    # Here we'll use random values for demonstration if not provided
    if impacts_plus is None or impacts_minus is None:
        np.random.seed(42)  # For reproducibility
        impacts_plus = {}
        impacts_minus = {}
        for name in param_names:
            impacts_plus[name] = np.random.uniform(0, 0.15)
            impacts_minus[name] = -np.random.uniform(0, 0.15)

    # Calculate total impacts for sorting
    if sort_by_impact:
        total_impacts = {name: max(abs(impacts_plus[name]), abs(impacts_minus[name])) for name in param_names}
        # Sort parameter names by impact magnitude
        param_names = sorted(param_names, key=lambda x: total_impacts[x], reverse=True)

    # Extract the POI best fit and uncertainty if not provided
    poi_name = param_of_interest.name
    if poi_best_fit is None:
        poi_best_fit = result.params[poi_name]["value"]

    if poi_uncertainty is None:
        poi_error = result.params[poi_name]["hesse"]['error']
        poi_uncertainty = (poi_error, poi_error)  # Symmetric uncertainty as default

    # Create figure with gridspec for the two panels
    fig = plt.figure(figsize=figsize)
    gs = gridspec.GridSpec(1, 2, width_ratios=[3, 1])

    # Left panel - pulls
    ax_pulls = plt.subplot(gs[0])

    # Right panel - impacts
    ax_impacts = plt.subplot(gs[1])

    # Initialize y positions
    y_pos = np.arange(n_params)

    # Plot pulls - grey dots with error bars and blue x markers
    for i, name in enumerate(param_names):
        pull = pulls[name]
        # Plot grey dot with error bar for the fit result
        ax_pulls.errorbar(pull, y_pos[i], xerr=1, fmt='o', color='grey', markersize=6)
        # Plot blue x for the pull
        ax_pulls.plot(pull, y_pos[i], 'x', color='blue', markersize=6)

    # Plot impacts
    for i, name in enumerate(param_names):
        plus_impact = impacts_plus[name]
        minus_impact = impacts_minus[name]

        # +1σ impact (pink/red)
        ax_impacts.barh(y_pos[i], plus_impact, height=0.6, color='#d6a6a6', alpha=0.7)

        # -1σ impact (blue)
        ax_impacts.barh(y_pos[i], minus_impact, height=0.6, color='#a6a6d6', alpha=0.7)

    # Set up the left panel (pulls)
    ax_pulls.set_yticks(y_pos)
    ax_pulls.set_yticklabels(param_names)
    ax_pulls.set_xlim(-2.5, 2.5)
    ax_pulls.axvline(x=0, color='grey', linestyle=':', alpha=0.7)
    ax_pulls.axvline(x=-1, color='grey', linestyle=':', alpha=0.5)
    ax_pulls.axvline(x=1, color='grey', linestyle=':', alpha=0.5)
    ax_pulls.set_xlabel(r'$(\hat{\theta}-\theta_0)/\sigma_i$')

    # Number the parameters
    for i, _ in enumerate(param_names):
        ax_pulls.text(-2.4, y_pos[i], f"{i+1}", va='center', ha='center', fontsize=9)

    # Set up the right panel (impacts)
    ax_impacts.set_yticks(y_pos)
    ax_impacts.set_yticklabels([])  # No labels on the right panel
    ax_impacts.axvline(x=0, color='grey', linestyle=':', alpha=0.7)
    max_impact = max(max(abs(v) for v in impacts_plus.values()),
                     max(abs(v) for v in impacts_minus.values()))
    impact_limit = max(0.15, max_impact * 1.2)  # Ensure we have reasonable limits
    ax_impacts.set_xlim(-impact_limit, impact_limit)
    ax_impacts.set_xlabel(r'$\Delta\hat{r}$')

    # Add POI best fit value and uncertainty at the top of the plot
    fig.text(0.7, 0.92, rf'$\hat{{r}} = {poi_best_fit:.1f}^{{+{poi_uncertainty[0]:.1f}}}_{{-{poi_uncertainty[1]:.1f}}}$',
             fontsize=14, ha='center')

    # Add legend to the bottom
    fig.text(0.07, 0.03, r'$\bullet$ Fit', fontsize=12, ha='left')
    fig.text(0.07, 0.01, r'$\times$ Pull', fontsize=12, ha='left', color='blue')
    fig.text(0.25, 0.03, r' +1σ Impact', fontsize=12, ha='left', color='#d6a6a6')
    fig.text(0.25, 0.01, r' -1σ Impact', fontsize=12, ha='left', color='#a6a6d6')

    # Add title
    fig.text(0.45, 0.94, title, fontsize=16, ha='center')

    # Adjust layout
    plt.subplots_adjust(wspace=0.02, left=0.3, right=0.95, bottom=0.1, top=0.9)

    # Save if requested
    if save_path:
        plt.savefig(save_path, dpi=300, bbox_inches='tight')

    return fig

# Example usage with random data for demonstration


def example_usage():
    """Demonstrate the function with random data"""
    import zfit
    from zfit import Parameter

    # Create dummy parameter of interest
    poi = Parameter("r", 1.3)

    # Create dummy nuisance parameters
    nuisance_params = [
        Parameter(f"CMS_scale_t_tautau_8TeV", 0),
        Parameter(f"CMS_eff_t_tt_8TeV", 0),
        Parameter(f"CMS_htt_tt_tauTau_1jet_high_highhiggs_8TeV_ZTT_bin_12", 0),
        Parameter(f"CMS_htt_tt_tauTau_1jet_high_highhiggs_8TeV_ZTT_bin_11", 0),
        Parameter(f"QCDscale_ggH1in", 0),
        Parameter(f"CMS_htt_QCDSyst_tauTau_vbf_8TeV", 0),
        Parameter(f"pdf_gg", 0),
        Parameter(f"CMS_htt_tt_tauTau_1jet_high_highhiggs_8TeV_ZTT_bin_13", 0),
        Parameter(f"CMS_htt_QCDSyst_tauTau_1jet_high_highhiggs_8TeV", 0),
        Parameter(f"CMS_htt_tt_tauTau_vbf_8TeV_ZTT_bin_6", 0),
    ]

    # Create dummy result
    class DummyResult:
        def __init__(self):
            self.params = {}
            for param in [poi] + nuisance_params:
                self.params[param.name] = {"value": np.random.normal(0, 1), "error": 1.0}
            # Set POI value and error
            self.params["r"] = {"value": 1.3, "error": 0.75}

    result = DummyResult()

    # Create pulls dictionary
    np.random.seed(42)
    pulls = {}
    for param in nuisance_params:
        pulls[param.name] = np.random.normal(0, 1)

    # Create impacts dictionaries
    impacts_plus = {}
    impacts_minus = {}
    for param in nuisance_params:
        impacts_plus[param.name] = np.random.uniform(0, 0.1)
        impacts_minus[param.name] = -np.random.uniform(0, 0.1)

    # Create the plot
    fig = create_cms_impact_plot(
        result=result,
        param_of_interest=poi,
        nuisance_params=nuisance_params,
        pulls=pulls,
        impacts_plus=impacts_plus,
        impacts_minus=impacts_minus,
        poi_best_fit=1.3,
        poi_uncertainty=(0.8, 0.7),
        title="CMS Internal"
    )

    plt.savefig("cms_impact_plot_example.png", dpi=300)
    plt.show()

# if __name__ == "__main__":
#     example_usage()

create_cms_impact_plot(result, param_of_interest=sig_yield, nuisance_params=nll.get_params() - {sig_yield}, sort_by_impact=True)

res = minimizer.minimize(loss).update_params()
# res.hesse()
# res.errors()
res

res.fmin- result.fmin

res.update_params()

What is ‘new_result’?

When profiling a likelihood, such as done in the algorithm used in errors, a new minimum can be found. If this is the case, this new minimum will be returned, otherwise new_result is None. Furthermore, the current result would be rendered invalid by setting the flag valid to False. Note: this behavior only applies to the zfit internal error estimator.

We can also access the covariance matrix of the parameters

result.covariance()

Bonus: serialization

Human-readable serialization is experimentally supported in zfit. zfit aims to follow the HS3 standard, which is currently in development.

Extended tutorial on serialization

model.to_dict()

gauss.to_yaml()

zfit.hs3.dumps(nll)

End of zfit

This is where zfit finishes and other libraries take over.

Beginning of hepstats

hepstats is a library containing statistical tools and utilities for high energy physics. In particular you do statistical inferences using the models and likelhoods function constructed in zfit.

Short example: let’s compute for instance a confidence interval at 68 % confidence level on the mean of the gaussian defined above.

from hepstats.hypotests import ConfidenceInterval
from hepstats.hypotests.calculators import AsymptoticCalculator
from hepstats.hypotests.parameters import POIarray

calculator = AsymptoticCalculator(input=result, minimizer=minimizer)

value = result.params[mu_sig]["value"]
error = result.params[mu_sig]["hesse"]["error"]

mean_scan = POIarray(mu_sig, np.linspace(value - 1.5 * error, value + 1.5 * error, 10))

ci = ConfidenceInterval(calculator, mean_scan)

ci.interval()

from utils import one_minus_cl_plot

ax = one_minus_cl_plot(ci)
ax.set_xlabel("mean")

There will be more of hepstats later.