<jupyter_start><jupyter_text>Example of upper limit computation.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import lauztat
from lauztat.parameters import POI
from lauztat.hypotests import UpperLimit
from lauztat.calculators import AsymptoticCalculator
from lauztat.config import Config
import tensorflow as tf
import zfit
from zfit import ztf
from zfit.core.loss import ExtendedUnbinnedNLL, UnbinnedNLL
from zfit.minimizers.minimizer_minuit import MinuitMinimizer<jupyter_output><empty_output><jupyter_text>ScenarioLet's assume you want to find a resonnance between 0.1 and 3.0 GeV/c$^2$ of invariant mass of a system.You have a signal region and a background only region:<jupyter_code>bounds = (0.1, 3.0)

#background only
np.random.seed(10)
tau = 2.0
beta = 1./tau
data = np.random.exponential(beta, 1000)

data_bkg = data[(data > bounds[0]) & (data < bounds[1])]
plt.hist(data_bkg, bins=100, histtype='step');
#signal + background

np.random.seed(0)
tau = -2.0
beta = -1/tau
data = np.random.exponential(beta, 300)
peak = np.random.normal(1.2, 0.1, 10)
data = np.concatenate((data,peak))
data = data[(data > bounds[0]) & (data < bounds[1])]<jupyter_output><empty_output><jupyter_text>Background fit, in background only region:The background model is an exponential.<jupyter_code>obs = zfit.Space('x', limits=bounds)
lambda_ = zfit.Parameter("lambda_",-2.0, -4.0, -0.5)
model_bkg = zfit.pdf.Exponential(obs=obs, lambda_=lambda_)
data_bkg_ = zfit.data.Data.from_numpy(obs=obs, array=data_bkg)
nll_bkg = UnbinnedNLL(model=[model_bkg], data=[data_bkg_], fit_range=[obs]) 
minimizer = MinuitMinimizer()
min_bkg = minimizer.minimize(loss=nll_bkg)
min_bkg.params[lambda_]
min_bkg.error();
lambda_mu = min_bkg.params[lambda_]["value"]
lambda_sigma = 0.0748696
print(lambda_mu, lambda_sigma)<jupyter_output>-2.022148383099551 0.0748696<jupyter_text>Signal + background fit in signal region:The total model is a sum of a gaussian and an exponential. The tau parameter of the exponential will be constrained to the value obtained in the background only region.<jupyter_code>mean = zfit.Parameter("mean", 1.2, 0.1, 2., floating=False)
sigma = zfit.Parameter("sigma", 0.1, floating=False)
Nsig = zfit.Parameter("Nsig", 1., -20., len(data))
Nbkg = zfit.Parameter("Nbkg", len(data), 0., len(data)*1.1)
signal = Nsig * zfit.pdf.Gauss(obs=obs, mu=mean, sigma=sigma)
background =  Nbkg * model_bkg
tot_model = signal + background
constraint = zfit.constraint.nll_gaussian(params=[lambda_], mu=[lambda_mu], sigma=[lambda_sigma])
data_ = zfit.data.Data.from_numpy(obs=obs, array=data)
nll = ExtendedUnbinnedNLL(model=[tot_model], data=[data_], fit_range=[obs], constraints=constraint) 
minimum = minimizer.minimize(loss=nll)
def plotfitresult(pdf, bounds, nbins, data):
    x = np.linspace(*bounds, num=1000)
    pdf = zfit.run(tot_model.pdf(x, norm_range=bounds) * tot_model.get_yield())
    _ = plt.plot(x, ((bounds[1] - bounds[0])/nbins)*(pdf), "-r")  
nbins = 100
plt.hist(data, bins=nbins, histtype='step', range=bounds);
plotfitresult(tot_model, bounds, nbins, data)
plt.xlabel("m [GeV/c$^2$]")
plt.ylabel("number of events")
plt.savefig("fit_upper_limit_ex.png")<jupyter_output><empty_output><jupyter_text>Upper limit:<jupyter_code>def lossbuilder(model, data):
    constraint = zfit.constraint.nll_gaussian(params=[lambda_], mu=[lambda_mu], sigma=[lambda_sigma])
    loss = ExtendedUnbinnedNLL(model=model, data=data, fit_range=[obs], constraints=constraint) 
    return loss

config = Config(tot_model, data_, lossbuilder, MinuitMinimizer())
config.bestfit
calc = AsymptoticCalculator(config)
poinull = POI(Nsig, value=np.linspace(0.0, 25, 20))
poialt = POI(Nsig, value=0)
ul_test = UpperLimit(poinull, poialt, calc, CLs=True, qtilde=False)
ul_test.upperlimit();
f, ax = plt.subplots(figsize=(10, 8))
ul_test.plot(ax=ax)
f.savefig("brazilian_plot.png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13915275994865212, 'fraction_numerical': 0.02798459563543004, 'mean_word_length': 6.699604743083004, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '477587', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1530, 'n_tokens_neox': 1436, 'n_words': 309, 'path': 'examples/notebooks/upperlimit_zfit_asy.ipynb', 'repo_name': 'marinang/statrise', 'size': '3895'}