Updated explicit with foil classes

Author: cdunn314

libra_toolbox/neutron_detection/activation_foils/explicit.py

@@ -1,43 +1,44 @@
 from settings import *
-# from .calculations import n93_number, delay_time
+import foils
+from .calculations import n93_number, delay_time
 import numpy as np
 import pandas as pd
 import os
 
 
-def get_foil_data(foil: dict, 
-                  filepath=os.path.join(os.path.dirname(os.path.abspath(__file__)), 'nuclide_data.xlsx'),
-                  xslib='EAF2010'):
+# def get_foil_data(foil: dict, 
+#                   filepath=os.path.join(os.path.dirname(os.path.abspath(__file__)), 'nuclide_data.xlsx'),
+#                   xslib='EAF2010'):
 
-    # Read in nuclide data
-    df = pd.read_excel(filepath, skiprows=1)
+#     # Read in nuclide data
+#     df = pd.read_excel(filepath, skiprows=1)
 
-    # Only get info for one nuclide
-    mask = df['Nuclide'] == foil['nuclide']
-
-    # Calculate number of nuclide atoms in foil
-    num_element = (foil['mass'] 
-                   / (df['Element_Atomic_Mass'][mask].item() * ureg.g / ureg.mol)
-                   * (6.022e23 * ureg.particle / ureg.mol)
-    )
+#     # Only get info for one nuclide
+#     mask = df['Nuclide'] == foil['nuclide']
+
+#     # Calculate number of nuclide atoms in foil
+#     num_element = (foil['mass'] 
+#                    / (df['Element_Atomic_Mass'][mask].item() * ureg.g / ureg.mol)
+#                    * (6.022e23 * ureg.particle / ureg.mol)
+#     )
 
-    foil['number'] = num_element * df['Natural_Abundance'][mask].item()
+#     foil['number'] = num_element * df['Natural_Abundance'][mask].item()
 
-    # Get density and mass attenuation coefficient
-    foil['density'] = df['Density'][mask].item() * ureg.g / ureg.cm**3
-    foil['mass_attenuation_coefficient'] = df['Mass_Attenuation_Coefficient'][mask].item() * (ureg.cm**2/(ureg.g))
+#     # Get density and mass attenuation coefficient
+#     foil['density'] = df['Density'][mask].item() * ureg.g / ureg.cm**3
+#     foil['mass_attenuation_coefficient'] = df['Mass_Attenuation_Coefficient'][mask].item() * (ureg.cm**2/(ureg.g))
 
-    # Get cross section for available reactions
-    for reaction in foil['reactions'].keys():
-        heading = reaction + '_' + xslib + '_xs'
-        if heading in df.keys():
-            foil['reactions'][reaction]['cross_section'] = df[heading][mask].item() * ureg.barn
+#     # Get cross section for available reactions
+#     for reaction in foil['reactions'].keys():
+#         heading = reaction + '_' + xslib + '_xs'
+#         if heading in df.keys():
+#             foil['reactions'][reaction]['cross_section'] = df[heading][mask].item() * ureg.barn
 
-    return foil
+#     return foil
 
 
 
-def get_chain(irradiations, decay_constant):
+def get_chain(irradiations, decay_constant=Nb92m_decay_constant):
     """
     Returns the value of
     (1 - exp(-\lambda * \Delta t_1)) * (1 - exp(-\lambda * \Delta t_2)) * ... * (1 - exp(-\lambda * \Delta t_n))
@@ -60,9 +61,34 @@ def get_chain(irradiations, decay_constant):
         result = 1 - result * np.exp(-decay_constant * delta_t)
     return result
 
-
-def get_neutron_flux(experiment: dict, irradiations: list, foil: dict, 
-                     reaction: str):
+def get_efficiency(energies=None, coeff=None,
+                   coeff_energy_bounds=None,
+                   coeff_type='total', geometric_eff=1.0):
+    """Calculates the total efficiency of a gamma detector
+    based on provided coefficients. These coefficients could
+    be for calculating the total efficiency (coeff_type='total')
+    or for calculating the intrinsic efficiency (coeff_type='intrinsic')"""
+
+    if coeff is None:
+        # These are the efficiencies reported for activation foil analysis
+        # of the BABY 100 mL runs
+        geometric_efficiency = 0.5
+        intrinsic_efficiency = 0.344917296922981
+        total_efficiency = geometric_eff * intrinsic_efficiency
+    else:
+        # Check that efficiency is being interpolated between the bounds of the fit curve
+        if energies.min() < coeff_energy_bounds.min() or energies.max() > coeff_energy_bounds.max():
+            raise Warning('Efficiency is being extrapolated according to efficiency fit curve bounds.')
+        if coeff_type.lower() is 'total':
+            total_efficiency = np.polyval(coeff, energies)
+        elif coeff_type.lower() is 'intrinsic':
+            intrinsic_efficiency = np.polyval(coeff, energies)
+            total_efficiency = geometric_eff * intrinsic_efficiency
+    total_efficiency *= ureg.count / ureg.particle
+    return total_efficiency
+
+
+def get_neutron_flux(experiment: dict, irradiations: list, foil: foils.Foil):
     """calculates the neutron flux during the irradiation
     Based on Equation 1 from:
     Lee, Dongwon, et al. "Determination of the Deuterium-Tritium (D-T) Generator 
@@ -77,9 +103,6 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: dict,
         pint.Quantity: neutron flux
     """
 
-    print('inside get_neutron_flux()')
-    decay_constant = foil['reactions'][reaction]['decay_constant']
-    foil = get_foil_data(foil)
 
     # time_between_generator_off_and_start_of_counting = delay_time(
     #     experiment["time_generator_off"], experiment["start_time_counting"]
@@ -88,44 +111,55 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: dict,
         experiment["start_time_counting"] - experiment["time_generator_off"]
     ).seconds * ureg.second
 
-    overall_efficiency = (
-        (experiment[reaction]['efficiency'] * foil['reactions'][reaction]['branching_ratio'])
-        * ureg.count
-        / ureg.particle
-    )
+    if 'total_eff_coeff' in experiment.keys():
+        total_efficiency = get_efficiency(experiment['total_eff_coeff'], foil.photon_energies,
+                                          coeff_energy_bounds=experiment['efficiency_bounds'],
+                                          coeff_type='total')
+    elif 'intrinsic_eff_coeff' in experiment.keys():
+        total_efficiency = get_efficiency(experiment['intrinsic_eff_coeff'], foil.photon_energies,
+                                          coeff_energy_bounds=experiment['efficiency_bounds'],
+                                          coeff_type='intrinsic')
+    else:
+        total_efficiency = get_efficiency()
+
+    #Spectroscopic Factor to account for the branching ratio and the
+    # total detection efficiency
+
+    f_spec = total_efficiency * foil.branching_ratio
 
-    print('efficiency', overall_efficiency)
-    number_of_decays_measured = experiment[reaction]["photon_counts"] / overall_efficiency
+    print('total efficiency', total_efficiency)
+    number_of_decays_measured = experiment["photon_counts"] / f_spec
     print('number of decays measured', number_of_decays_measured)
-    print('number', foil['number'])
-    print('cross section', foil['reactions'][reaction]['cross_section'])
+    print('number', foil.atoms)
+    print('cross section', foil.cross_section)
 
 
     flux = (
         number_of_decays_measured
-        / foil['number']
-        / foil['reactions'][reaction]['cross_section']
+        / foil.atoms
+        / foil.cross_section
     )
 
-    f_time = (get_chain(irradiations, foil['reactions'][reaction]['decay_constant'])
-                * np.exp( -decay_constant
+    f_time = (get_chain(irradiations, foil.decay_constant)
+                * np.exp( -foil.decay_constant
                      * time_between_generator_off_and_start_of_counting) 
-                * (1 - np.exp( -decay_constant * experiment['real_count_time']))
+                * (1 - np.exp( -foil.decay_constant * experiment['real_count_time']))
                 * (experiment['live_count_time'] / experiment['real_count_time'])
-                / decay_constant
+                / foil.decay_constant
     )
 
     print('att coeff', foil['mass_attenuation_coefficient'])
     print('density', foil['density'])
     print('thickness', foil['thickness'])
 
+    # Correction factor of gamma-ray self-attenuation in the foil
     f_self = ( (1 - 
-                    np.exp(-foil['mass_attenuation_coefficient']
-                           * foil['density']
-                           * foil['thickness']))
-                / (foil['mass_attenuation_coefficient'] 
-                   * foil['density']
-                   * foil['thickness'])
+                    np.exp(-foil.mass_attenuation_coefficient
+                           * foil.density
+                           * foil.thickness))
+                / (foil.mass_attenuation_coefficient
+                   * foil.density
+                   * foil.thickness)
     ).to('dimensionless')
 
     print('flux: ', flux.to_reduced_units())

libra_toolbox/neutron_detection/activation_foils/foils.py

@@ -70,7 +70,7 @@ def __init__(self, mass: float, thickness: float, name=''):
         self.get_atoms()
 
         # Should update mass attenuation coefficient to be photon energy-dependent
-        self.mass_attenuation_coefficient = np.array([0.05866]) * ureg.g / ureg.cm**2   # at 1 MeV
+        self.mass_attenuation_coefficient = np.array([0.06120]) * ureg.cm**2 / ureg.g   # at 1 MeV
 
 class Zirconium(Foil):
     def __init__(self, mass: float, thickness: float, name=''):
@@ -96,8 +96,8 @@ def __init__(self, mass: float, thickness: float, name=''):
         self.abundance = np.array([0.5145,
                                    0.5145])
         self.get_atoms()
-        
-        # Should update mass attenuation coefficient to be photon energy-dependent
-        self.mass_attenuation_coefficient = np.array([0.0581,
-                                                      0.0581]) * ureg.g / ureg.cm**2  # at 1 MeV
+
+        #From NIST Xray Mass Attenuation Coefficients Table 3
+        self.mass_attenuation_coefficient = np.array([0.08590, 
+                                                      0.06156]) * ureg.cm**2 / ureg.g  # at 1 MeV