# Examples¶

## Composing PDFs¶

#  Copyright (c) 2020 zfit

import zfit

# create space
obs = zfit.Space("x", limits=(-10, 10))

# parameters
mu = zfit.Parameter("mu", 1., -4, 6)
sigma = zfit.Parameter("sigma", 1., 0.1, 10)
lambd = zfit.Parameter("lambda", -1., -5., 0)
frac = zfit.Parameter("fraction", 0.5, 0., 1.)

# pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)

sum_pdf = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)


## Custom 3D Functor¶

#  Copyright (c) 2019 zfit

# Functors are PDfs and Functions that depend on other PDFs or Functions. They can be used to define
# in a custom way combinations of PDFs or wrapping a single PDF.
# An example would be to create the sum of two PDFs. Of course this is already implemented in zfit
# as the functor SumPDF([pdf1, pdf2], fracs=...). For advanced uses, you
# can define your own functor as demonstrated below.

import zfit
import numpy as np

class CombinePolynomials(zfit.pdf.BaseFunctor):
"""Example of a functor pdf that adds three polynomials.

DEMONSTRATION PURPOSE ONLY, DO **NOT** USE IN REAL CASE.
"""

def __init__(self, angle1, angle2, angle3, name="3dPolynomial"):
pdfs = [angle1, angle2, angle3]
super().__init__(pdfs=pdfs, name=name)

_N_OBS = 3

def _unnormalized_pdf(self, x):
cosangle1_data, cosangle2_data, angle3_data = x.unstack_x()
cosangle1, cosangle2, angle3 = self.pdfs

pdf = cosangle1.pdf(cosangle1_data) + cosangle2.pdf(cosangle2_data) + angle3.pdf(angle3_data)
return pdf

def create_angular():
# create parameters
# c00 = zfit.Parameter(...)
# ... and so on

cosangle1_space = zfit.Space('cos angle 1', limits=(-1, 1))
cosangle2_space = zfit.Space('cos angle 2', limits=(-1, 1))
angle3_space = zfit.Space('angle 3', limits=(-np.pi, np.pi))

# this part could also be moved inside the __init__, but then the init would need to take all
# the coefficients and the spaces as arguments
cosangle1_pdf = zfit.pdf.Chebyshev(obs=cosangle1_space, coeffs=[c00, c01, c02])
cosangle2_pdf = zfit.pdf.Chebyshev(obs=cosangle2_space, coeffs=[c10, c11, c12])
angle3_pdf = zfit.pdf.Chebyshev(obs=angle3_space, coeffs=[c20, c21, c22])

return CombinePolynomials(angle1=cosangle1_pdf, angle2=cosangle2_pdf, angle3=angle3_pdf)


## Custom PDF simple¶

#  Copyright (c) 2020 zfit

import zfit
from zfit import z

class CustomPDF(zfit.pdf.ZPDF):
"""1-dimensional PDF implementing the exp(alpha * x) shape."""
_PARAMS = ['alpha']  # specify which parameters to take

def _unnormalized_pdf(self, x):  # implement function
data = z.unstack_x(x)
alpha = self.params['alpha']

return z.exp(alpha * data)

obs = zfit.Space("obs1", limits=(-4, 4))

custom_pdf = CustomPDF(obs=obs, alpha=0.2)

integral = custom_pdf.integrate(limits=(-1, 2))
sample = custom_pdf.sample(n=1000)
prob = custom_pdf.pdf(sample)


#  Copyright (c) 2020 zfit

import tensorflow as tf

import zfit
from zfit import z

class CustomPDF2D(zfit.pdf.BasePDF):
"""My custom, 2 dimensional pdf. The axes are: Energy, Momentum.

"""

def __init__(self, param1, param2, param3, obs, name="CustomPDF", ):
# we can now do complicated stuff here if needed
# only thing: we have to specify explicitly here what is which parameter
params = {'super_param': param1,  # we can change/compose etc parameters
'param2': param2, 'param3': param3}
super().__init__(obs, params, name=name)

def _unnormalized_pdf(self, x):
energy, momentum = x.unstack_x()
param1 = self.params['super_param']
param2 = self.params['param2']
param3 = self.params['param3']

# just a fantasy function
probs = param1 * tf.cos(energy ** 2) + tf.math.log(param2 * momentum ** 2) + param3
return probs

# define the integral function
def integral_full(x, limits, norm_range, params, model):
lower, upper = limits.limit1d
param1 = params['super_param']
param2 = params['param2']
param3 = params['param3']

lower = z.convert_to_tensor(lower)
upper = z.convert_to_tensor(upper)

# calculate the integral here, dummy integral, wrong!
integral = param1 * param2 * param3 + z.reduce_sum([lower, upper])
return integral

# define the space over which it is defined. Here, we use the axes
lower_full = (-10, zfit.Space.ANY_LOWER)
upper_full = (10, zfit.Space.ANY_UPPER)
integral_full_limits = zfit.Space(axes=(0, 1),
limits=(lower_full, upper_full))

CustomPDF2D.register_analytic_integral(func=integral_full,
limits=integral_full_limits)

# define the partial integral function
def integral_axis1(x, limits, norm_range, params, model):
data_0 = x.unstack_x()  # data from axis 0

param1 = params['super_param']
param2 = params['param2']
param3 = params['param3']

lower, upper = limits.limit1d
lower = z.convert_to_tensor(lower)  # the limits are now 1-D, for axis 1
upper = z.convert_to_tensor(upper)

# calculate the integral here, dummy integral
integral = data_0 ** 2 * param1 * param2 * param3 + z.reduce_sum([lower, upper])
return integral

# define the space over which it is defined. Here, we use the axes
lower_axis1 = ((-5,),)
upper_axis1 = ((zfit.Space.ANY_UPPER,),)
integral_axis1_limits = zfit.Space(axes=(1,),
limits=(lower_axis1, upper_axis1))

CustomPDF2D.register_analytic_integral(func=integral_axis1,
limits=integral_axis1_limits)

if __name__ == '__main__':
obs = zfit.Space('obs1', (-10, 10)) * zfit.Space('obs2', (-3, 5))
pdf = CustomPDF2D(1, 2, 3, obs=obs)
sample = pdf.sample(n=1000)


## Custom 3D PDF¶

#  Copyright (c) 2019 zfit

import zfit
from zfit import z
import tensorflow as tf

# IMPORTANT! The communication of which axis corresponds to which data point happens here. So the user knows now that
# he should create this pdf with a space in the obs (x, y, z).

class CustomPDF(zfit.pdf.ZPDF):
"""3-dimensional PDF calculating doing something fancy. Takes x, y, z as data."""
_PARAMS = ['alpha', 'beta']  # specify which parameters to take
_N_OBS = 3

def _unnormalized_pdf(self, x):  # implement function
x, y, z = x.unstack_x()
alpha = self.params['alpha']
beta = self.params['beta']
x_new = tf.math.cos(alpha) * x
y_new = tf.math.sinh(beta) * y
z_new = z + 4.2
return x_new ** 2 + y_new ** 2 + z_new ** 2

xobs = zfit.Space('xobs', (-4, 4))
yobs = zfit.Space('yobs', (-3, 3))
zobs = zfit.Space('z', (-2, 2))
obs = xobs * yobs * zobs

alpha = zfit.Parameter("alpha", 0.2)  # floating
beta = zfit.Parameter("beta", 0.4, floating=False)  # non-floating
custom_pdf = CustomPDF(obs=obs, alpha=alpha, beta=beta)

integral = custom_pdf.integrate(limits=obs)  # = 1 since normalized
sample = custom_pdf.sample(n=1000)  # DO NOT USE THIS FOR TOYS!
prob = custom_pdf.pdf(sample)  # DO NOT USE THIS FOR TOYS!

integral_np, sample_np, prob_np = [integral.numpy(), sample.numpy(), prob.numpy()]


## Multidimensional preprocess fit¶

#  Copyright (c) 2020 zfit

import numpy as np

import zfit

# create space
xobs = zfit.Space('xobs', (-4, 4))
yobs = zfit.Space('yobs', (-3, 5))
zobs = zfit.Space('z', (-2, 4))
obs = xobs * yobs * zobs

# parameters
mu1 = zfit.Parameter("mu1", 1., -4, 6)
mu23 = zfit.Parameter("mu_shared", 1., -4, 6)
sigma12 = zfit.Parameter("sigma_shared", 1., 0.1, 10)
sigma3 = zfit.Parameter("sigma3", 1., 0.1, 10)

# model building, pdf creation
gauss_x = zfit.pdf.Gauss(mu=mu1, sigma=sigma12, obs=xobs)
gauss_y = zfit.pdf.Gauss(mu=mu23, sigma=sigma12, obs=yobs)
gauss_z = zfit.pdf.Gauss(mu=mu23, sigma=sigma3, obs=zobs)

product_gauss = zfit.pdf.ProductPDF([gauss_x, gauss_y, gauss_z])

# OR create directly from your 3 dimensional pdf as
# model = MyPDF(obs=obs, param1=..., param2,...)

# data
normal_np = np.random.normal(loc=[2., 2.5, 2.5], scale=[3., 3, 1.5], size=(10000, 3))
data_raw = zfit.Data.from_numpy(obs=obs, array=normal_np)  # or from anywhere else, e.g. root

df = data_raw.to_pandas()
# preprocessing here, rename things. Match column names with the observable names "xobs", "yobs", "z" (they have to be
# contained, more columns in the df is not a problem)
data = zfit.Data.from_pandas(df, obs=obs)

# create NLL
nll = zfit.loss.UnbinnedNLL(model=product_gauss, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll)
print(result.params)

# do the error calculations, here with minos
param_errors, _ = result.errors()
print(param_errors)


## Plot signal background¶

#  Copyright (c) 2020 zfit

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

# create space
obs = zfit.Space("x", limits=(-10, 10))

# parameters
mu = zfit.Parameter("mu", 1., -4, 6)
sigma = zfit.Parameter("sigma", 1., 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)
frac = zfit.Parameter("fraction", 0.3, 0, 1)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)
model = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)

# data
n_sample = 10000

exp_data = exponential.sample(n=n_sample * (1 - frac)).numpy()

gauss_data = gauss.sample(n=n_sample * frac).numpy()

data = model.create_sampler(n_sample, limits=obs)
data.resample()

mu.set_value(0.5)
sigma.set_value(1.2)
lambd.set_value(-0.05)
frac.set_value(0.07)

# plot the data
data_np = data[:, 0].numpy()
color = 'black'
n_bins = 50

linewidth = 2.5
plot_scaling = n_sample / n_bins * obs.area()

x = np.linspace(-10, 10, 1000)

# plot the pdf BEFORE fitting
plt.figure()
plt.title("Before fitting")
# plot the data
plt.hist(data_np, color=color, bins=n_bins, histtype="stepfilled", alpha=0.1)
plt.hist(data_np, color=color, bins=n_bins, histtype="step")
# plot the pdfs
y = model.pdf(x).numpy()
y_gauss = (gauss.pdf(x) * frac).numpy()  # notice the frac!
y_exp = (exponential.pdf(x) * (1 - frac)).numpy()  # notice the frac!

plt.plot(x, y * plot_scaling, label="Sum - Model", linewidth=linewidth * 2)
plt.plot(x, y_gauss * plot_scaling, '--', label="Gauss - Signal", linewidth=linewidth)
plt.plot(x, y_exp * plot_scaling, '--', label="Exponential - Background", linewidth=linewidth)
plt.xlabel("Physical observable")
plt.legend()

# create NLL
nll = zfit.loss.UnbinnedNLL(model=model, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll)

# do the error calculations, here with minos
param_errors, _ = result.errors()

plt.figure()
plt.title("After fitting")
# plot the data
plt.hist(data_np, color=color, bins=n_bins, histtype="stepfilled", alpha=0.1)
plt.hist(data_np, color=color, bins=n_bins, histtype="step")

y = model.pdf(x).numpy()  # rerun now after the fitting
y_gauss = (gauss.pdf(x) * frac).numpy()
y_exp = (exponential.pdf(x) * (1 - frac)).numpy()

plt.plot(x, y * plot_scaling, label="Sum - Model", linewidth=linewidth * 2)
plt.plot(x, y_gauss * plot_scaling, '--', label="Gauss - Signal", linewidth=linewidth)
plt.plot(x, y_exp * plot_scaling, '--', label="Exponential - Background", linewidth=linewidth)
plt.xlabel("Physical observable")
plt.legend()

plt.show()


## Simple fit¶

#  Copyright (c) 2020 zfit

import numpy as np
import zfit

# create space
obs = zfit.Space("x", limits=(-2, 3))

# parameters
mu = zfit.Parameter("mu", 1.2, -4, 6)
sigma = zfit.Parameter("sigma", 1.3, 0.1, 10)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)

# data
normal_np = np.random.normal(loc=2., scale=3., size=10000)
data = zfit.Data.from_numpy(obs=obs, array=normal_np)

# create NLL
nll = zfit.loss.UnbinnedNLL(model=gauss, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll)

# do the error calculations, here with minos
param_errors, _ = result.errors()


## Signal-background-mass fit¶

#  Copyright (c) 2020 zfit

import zfit

# create space
obs = zfit.Space("x", limits=(-10, 10))

# parameters
mu = zfit.Parameter("mu", 1., -4, 6)
sigma = zfit.Parameter("sigma", 1., 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)
frac = zfit.Parameter("fraction", 0.3, 0, 1)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)
model = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)

# data
n_sample = 10000

data = model.create_sampler(n_sample, limits=obs)
data.resample()

# set the values to a start value for the fit
mu.set_value(0.5)
sigma.set_value(1.2)
lambd.set_value(-0.05)
frac.set_value(0.07)

# create NLL
nll = zfit.loss.UnbinnedNLL(model=model, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll)
print(result)

# do the error calculations, here with minos
param_hesse = result.hesse()
param_errors, new_result = result.errors()
print(result.params)


## Signal-background-mass fit extended¶

#  Copyright (c) 2020 zfit

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

n_bins = 50

# create space
obs = zfit.Space("x", limits=(-10, 10))

# parameters
mu = zfit.Parameter("mu", 1., -4, 6)
sigma = zfit.Parameter("sigma", 1., 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)

n_bkg = zfit.Parameter('n_bkg', 20000)
n_sig = zfit.Parameter('n_sig', 1000)
gauss_extended = gauss.create_extended(n_sig)
exp_extended = exponential.create_extended(n_bkg)
model = zfit.pdf.SumPDF([gauss_extended, exp_extended])

# data
# n_sample = 10000

data = model.create_sampler(n=21200)
data.resample()

# set the values to a start value for the fit
mu.set_value(0.5)
sigma.set_value(1.2)
lambd.set_value(-0.05)

# create NLL
nll = zfit.loss.ExtendedUnbinnedNLL(model=model, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit(verbosity=7)
result = minimizer.minimize(nll)
print(result.params)
# do the error calculations, here with hesse, than with minos
param_hesse = result.hesse()
param_errors, _ = result.errors()  # this returns a new FitResult if a new minimum was found
print(result.valid)  # check if the result is still valid


## Simultanous fit¶

#  Copyright (c) 2019 zfit

# create space
import zfit
import numpy as np

obs = zfit.Space("x", limits=(-10, 10))

# parameters
mu_shared = zfit.Parameter("mu_shared", 1., -4, 6)  # mu is a shared parameter
sigma1 = zfit.Parameter("sigma_one", 1., 0.1, 10)
sigma2 = zfit.Parameter("sigma_two", 1., 0.1, 10)

# model building, pdf creation
gauss1 = zfit.pdf.Gauss(mu=mu_shared, sigma=sigma1, obs=obs)
gauss2 = zfit.pdf.Gauss(mu=mu_shared, sigma=sigma2, obs=obs)

# data
normal_np = np.random.normal(loc=2., scale=3., size=10000)
data1 = zfit.Data.from_numpy(obs=obs, array=normal_np)

# data
normal_np = np.random.normal(loc=2., scale=4., size=10000)
data2 = zfit.Data.from_numpy(obs=obs, array=normal_np)

# create simultaenous loss, two possibilities
nll_simultaneous = zfit.loss.UnbinnedNLL(model=[gauss1, gauss2],
data=[data1, data2])
# OR, equivalently
nll1 = zfit.loss.UnbinnedNLL(model=gauss1, data=data1)
nll2 = zfit.loss.UnbinnedNLL(model=gauss2, data=data2)
nll_simultaneous2 = nll1 + nll2

minimizer = zfit.minimize.Minuit()
minimizer.minimize(nll_simultaneous2)


## Spaces¶

#  Copyright (c) 2020 zfit

import zfit

# Addition of limit with the same observables
simple_limit1 = zfit.Space(obs='obs1', limits=(-5, 1))
simple_limit2 = zfit.Space(obs='obs1', limits=(3, 7.5))

# OR equivalently

# multiplication of limits with different observables
first_limit_lower = (-5, 6)
first_limit_upper = (5, 10)

second_limit_lower = (7, 12)
second_limit_upper = (9, 15)

space1 = zfit.Space(obs=['obs1', 'obs2'], limits=(first_limit_lower, first_limit_upper))
space2 = zfit.Space(obs=['obs3', 'obs4'], limits=(second_limit_lower, second_limit_upper))

space4 = space1 * space2

assert space4.obs == ('obs1', 'obs2', 'obs3', 'obs4')
assert space4.n_obs == 4