import numpy as np
import matplotlib.pyplot as plt
from iminuit import Minuit

xmin, xmax = 2.6, 3.6

def model(x, ns, nb, mean, sigma, c1):
    
    linear = (1. + c1*x)/((xmax-xmin) + c1*(xmax**2-xmin**2)/2.)
    
    gaussian = 1./(2.*np.pi)**0.5/sigma * \
        np.exp(-0.5*((x-mean)/sigma)**2)

    return ns*gaussian + nb*linear

def fcn(ns, nb, mean, sigma, c1):
    
    L = model(evt, ns, nb, mean, sigma, c1)
    if np.any(L<=0.): return 1E100

    return 2.*(ns+nb)-2.*np.log(L).sum()

values = np.zeros(1000)
errors = np.zeros(1000)

for idx in range(len(values)):
    S = np.random.randn(np.random.poisson(200.))*0.033+3.092
    B = np.random.rand(np.random.poisson(400.))+2.6
    evt = np.hstack([S,B])

    m = Minuit(fcn, ns=200., nb=400., mean=3.092, sigma=0.033, c1=0.)
    m.migrad()
    
    values[idx] = m.values['ns']
    errors[idx] = m.errors['ns']
    
pull = ((values-200.)/errors)[errors!=0.]

print('Pull mean:',pull.mean())
print('Pull width:',pull.std())

fig = plt.figure(figsize=(6,6), dpi=80)
plt.hist(pull, range=(-5.,5.), bins=50)
plt.grid()
plt.show()