First, import the required documents:
In [ ]:
import numpy as np
from becquerel.tools.isotope import Isotope
from becquerel.tools.isotope_qty import IsotopeQuantity, NeutronIrradiation
import datetime
from becquerel import Spectrum
import efficiencies as ef
import pottery_efficiency as pe
from bs4 import BeautifulSoup
import urllib.request
import math
import NAA_Isotopes as na
import spreadsheets as sp
The python file efficiencies is a script in charge of generating the efficiency spectra needed to calculate the efficiency values for the variety of energy peaks used for analysis. The efficiency spectra was created using a variety of sources whose imported spectra.
NAA_Isotopes is a script with a class class containing dictionaries of data for every isotope that is observed in neutron activation analysis (of marine life). This is used to identify the isotope that belongs to an observed energy peak
spreadsheets is a python script containing a class "sheet" is contains the class members: isotope name, energy, branching ratio, counts, and uncertainty that were observed in the spectra. These values were extracted using PeakEasy instead of code generated so the peaks are meant to be accurate and set as default to use as comparison for the values generated by Isotope_Composition_Calculator.py.
Load the spectra for observation as well as the background spectra for background subtraction For now, the spectra must be inputed into the code manually.
In [ ]:
spec_S1 = Spectrum.from_file('/Users/jackiegasca/Documents/spectras/Sample4_24h_C.Spe')
back_spec = Spectrum.from_file('/Users/jackiegasca/Documents/2017.5.1_long_background.Spe')
Convert the bins into energy using the Spectrum attribute from the becquerel package:
In [ ]:
back_spec_ener = back_spec.energies_kev[0:len(back_spec)]
spec_S1_ener_spec = spec_S1.energies_kev[0:len(spec_S1)]
Instead of having the code search for the energy peaks and identifying the number of counts, the isotope list is assigned from the corresponding spreadsheets.py item:
In [ ]:
isotope_list = sp.Sample3_24h.element
The irradiation info is manually inputed since there is no parser for the file containing this information:
In [ ]:
irr_start = '2017-04-27 14:02:00'
irr_stop = '2017-04-27 14:12:00'
flux = 3.1e11
N_0 = 6.02e23
IsotopeActivity() is a function created to convert the counts into activity. (Note: an edit is needed to do background sobtraction).
In [ ]:
def IsotopeActivity():
iso_name = sp.Sample4_24h.element
iso_energy = sp.Sample4_24h.energy
back_cps = []
for i in range(len(iso_energy)):
E = iso_energy[i]
back_peak = (np.abs(back_spec_ener - E)).argmin()
iso_cps = sp.Sample4_24h.counts
iso_cps[:] = [x / spec_S1.livetime for x in iso_cps]
iso_br = sp.Sample4_24h.br
#return(iso_name, iso_energy, iso_cps, iso_br)
isotope_activities = []
for j in range(len(iso_name)):
ef_en = (np.abs(ef.x - iso_energy[j])).argmin()
activ = iso_cps[j] / (ef.fit[ef_en] * iso_br[j])
isotope_activities.append(activ)
return(iso_name, iso_energy, iso_cps, iso_br, isotope_activities)
The function then needs to be called to generate the data needed to proceed with the concentrtion calculations. The values are assigned to their corresponding list as follows:
In [ ]:
lists = IsotopeActivity()
iso_cps = lists[2]
iso_name = lists[0]
print(iso_name)
iso_energy = lists[1]
iso_br = lists[3]
isotope_activities = lists[4]
print(iso_energy)
print(type(isotope_activities[0]))
The function Remover was created to remove any negative activities in the activity list. It must be run a couple of times since it does not remove all the negative activities the first time:
In [ ]:
def Remover():
for i in iso_cps:
if i <= 0:
n = iso_cps.index(i)
iso_name.remove(iso_name[n])
iso_energy.remove(iso_energy[n])
iso_br.remove(iso_br[n])
isotope_activities.remove(isotope_activities[n])
iso_cps.remove(iso_cps[n])
else:
pass
return(iso_name, iso_cps, iso_energy, iso_br, isotope_activities)
lists = Remover()
lists = Remover()
lists = Remover()
The function Concentration() calcualtes the composition of the isotope in the iso_name list using it's activity and branching ratio at the energy peak. The first part of the function converts the isotopes name into a format that is compatible with the becquerel code. The function urlcreator(abb, A_0) takes the name of the element and the preirradiated nuclide number to search through the jaea library to find the cross section of the preirradiated isotope. It returns a link with the cross section information. After the link is accessed by the webscraper, tabledata(bslink) is a function that accesses the link and converts the data from the web into a table that can be parsed. It looks for the thermal neutron cross section and returns that number. The if statement towards the bottom converts the 'k', 'm' and '&' values into their respective scientific notation. The function returns the approproate cross section.
When the cross section is returned, Concentration() uses the tools from becquerel to calculate the concentration, which is returned at the end of the function.
In [ ]:
def Concentration():
conc = []
for i in range(len(iso_name)):
c = iso_name[i].split('_')
abb = c[0]
A = int(c[1])
A_0 = A - 1
iso_2 = '{0}-{1}'.format(abb, A_0)
iso_1 = '{0}-{1}'.format(abb, A)
nuclide = Isotope(iso_1)
def urlcreator(abb, A_0):
A_num = str(A_0)
if len(A_num) == 1:
A_num = '00' + A_num
elif len(A_num) == 2:
A_num = '0' + A_num
else:
A_num = A_num
url = 'http://wwwndc.jaea.go.jp/cgi-bin/Tab80WWW.cgi?/data' \
+ '/JENDL/JENDL-4-prc/intern/' + abb + A_num + '.intern'
html = urllib.request.urlopen(url)
bslink = BeautifulSoup(html, 'lxml')
return(bslink)
bslink = urlcreator(abb, A_0)
def tabledata(bslink):
'''extracts data from the jaea website'''
table = bslink.table
table_rows = table.find_all('tr')
for tr in table_rows:
td = tr.find_all('td')
row = [i.text for i in td]
if len(row) == 7:
if row[0] == 'total ':
x_sec = row[1]
x_sec_s = x_sec.split(' ')
x_val = float(x_sec_s[0])
barn = x_sec_s[1]
if barn[1] == 'k':
x_val = 10**(3) * x_val
return(x_val)
elif barn[1] == 'm':
x_val = 10**(-3) * x_val
return(x_val)
elif barn[1] == '&':
x_val = 10**(-6) * x_val
return(x_val)
else:
x_val = x_val
return(x_val)
else:
pass
else:
pass
return(x_val)
x_val = tabledata(bslink)
print(x_val)
quantity = IsotopeQuantity(nuclide, date=spec_S1.start_time, bq=isotope_activities[i])
irrad_quan = quantity.bq_at(irr_stop)
irrad_act = IsotopeQuantity(nuclide, date=irr_stop, bq=irrad_quan)
ni = NeutronIrradiation(irr_start, irr_stop, n_cm2_s=flux)
init_comp = ni.activate(x_val, initial=Isotope(iso_2), activated=irrad_act)
print(init_comp)
Finally, the function must be called to generate the results, which are then returned onto terminal.
In [ ]:
u = Concentration()
print(u)