Added Experiment class

Author: cdunn314

libra_toolbox/neutron_detection/activation_foils/explicit.py

@@ -1,55 +1,24 @@
 from settings import *
 import foils
-from .calculations import n93_number, delay_time
+from calculations import n93_number, delay_time
 import numpy as np
 import pandas as pd
 import os
+import warnings
 
 
-# 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)
-    
-#     # 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()
-
-#     # 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
+def get_chain(irradiations, decay_constant):
 
-#     return foil
-
-
-
-def get_chain(irradiations, decay_constant=Nb92m_decay_constant):
-    """
-    Returns the value of
+    """ Returns the value of
     (1 - exp(-\lambda * \Delta t_1)) * (1 - exp(-\lambda * \Delta t_2)) * ... * (1 - exp(-\lambda * \Delta t_n))
     where \Delta t_i is the time between the end of the i-th period (rest or irradiation) and the start of the next one
 
     Args:
         irradiations (list): list of dictionaries with keys "t_on" and "t_off" for irradiations
 
     Returns:
-        float or pint.Quantity: the value of the chain
-    """
+        float or pint.Quantity: the value of the chain """
+    
     result = 1
     periods = [{"start": irradiations[0]["t_on"], "end": irradiations[0]["t_off"]}]
     for irr in irradiations[1:]:
@@ -61,6 +30,7 @@ def get_chain(irradiations, decay_constant=Nb92m_decay_constant):
         result = 1 - result * np.exp(-decay_constant * delta_t)
     return result
 
+
 def get_efficiency(energies=None, coeff=None,
                    coeff_energy_bounds=None,
                    coeff_type='total', geometric_eff=1.0):
@@ -72,23 +42,26 @@ def get_efficiency(energies=None, coeff=None,
     if coeff is None:
         # These are the efficiencies reported for activation foil analysis
         # of the BABY 100 mL runs
-        geometric_efficiency = 0.5
+        geometric_eff = 0.5
         intrinsic_efficiency = 0.344917296922981
         total_efficiency = geometric_eff * intrinsic_efficiency
+        warnings.warn('Using NaI efficiency from BABY 100 mL runs, which may not reflect detector efficiency accurately.')
     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)
+        ergs = energies.to(ureg.keV).magnitude
+        erg_bounds = coeff_energy_bounds.to(ureg.keV).magnitude
+        if np.min(ergs) < np.min(erg_bounds) or np.max(ergs) > np.max(erg_bounds):
+            warnings.warn('Efficiency is being extrapolated according to efficiency fit curve bounds.')
+        if coeff_type.lower()=='total':
+            total_efficiency = np.polyval(coeff, ergs)
+        elif coeff_type.lower()=='intrinsic':
+            intrinsic_efficiency = np.polyval(coeff, ergs)
             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):
+def get_neutron_flux(experiment: foils.Experiment, 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 
@@ -108,15 +81,17 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: foils.Foil):
     #     experiment["time_generator_off"], experiment["start_time_counting"]
     # )
     time_between_generator_off_and_start_of_counting = (
-        experiment["start_time_counting"] - experiment["time_generator_off"]
+        experiment.start_time_counting - experiment.time_generator_off
     ).seconds * ureg.second
 
-    if 'total_eff_coeff' in experiment.keys():
-        total_efficiency = get_efficiency(experiment['total_eff_coeff'], foil.photon_energies,
-                                          coeff_energy_bounds=experiment['efficiency_bounds'],
+    if experiment.total_eff_coeff is not None:
+        total_efficiency = get_efficiency(energies=foil.photon_energies,
+                                          coeff=experiment.total_eff_coeff,
+                                          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,
+    elif experiment.intrinsic_eff_coeff is not None:
+        total_efficiency = get_efficiency(energies=foil.photon_energies,
+                                          coeff=experiment.intrinsic_eff_coeff,
                                           coeff_energy_bounds=experiment['efficiency_bounds'],
                                           coeff_type='intrinsic')
     else:
@@ -128,10 +103,11 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: foils.Foil):
     f_spec = total_efficiency * foil.branching_ratio
 
     print('total efficiency', total_efficiency)
-    number_of_decays_measured = experiment["photon_counts"] / f_spec
+    number_of_decays_measured = experiment.photon_counts / f_spec
+
     print('number of decays measured', number_of_decays_measured)
-    print('number', foil.atoms)
-    print('cross section', foil.cross_section)
+    # print('number', foil.atoms)
+    # print('cross section', foil.cross_section)
 
 
     flux = (
@@ -143,14 +119,11 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: foils.Foil):
     f_time = (get_chain(irradiations, foil.decay_constant)
                 * np.exp( -foil.decay_constant
                      * time_between_generator_off_and_start_of_counting) 
-                * (1 - np.exp( -foil.decay_constant * experiment['real_count_time']))
-                * (experiment['live_count_time'] / experiment['real_count_time'])
+                * (1 - np.exp( -foil.decay_constant * experiment.real_count_time))
+                * (experiment.live_count_time / experiment.real_count_time)
                 / 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 - 
@@ -162,23 +135,23 @@ def get_neutron_flux(experiment: dict, irradiations: list, foil: foils.Foil):
                    * foil.thickness)
     ).to('dimensionless')
 
-    print('flux: ', flux.to_reduced_units())
-    print('f_time', f_time)
-    print('f_self', f_self)
-
     flux /= (f_time * f_self)
-    
 
+    # print('flux: ', flux.to_reduced_units())
+    # print('f_time', f_time)
+    # print('f_self', f_self)
 
     # convert n/cm2/s to n/s
-    area_of_sphere = 4 * np.pi * experiment["distance_from_center_of_target_plane"] ** 2
+    area_of_sphere = 4 * np.pi * experiment.distance_from_center_of_target_plane ** 2
 
     flux *= area_of_sphere
 
+    flux = flux.to(1/ureg.s)
+
     return flux
 
 
-def get_neutron_flux_error(experiment: dict):
+def get_neutron_flux_error(experiment: foils.Experiment, foil: foils.Foil):
     """
     Returns the uncertainty of the neutron flux as a pint.Quantity
 
@@ -189,7 +162,7 @@ def get_neutron_flux_error(experiment: dict):
     Returns:
         pint.Quantity: uncertainty of the neutron flux
     """ 
-    error_counts = experiment["photon_counts_uncertainty"] / experiment["photon_counts"]
+    error_counts = experiment.photon_counts_uncertainty / experiment.photon_counts
     error_mass = 0.0001 * ureg.g / experiment["foil_mass"]
     error_geometric_eff = 0.025 / geometric_efficiency
     error_intrinsic_eff = 0.025 / nal_gamma_efficiency

libra_toolbox/neutron_detection/activation_foils/foils.py

@@ -1,6 +1,55 @@
 import numpy as np
+import datetime
+from zoneinfo import ZoneInfo
 from settings import ureg
 
+
+def convert_to_datetime(time:str, tzinfo=ZoneInfo('America/New_York')):
+    if isinstance(time, str):
+        datetime_obj = datetime.datetime.strptime(time, "%m/%d/%Y %H:%M:%S")
+    elif isinstance(time, datetime.datetime):
+        datetime_obj = time
+    else:
+        raise ValueError('Time is neither a string nor a datetime.datetime object.')
+    # Check if timezone info is included
+    if datetime_obj.tzinfo is None:
+        # If not, use New York time as default
+        datetime_obj = datetime_obj.replace(tzinfo=tzinfo)
+    return datetime_obj
+
+
+class Experiment:
+    def __init__(self, time_generator_off,
+                 start_time_counting,
+                 distance_from_center_of_target_plane: ureg.Quantity,
+                 real_count_time: ureg.Quantity,
+                 name=None, generator=None, run=None,
+                 live_count_time=None,
+                 tzinfo=ZoneInfo('America/New_York')):
+        # Ensure times are all datetime objects with timezones
+        self.time_generator_off = convert_to_datetime(time_generator_off)
+        self.start_time_counting = convert_to_datetime(start_time_counting)
+        self.real_count_time = real_count_time
+        self.distance_from_center_of_target_plane = distance_from_center_of_target_plane
+        self.name = name
+        self.generator = generator
+        self.run = run
+
+        if live_count_time:
+            self.live_count_time = live_count_time
+        else:
+            self.live_count_time = real_count_time
+        self.total_eff_coeff = None
+        self.intrinsic_eff_coeff = None
+        self.geometric_efficiency = 1
+        self.efficiency_bounds = np.array([])
+        self.photon_counts = None
+        self.photon_counts_uncertainty = None
+
+    
+
+
+
 class Foil:
     def __init__(self, mass: float, thickness: float, name=''):
         """
@@ -56,7 +105,7 @@ def __init__(self, mass: float, thickness: float, name=''):
         self.branching_ratio = np.array([0.9915])
         # branching ratio used in BABY 100 mL analysis
         # self.branching_ratio = [0.999]
-        self.photon_energies = np.array([934.44])
+        self.photon_energies = np.array([934.44]) * ureg.keV
 
         # Cross section from EAF-2010 at 14 MeV
         self.cross_section = np.array([0.4729]) * ureg.barn
@@ -81,14 +130,17 @@ def __init__(self, mass: float, thickness: float, name=''):
 
         self.reaction = ["Zr90(n,2n)Zr89", 
                          "Zr90(n,2n)Zr89"]
-        self.branching_ratio = np.array([0.9904, 
-                                         0.45])
-        self.photon_energies = np.array([909.15, 
-                                         511]) * ureg.keV
+        self.branching_ratio = np.array([0.45, 
+                                         0.9904])
+        self.photon_energies = np.array([511, 
+                                         909.15]) * ureg.keV
 
         # Cross sections from EAF-2010 at 14 MeV
         self.cross_section = np.array([0.5068,
                                        0.5068]) * ureg.barn
+        # Cross sections from ENDF/B-VIII.0 at 14.1 MeV
+        # self.cross_section = np.array([0.5382,
+        #                                0.5382]) * ureg.barn
 
         self.density = 6.505 * ureg.g / ureg.cm**3 
         self.atomic_mass = 91.222 * ureg.g / ureg.mol
@@ -98,6 +150,27 @@ def __init__(self, mass: float, thickness: float, name=''):
         self.get_atoms()
 
         #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
+        self.mass_attenuation_coefficient = np.array([0.06156,
+                                                      0.08590]) * ureg.cm**2 / ureg.g  # at 1 MeV
+
+
+def convert_dict_to_foils(data:dict):
+    if data["element"] == "Nb":
+        foil = Niobium(mass=data["foil_mass"],
+                       thickness=0.01 * ureg.inch,
+                       name=data["foil_name"])
+    elif data["element"] == "Zr":
+        foil = Zirconium(mass=data["foil_mass"],
+                         thickness=0.005 * ureg.inch,
+                         name=data["foil_name"])
+        
+    # Clean up data dictionary to not include foil class arguments
+    experiment_dict = {}
+    for key in data.keys():
+        if key not in ["element", "foil_mass", "foil_name"]:
+            experiment_dict[key] = data[key]
+    experiment = Experiment(**experiment_dict)
+
+    return experiment, foil
+