<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Pb 212 - 238
Ac 228 - 911
Ti 208 - 2614 keV
Here's that one pdf I found for reference gamma energies
Step1:  Now we find the peak at Pb-212 ~ 238 keV
Step2:  This is good enough for now but we can fix l8tr if needed
Step3:  expected
Step4:  Expected
Step5:  This one's a lot noisier than the others.
Step6:  I didn't like how the channel dot wasn't on top of the peaks so I just changed the x position to be within 10 indices of the mean index.
Step7:  Okay but what about the peak at ~3700?
Step8:  Same here I didn't like how the fit doesnt touch the top of the peaks but irl it makes only aesthetic difference since the offset is so little.
Step9:  Ok now with all the values
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
csv = np.genfromtxt('thorium_test_2019-02-19_D3S.csv', delimiter= ",").T
summed = np.sum(csv[:-1], axis=1) # gets rid of last value
plt.plot(summed)
plt.yscale('log')
plt.show()
Explanation: Pb 212 - 238
Ac 228 - 911
Ti 208 - 2614 keV
Here's that one pdf I found for reference gamma energies: https://www.fmf.uni-lj.si/~jazbinsek/Praktikum5/gammaenergies.pdf
scp dosenet@dosenet.dhcp.lbl.gov:data/thorium_test_2019-02-19_D3S.csv .
First we import data and create summed plot:
End of explanation
Pb_shift = 250
Pb_range = 250
Pb_sample = summed[Pb_shift: Pb_shift + Pb_range]
plt.plot(Pb_sample)
plt.show()
Pb_mu = np.mean(Pb_sample)
Pb_sig = np.std(Pb_sample)
def func(x, a, m, s, c):
    return a * np.exp(-(x - m)**2 / (2 * s**2)) + c                                      
Pb_xdata = range(Pb_range)
trydata = func(Pb_sample, np.max(Pb_sample), Pb_mu, Pb_sig, np.max(Pb_sample) + 50)
p0 = [4000,140,75,1500]
Pb_popt, Pb_pcov = curve_fit(func, Pb_xdata, Pb_sample, p0)
print(Pb_popt)
plt.plot(Pb_xdata, Pb_sample)
plt.plot(Pb_xdata, func(Pb_xdata, *Pb_popt))
plt.plot(int(Pb_popt[1]), func(Pb_xdata, *Pb_popt)[int(Pb_popt[1])], 'ro')
Pb_channel = Pb_shift + int(Pb_popt[1])
plt.show()
Explanation: Now we find the peak at Pb-212 ~ 238 keV:
End of explanation
plt.plot(summed)
plt.plot(Pb_channel, summed[Pb_channel], 'r.')
plt.yscale('log')
print(Pb_channel)
plt.show()
Explanation: This is good enough for now but we can fix l8tr if needed
End of explanation
Ac_shift = 830
Ac_range = 230
Ac_sample = summed[Ac_shift: Ac_shift + Ac_range]
plt.plot(Ac_sample)
plt.show()
Ac_mu = np.mean(Ac_sample)
Ac_sig = np.std(Ac_sample)
Ac_xdata = range(Ac_range)
Ac_trydata = func(Ac_sample, np.max(Ac_sample), Ac_mu, Ac_sig, np.max(Ac_sample) + 50)
Ac_p0 = [700,100,25,300]
Ac_popt, Ac_pcov = curve_fit(func, Ac_xdata, Ac_sample, Ac_p0)
print(Ac_popt)
plt.plot(Ac_xdata, Ac_sample)
plt.plot(Ac_xdata, func(Ac_xdata, *Ac_popt))
plt.plot(int(Ac_popt[1]), func(Ac_xdata, *Ac_popt)[int(Ac_popt[1])], 'ro')
Ac_channel = Ac_shift + int(Ac_popt[1])
plt.show()
plt.plot(summed)
plt.plot(Ac_channel, summed[Ac_channel], 'r.')
plt.yscale('log')
print(Ac_channel)
plt.show()
Explanation: expected: 238 keV, channel: 390
Now for Ac-228:
End of explanation
Ti_shift = 2200
Ti_range = 500
Ti_sample = summed[Ti_shift: Ti_shift + Ti_range]
plt.plot(Ti_sample)
plt.show()
Explanation: Expected: 911 keV, channel: 925
Finally, Ti-208:
End of explanation
Ti_mu = np.mean(Ti_sample)
Ti_sig = np.std(Ti_sample)
Ti_xdata = range(Ti_range)
Ti_trydata = func(Ti_sample, np.max(Ti_sample), Ti_mu, Ti_sig, np.max(Ti_sample) + 50)
Ti_p0 = [80,200,20,50]
Ti_popt, Ti_pcov = curve_fit(func, Ti_xdata, Ti_sample, Ti_p0)
print(Ti_popt)
plt.plot(Ti_xdata, Ti_sample)
plt.plot(Ti_xdata, func(Ti_xdata, *Ti_popt))
plt.plot(int(Ti_popt[1]), func(Ti_xdata, *Ti_popt)[int(Ti_popt[1])], 'ro')
maxish = max(Ti_sample[int(Ti_popt[1]) - 5: int(Ti_popt[1]) + 5])
x_maxish = np.argmax(Ti_sample[int(Ti_popt[1]) - 5: int(Ti_popt[1]) + 5]) + int(Ti_popt[1]) - 5
plt.plot(x_maxish, maxish, 'r^')
Ti_channel = Ti_shift + x_maxish
plt.show()
Explanation: This one's a lot noisier than the others.
End of explanation
plt.plot(summed)
plt.plot(Ti_channel, summed[Ti_channel], 'r.')
plt.yscale('log')
print(Ti_channel)
plt.show()
Explanation: I didn't like how the channel dot wasn't on top of the peaks so I just changed the x position to be within 10 indices of the mean index.
End of explanation
z_shift = 3525
z_range = 300
z_sample = summed[z_shift: z_shift + z_range]
plt.plot(z_sample)
plt.show()
z_mu = np.mean(z_sample)
z_sig = np.std(z_sample)
z_xdata = range(z_range)
z_trydata = func(z_sample, np.max(z_sample), z_mu, z_sig, np.max(z_sample) + 50)
z_p0 = [900,180,20,15]
z_popt, z_pcov = curve_fit(func, z_xdata, z_sample, z_p0)
print(z_popt)
plt.plot(z_xdata, z_sample)
plt.plot(z_xdata, func(z_xdata, *z_popt))
plt.plot(int(z_popt[1]), func(z_xdata, *z_popt)[int(z_popt[1])], 'ro')
z_maxish = max(z_sample[int(z_popt[1]) - 5: int(z_popt[1]) + 5])
x_z_maxish = np.argmax(z_sample[int(z_popt[1]) - 5: int(z_popt[1]) + 5]) + int(z_popt[1]) - 5
plt.plot(x_z_maxish, z_maxish, 'r^')
z_channel = z_shift + x_z_maxish
plt.show()
Explanation: Okay but what about the peak at ~3700?
End of explanation
plt.plot(summed)
plt.plot(z_channel, summed[z_channel], 'r.')
plt.yscale('log')
print(z_channel)
plt.show()
Explanation: Same here I didn't like how the fit doesnt touch the top of the peaks but irl it makes only aesthetic difference since the offset is so little.
End of explanation
plt.plot(summed)
plt.plot(Pb_channel, summed[Pb_channel], 'r.')
plt.plot(Ac_channel, summed[Ac_channel], 'r.')
plt.plot(Ti_channel, summed[Ti_channel], 'r.')
plt.plot(z_channel, summed[z_channel], 'r.')
plt.yscale('log')
print(Pb_channel, Ac_channel, Ti_channel, z_channel)
print(summed[Pb_channel], summed[Ac_channel], summed[Ti_channel], summed[z_channel])
plt.show()
Explanation: Ok now with all the values:
From left to right: Pb-212, Ac-228, Ti-208, and the 'mystery z-peak'.
End of explanation
<END_TASK>