In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
In [2]:
Test1, Test2 = np.empty([8, 4096]), np.empty([8, 4096])
for i in range(1, 9):
t1 = "RPi_data/Test_1_p1_g" + str(i) + "_2019-05-28_D3S.csv"
t2 = "RPi_data/Test_2_p1_g" + str(i) + "_2019-05-28_D3S.csv"
Test1[i-1] = np.sum(np.genfromtxt(t1, delimiter= ",").T, axis=1)
Test2[i-1] = np.sum(np.genfromtxt(t2, delimiter= ",").T, axis=1)
In [3]:
for i in range(1, 9):
subplot = int(str(42)+str(i))
#print(subplot)
plt.subplot(subplot)
if i % 2 == 1:
det = i//2 + 1
plt.plot(Test1[det-1])
plt.yscale('log')
plt.title("det " + str(det) + " Thorium")
else:
det = i//2
plt.plot(Test2[det-1], 'g')
plt.yscale('log')
plt.title("det " + str(det) + " Uranium")
plt.subplots_adjust(top=1, bottom=-0.2, left=0, right=1.25, hspace=0.75, wspace=0.35)
plt.show()
In [4]:
for i in range(1, 9):
subplot = int(str(42)+str(i))
#print(subplot)
plt.subplot(subplot)
if i % 2 == 1:
det = i//2 + 1 + 4
plt.plot(Test1[det-1])
plt.yscale('log')
plt.title("det " + str(det) + " Thorium")
else:
det = i//2 + 4
plt.plot(Test2[det-1], 'g')
plt.yscale('log')
plt.title("det " + str(det) + " Uranium")
plt.subplots_adjust(top=1, bottom=-0.2, left=0, right=1.25, hspace=0.75, wspace=0.35)
plt.show()
In [7]:
%%writefile tryPeakfinder.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
'''
This is a helper fuction that looks at each index and checks if it is a peak.
REMOVED:::::Does not look at values under 1/4 of np.average(data):::::::
'''
def checkShape(i, data, r, e):
sweep = [data[i + dx] for dx in range(-r, r+1)]
prev=sweep[r]
if not prev == max(sweep):# or prev < np.average(data)/4:
return False
# if not prev > np.average(sweep) * 1.5:
# return False
e = e * 2
# ^because the code checks r indices to the left and right
for k in range(1, r+1):
if e < 0:
#print(e)
return False
if sweep[r-k] > prev:
e = e - 1
prev = sweep[r-k]
prev=sweep[r]
for k in range(1, r+1):
if e < 0:
return False
if sweep[r+k] > prev:
e = e - 1
prev = sweep[r+k]
return e >= 0
'''
Takes in a summed peak count, a peak range, and an error allowance and returns possible peaks.
Peak range is the number of values the function will look at on either side
Error allowance is the number of values within the peak range that are allowed to not fit a downwards slope
'''
def sweepLeft(data, r, e):
peaks = []
index = r
while index < len(data) - r:
if checkShape(index, data, r, e):
peaks.append(index)
index = index + r - e//2
else:
index += 1
return peaks
for t in range(1, 3):
for d in range(1, 9):
testnum = t
detnum = d
print(t, d)
# testnum = int(input("Enter test number (1, 2): "))
# detnum = int(input("Enter det number (1, 8): "))
gentext = "RPi_data/Test_" + str(testnum) + "_p1_g" + str(detnum) + "_2019-05-28_D3S.csv"
csv = np.genfromtxt(gentext, delimiter= ",").T
summed = np.sum(csv, axis=1)
peakRange = 60 #int(input("Enter a peak range: "))
errAllo = 50 #int(input("Enter an error allowance: "))
ldots = sweepLeft(summed, peakRange, errAllo)
print("returned peaks:", ldots)
print("len peaklist:", len(ldots))
#print(len(ldots))
#print(np.average(summed)/4)
x=np.arange(len(summed))
plt.plot(summed)
#plt.plot(x, np.average(summed)/4 + 0*x)
plt.plot(ldots, summed[ldots], 'ro')
plt.yscale('log')
plt.show()
In [5]:
%run tryPeakfinder.py
In [7]:
# ^the smart way to have done that would've been to create a class for each measurement
# with the peaklist and pyplot graph as class variables to access
# maybe I'll go back and make my work more respectable if I ever need to present this
For thorium, we care about these three peaks:
In English, they look like the highest peak, the next major peak that rises up, and the final major peak.
For uranium, we want the first five "major peaks".
It is easiest to manually pick them out.
In [37]:
t1 = [396, 969, 3786]
t2 = [407, 963, 3758]
t3 = [228, 675, 2782]
t4 = [407, 971, 3750]
t5 = [425, 982, 3770]
t6 = [429, 984, 3750]
t7 = [421, 978, 3792]
t8 = [417, 971, 3727]
u1 = [117, 303, 582, 1005, 1801]
u2 = [143, 327, 596, 1003, 1787]
u3 = [110, 223, 416, 708, 1264]
u4 = [138, 325, 606, 1022, 1815]
u5 = [143, 345, 614, 1034, 1834]
u6 = [157, 345, 614, 1026, 1806]
u7 = [148, 333, 612, 1028, 1797]
u8 = [143, 328, 596, 1011, 1793]
In [38]:
c1 = t1 + u1
c2 = t2 + u2
c3 = t3 + u3
c4 = t4 + u4
c5 = t5 + u5
c6 = t6 + u6
c7 = t7 + u7
c8 = t8 + u8
en = [238.6, 583.1, 2614.7] + [92.6, 185.7, 352, 609.3, 1120.3]
In [39]:
def polyfit(x, b, m, r):
return r * x*x + m*x + b
In [41]:
plt.subplot(421)
plt.plot(c1, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c1, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.grid(True)
plt.plot(polyfit(range(max(c1)), *xpopt), 'g')
plt.show()
plt.subplot(422)
plt.plot(c2, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c2, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c2)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(423)
plt.plot(c3, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c3, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c3)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(424)
plt.plot(c4, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c4, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c4)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(425)
plt.plot(c5, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c5, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c5)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(426)
plt.plot(c6, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c6, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c6)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(427)
plt.plot(c7, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c7, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c7)), *xpopt))
plt.grid(True)
plt.show()
plt.subplot(428)
plt.plot(c8, en, 'r.')
p0 = [.6, 1, 2]
xpopt, xpcov = curve_fit(polyfit, c8, en, p0)
print("parameters:", xpopt)
print("uncertainties:", np.sqrt(np.diag(xpcov)))
print("values: [b,m,r]")
plt.plot(polyfit(range(max(c8)), *xpopt))
plt.grid(True)
plt.show()
Everything looks about right!
In [ ]: