notebook.community

Edit and run



In [173]:

    
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit  # import the curve fitting function
%matplotlib inline



In [174]:

    
CS137Peaks = np.array([416.77])
CS137Energy = np.array([662]) #KeV



In [175]:

    
BA133Peaks = np.array([224.5,189.14,55.23])
BA133Energy = np.array([356,303,81])



In [176]:

    
CO60Peaks = np.array([809.59,714.21])
CO60Energy = np.array([1333,1173])



In [177]:

    
Peaks = np.hstack([CS137Peaks,BA133Peaks,CO60Peaks])



In [178]:

    
Peaks









    Out[178]:





array([ 416.77,  224.5 ,  189.14,   55.23,  809.59,  714.21])



In [179]:

    
Energy = np.hstack([CS137Energy,BA133Energy,CO60Energy])



In [180]:

    
Energy









    Out[180]:





array([ 662,  356,  303,   81, 1333, 1173])



In [181]:

    
UnknownPeaks = np.array([786.45,323.39,123.85])



In [182]:

    
plt.figure(figsize=(10,6));
plt.scatter(Peaks,Energy);

plt.xlabel('MCA Number',fontsize=20);
plt.ylabel('Energy (keV)',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
#plt.savefig('Sample')









    Out[182]:





(array([    0.,   200.,   400.,   600.,   800.,  1000.,  1200.,  1400.]),
 <a list of 8 Text yticklabel objects>)



In [183]:

    
def myfun(N,a,b,c):
    ans = a + b*N + c*N**2  # this is y, "the function to be fit"
    return ans



In [184]:

    
p0 = [-2,1,0]



In [185]:

    
xlots = np.linspace(0,910)  # need lots of data points for smooth curve
yfit = np.zeros((len(Peaks),xlots.size))

plsq, pcov = curve_fit(myfun, Peaks, Energy, p0)  # curve fit returns p and covariance matrix
# these give the parameters and the uncertainties
a = plsq[0]
ea = np.sqrt(pcov[0,0])
b = plsq[1]
eb = np.sqrt(pcov[1,1])
c = plsq[2]
ec = np.sqrt(pcov[2,2])

yfit = myfun(xlots,plsq[0],plsq[1],plsq[2])  # use fit results for a, b, c
    
print('a = %.7f +/- %.7f' % (plsq[0], np.sqrt(pcov[0,0])))
print('b = %.7f +/- %.7f' % (plsq[1], np.sqrt(pcov[1,1])))
print('c = %.7f +/- %.7f' % (plsq[2], np.sqrt(pcov[2,2])))









    



a = -2.3237904 +/- 6.9284232
b = 1.5636096 +/- 0.0420454
c = 0.0001082 +/- 0.0000461

1



In [186]:

    
plt.figure(figsize=(10,6));
plt.scatter(Peaks,Energy);
plt.xlim(0,850)
plt.ylim(0,1600)
plt.xlabel('x (mm)');
plt.ylabel('y (mm)');
plt.plot(xlots,yfit);
plt.legend(['data','Fit'],loc='lower right');
plt.text(20,1400,'a = %.3f +/- %.3f keV' % (plsq[0], np.sqrt(pcov[0,0])),size=17)
plt.text(20,1250,'b = %.3f +/- %.3f keV MCA$^{-1}$' % (plsq[1], np.sqrt(pcov[1,1])),size=17)
plt.text(20,1100,'c = (%.4f +/- %.4f)$\cdot 10^{-3}$ keV MCA$^{-2}$' % (plsq[2]*1e3, np.sqrt(pcov[2,2])*1e3),size=17)
plt.xlabel('MCA Number',fontsize=20);
plt.ylabel('Energy (keV)',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
plt.savefig('Linear')

2



In [187]:

    
E_calc = np.array([myfun(entry,a,b,c) for entry in Peaks])



In [188]:

    
delta_E = np.mean(np.array([abs(E_calc[i] - Energy[i]) for i in np.arange(len(E_calc))]))
delta_E









    Out[188]:





3.6514411763655139

3



In [189]:

    
CS137Energy[0]









    Out[189]:





662



In [190]:

    
E_calc[0]









    Out[190]:





668.131030809398



In [191]:

    
myfun(416.77+22.48/2,a,b,c) - myfun(416.77-22.48/2,a,b,c)









    Out[191]:





37.176875180458865



In [192]:

    
myfun(291,a,b,c)









    Out[192]:





461.84675487365394



In [193]:

    
myfun(130.25,a,b,c)









    Out[193]:





203.17151562419369

6



In [194]:

    
#All in channels
CSfwhm = np.array([22.48])

BAfwhm = np.array([17.03,14,6.2])

COfwhm = np.array([32.1, 30.71])



In [195]:

    
fwhmChannels = np.hstack([CSfwhm,BAfwhm,COfwhm])



In [196]:

    
fwhmEnergy = np.array([abs(myfun(Peaks[i] + fwhmChannels[i]/2,a,b,c)
                           - myfun(Peaks[i] - fwhmChannels[i]/2,a,b,c)) for i in np.arange(len(Peaks))])



In [197]:

    
fwhmEnergy









    Out[197]:





array([ 37.17687518,  27.45540851,  22.46340688,   9.76846181,
        55.81419842,  52.76362264])



In [198]:

    
resolution = np.array([fwhmEnergy[i]/myfun(Peaks[i],a,b,c) for i in np.arange(len(Peaks))])



In [199]:

    
EPeaks = myfun(Peaks,a,b,c)



In [200]:

    
plt.figure(figsize=(10,6));
plt.scatter(EPeaks,resolution);

plt.xlabel('Channel',fontsize=20);
plt.ylabel('Resolution (keV)',fontsize = 20);
plt.xticks(size = 13);
plt.yticks(size = 13);
#plt.savefig('Sample')









    Out[200]:





(array([ 0.02,  0.04,  0.06,  0.08,  0.1 ,  0.12,  0.14]),
 <a list of 7 Text yticklabel objects>)

7



In [201]:

    
def myfun2(E,a,x):
    ans = a*E**x  # this is y, "the function to be fit"
    return ans



In [202]:

    
p02 = [.7,-1/2]



In [203]:

    
EPeaks









    Out[203]:





array([  668.13103081,   354.15848584,   297.28708713,    84.36433394,
        1334.45896011,  1169.60010483])



In [204]:

    
xlots2 = np.linspace(0,1400)  # need lots of data points for smooth curve
yfit2 = np.zeros((len(Peaks),xlots.size))

plsq, pcov = curve_fit(myfun2, EPeaks, resolution, p02)  # curve fit returns p and covariance matrix
# these give the parameters and the uncertainties
aa = plsq[0]
eaa = np.sqrt(pcov[0,0])
x = plsq[1]
ex = np.sqrt(pcov[1,1])


yfit2 = myfun2(xlots2,plsq[0],plsq[1])  # use fit results for a, b, c
    
print('a = %.7f +/- %.7f' % (plsq[0], np.sqrt(pcov[0,0])))
print('x = %.7f +/- %.7f' % (plsq[1], np.sqrt(pcov[1,1])))









    



a = 0.5579101 +/- 0.0663007
x = -0.3511851 +/- 0.0214736



In [205]:

    
plt.figure(figsize=(10,6));
plt.scatter(EPeaks,resolution);
plt.xlim(0,1400)
plt.ylim(0,.2)
plt.xlabel('x (mm)');
plt.ylabel('y (mm)');
plt.plot(xlots2,yfit2);
plt.legend(['fit','data'],loc='best');
plt.text(400,.14,'a = %.3f +/- %.3f' % (plsq[0], np.sqrt(pcov[0,0])),size=17)
plt.text(400,.12,'x = %.3f +/- %.3f' % (plsq[1], np.sqrt(pcov[1,1])),size=17)
plt.xlabel('Energy (keV)',fontsize=20);
plt.ylabel('Resolution',fontsize = 20);

plt.xticks(size = 13);
plt.yticks(size = 13);
plt.savefig('Resolution')

8



In [206]:

    
myfun(786.45,a,b,c)









    Out[206]:





1294.281973868



In [207]:

    
myfun(323.39,a,b,c)









    Out[207]:





514.6447100556918



In [208]:

    
myfun(123.85,a,b,c)









    Out[208]:





192.98849977507311