notebook.community

Edit and run 





In [1]:



    


import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import curve_fit

Models





In [100]:



    


def gauss(x, a, center, sigma, offset):
    return a * np.exp(-(x-center)**2/sigma**2) + offset
def twogauss(x, a1, center1, sigma1, offset1, a2, center2, sigma2, offset2):
    return gauss(x, a1, center1, sigma1, offset1) + gauss(x, a2, center2, sigma2, offset2)
def gausslinear(x, a, center, sigma, offset, slope, intercept):
    return gauss(x, a, center, sigma, offset) + slope * x + intercept
def twogausslinear(x, a, center, sigma, offset, a1, center1, sigma1, offset1, slope, intercept):
    return gausslinear(x, a, center, sigma, offset, slope, intercept) + gauss(x, a1, center1, sigma1, offset1)

Single Gaussian Fit





In [34]:



    


number_of_events = 1000
bin_range = [-10, 10]
bin_width = 0.3



    











In [35]:



    


bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.



    











In [69]:



    


line_strength = [0.4, 0.6]
line = np.random.randn(number_of_events*line_strength[0])
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
photons = np.append(line, background)



    











In [70]:



    


hist, bin_edges = np.histogram(photons, bins = bins)



    











In [71]:



    


ydata = hist
xdata = bin_mids[:-1]
plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)



    


















    
Out[71]:








[<matplotlib.lines.Line2D at 0x108103cd0>]

Fit





In [72]:



    


ydata = hist
xdata = bin_mids[:-1]
fit_guess = [1, 0, 1, 0, 1, 1]
popt, pcov = curve_fit(gausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=fit_guess)



    











In [73]:



    


popt



    


















    
Out[73]:








array([  5.07269116e+01,   2.44597480e-02,  -1.42691635e+00,
         4.14966363e+02,  -8.54953992e-02,  -4.06930553e+02])

















In [77]:



    


plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
#plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, gausslinear(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5]), label='Fit')
plt.legend()



    


















    
Out[77]:








<matplotlib.legend.Legend at 0x10b832190>










    
























In [75]:



    


perr = np.sqrt(np.diag(pcov))



    











In [76]:



    


perr



    


















    
Out[76]:








array([  2.87543297e+00,   5.30219604e-02,   7.19252050e-02,
         2.36600951e+07,   5.81707460e-02,   2.36600951e+07])

















In [101]:



    


number_of_trials = 10000
centers = np.zeros(number_of_trials)
centers_errors = np.zeros(number_of_trials)

def random_gauss_fits(number_of_trials):
    
    line_strength = [0.4, 0.6]

    for i in np.arange(0, number_of_trials):
        line = np.random.randn(number_of_events*line_strength[0])
        background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[1])
        photons = np.append(line, background)
        hist, bin_edges = np.histogram(photons, bins = bins)
        fit_guess = [1, 0, 1, 0, 1, 1]
        popt, pcov = curve_fit(gausslinear, bin_mids[:-1], hist, sigma=np.sqrt(ydata)+1, p0=fit_guess)
        centers[i] = popt[1]
    return centers



    











In [102]:



    


centers = random_gauss_fits(number_of_trials)

Distribution of Fit Centers





In [110]:



    


plt.hist(centers, bins=20, label='$\sigma$ = ' + "%0.3f" % centers.std(), range=[-0.5,0.5])
plt.xlabel("Fit center")
plt.xlim(-0.5, 0.5)
plt.legend()



    


















    
Out[110]:








<matplotlib.legend.Legend at 0x10d7b0d50>

Two Gaussians





In [130]:



    


number_of_events = 1000
line_strength = np.array([0.4, 0.4, 0.2])
line_positions = [0, 10]
bin_range = [-10, 20]
bin_width = 0.3



    











In [131]:



    


bins = np.arange(bin_range[0], bin_range[1], bin_width)
bin_mids = np.arange(bin_range[0], bin_range[1], bin_width) + bin_width/2.



    











In [132]:



    


line1 = np.random.randn(number_of_events*line_strength[0]) + line_positions[0]
line2 = np.random.randn(number_of_events*line_strength[1]) + line_positions[1]
background = np.random.uniform(bin_range[0], bin_range[1], size=number_of_events*line_strength[2])
photons = np.append(line1, np.append(line2, background))



    











In [133]:



    


hist, bin_edges = np.histogram(photons, bins = bins)



    











In [134]:



    


plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.plot(bin_mids[:-1], hist)



    


















    
Out[134]:








[<matplotlib.lines.Line2D at 0x10db8b250>]










    
























In [10]:



    


popt, pcov = curve_fit(twogausslinear, xdata, ydata, sigma=np.sqrt(ydata)+1, p0=[1, 0, 1, 0, 1, 10, 1, 0])



    











In [11]:



    


popt



    


















    
Out[11]:








array([  9.50964076e+01,  -2.44937032e-02,   1.38935303e+00,
         7.99067997e+01,   2.59030560e+01,   1.00315649e+01,
         1.24359943e+00,  -7.99520670e+01])

















In [12]:



    


pcov



    


















    
Out[12]:








array([[  6.11554534e+00,  -1.29398291e-03,  -3.47002223e-02,
          1.78619449e+04,   1.63665243e-03,   7.46230510e-06,
         -8.12015335e-04,  -1.78619393e+04],
       [ -1.29398294e-03,   5.10982021e-04,   2.03100396e-05,
          1.40794303e+03,   4.70732715e-06,   2.57796229e-08,
         -2.19259263e-06,  -1.40794302e+03],
       [ -3.47002223e-02,   2.03100384e-05,   6.02941493e-04,
          3.11669352e+03,  -7.98020166e-05,  -3.51174601e-07,
          4.00132563e-05,  -3.11669382e+03],
       [  1.78626616e+04,   1.40793974e+03,   3.11665515e+03,
         -3.78035305e+13,  -8.80817424e+03,  -1.60122181e+02,
          3.98058125e+02,   3.78035307e+13],
       [  1.63665261e-03,   4.70732710e-06,  -7.98020240e-05,
         -8.80819344e+03,   2.25888267e+00,   6.86695748e-03,
         -4.48027080e-02,   8.80819480e+03],
       [  7.46230819e-06,   2.57796131e-08,  -3.51174754e-07,
         -1.60122287e+02,   6.86695748e-03,   2.23965485e-03,
         -3.95895585e-04,   1.60122294e+02],
       [ -8.12015343e-04,  -2.19259267e-06,   4.00132565e-05,
          3.98060051e+02,  -4.48027080e-02,  -3.95895585e-04,
          2.76672551e-03,  -3.98060704e+02],
       [ -1.78626560e+04,  -1.40793971e+03,  -3.11665542e+03,
          3.78035307e+13,   8.80817560e+03,   1.60122188e+02,
         -3.98058777e+02,  -3.78035310e+13]])

















In [13]:



    


plt.bar(bin_edges[:-1], hist, width = bin_width);
plt.xlim(min(bin_edges), max(bin_edges));
plt.ylim(0)
plt.plot(bin_mids[:-1], hist)
plt.plot(xdata, func(xdata, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], popt[7]), label='Fit')
plt.legend()



    


















    
Out[13]:








<matplotlib.legend.Legend at 0x10780f490>










    
























In [17]:



    


pcov



    


















    
Out[17]:








array([[  6.11554534e+00,  -1.29398291e-03,  -3.47002223e-02,
          1.78619449e+04,   1.63665243e-03,   7.46230510e-06,
         -8.12015335e-04,  -1.78619393e+04],
       [ -1.29398294e-03,   5.10982021e-04,   2.03100396e-05,
          1.40794303e+03,   4.70732715e-06,   2.57796229e-08,
         -2.19259263e-06,  -1.40794302e+03],
       [ -3.47002223e-02,   2.03100384e-05,   6.02941493e-04,
          3.11669352e+03,  -7.98020166e-05,  -3.51174601e-07,
          4.00132563e-05,  -3.11669382e+03],
       [  1.78626616e+04,   1.40793974e+03,   3.11665515e+03,
         -3.78035305e+13,  -8.80817424e+03,  -1.60122181e+02,
          3.98058125e+02,   3.78035307e+13],
       [  1.63665261e-03,   4.70732710e-06,  -7.98020240e-05,
         -8.80819344e+03,   2.25888267e+00,   6.86695748e-03,
         -4.48027080e-02,   8.80819480e+03],
       [  7.46230819e-06,   2.57796131e-08,  -3.51174754e-07,
         -1.60122287e+02,   6.86695748e-03,   2.23965485e-03,
         -3.95895585e-04,   1.60122294e+02],
       [ -8.12015343e-04,  -2.19259267e-06,   4.00132565e-05,
          3.98060051e+02,  -4.48027080e-02,  -3.95895585e-04,
          2.76672551e-03,  -3.98060704e+02],
       [ -1.78626560e+04,  -1.40793971e+03,  -3.11665542e+03,
          3.78035307e+13,   8.80817560e+03,   1.60122188e+02,
         -3.98058777e+02,  -3.78035310e+13]])

















In [14]:



    


perr = np.sqrt(np.diag(pcov))



    


















    






/Users/schriste/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: invalid value encountered in sqrt
  if __name__ == '__main__':

















In [15]:



    


perr



    


















    
Out[15]:








array([ 2.47296287,  0.02260491,  0.02455487,         nan,  1.50295797,
        0.04732499,  0.05259967,         nan])

















In [16]:



    


x = 
plt.plot()



    


















    






  File "<ipython-input-16-ebb8e0862158>", line 1
    x =
        ^
SyntaxError: invalid syntax

















In [ ]:

Content source: ehsteve/ipython-notebooks

Similar notebooks:

notebook.community | gallery | about