Merge pull request #37 from aidancrilly/LOS_attenuation

Line of sight attenuation model

example/NeSST Guide.ipynb

@@ -45,7 +45,10 @@
     "* [Additional Materials](#third-bullet)\n",
     "* [Scattering Ion Kinematics](#fourth-bullet)\n",
     "* [DT Fitting Function](#fifth-bullet)\n",
-    "* [Time of Flight synthetic diagnostic](#sixth-bullet)"
+    "* [Time of Flight synthetic diagnostic](#sixth-bullet)\n",
+    "    * [Instrument Response Functions](#seventh-bullet)\n",
+    "    * [Scintillator Sensitivity](#eighth-bullet)\n",
+    "    * [Line of Sight Attenuation](#ninth-bullet)"
    ]
   },
   {
@@ -1053,6 +1056,8 @@
     "\n",
     "Where $E_n$ is the neutron energy, $t_d$ is the detected time, $t_a$ is the arrival time, $S$ is the sensitivity function and $R$ is the instrument response function. For further details see [Hatarik](https://pubs.aip.org/aip/jap/article/118/18/184502/140977/Analysis-of-the-neutron-time-of-flight-spectra) and [Mohamed](https://pubs.aip.org/aip/jap/article/128/21/214501/1078858) papers.\n",
     "\n",
+    "## Instrument Response Function <a class=\"anchor\" id=\"seventh-bullet\"></a>:\n",
+    "\n",
     "Below shows a basic usage of the NeSST time of flight class. If modelling a real detector, basic instrument response terms should be replaced by correct detector specific reponses."
    ]
   },
@@ -1158,7 +1163,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Sensitivity functions:\n",
+    "## Sensitivity functions <a class=\"anchor\" id=\"eighth-bullet\"></a>:\n",
     "\n",
     "Organic scintillators are commonly used for fast neutron detection. The light output function, L(E), of the scintillator can be used to calculate the neutron sensitivity. If we assume knocked on protons dominate the scintillation, the sensitivity is given by:\n",
     "\n",
@@ -1215,6 +1220,32 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "## Line of sight attenuation <a class=\"anchor\" id=\"ninth-bullet\"></a>:\n",
+    "\n",
+    "NeSST also has the capability to create line-of-sight/beamline attenuation models based on nuclear evaluated data. This will require you download the evaluated data for the relevant materials, for example from [ENDF B-VIII.0](https://www.nndc.bnl.gov/endf-b8.0/download.html). Below is an example of how you might at 20m of air to the total sensitivity function:\n",
+    "\n",
+    "```python\n",
+    "Nitrogen = nst.time_of_flight.LOS_material_component(number_fraction=0.78,A=14.0,ENDF_file='n-007_N_014.endf')\n",
+    "Oxygen   = nst.time_of_flight.LOS_material_component(number_fraction=0.21,A=16.0,ENDF_file='n-008_O_016.endf')\n",
+    "Argon    = nst.time_of_flight.LOS_material_component(number_fraction=0.01,A=40.0,ENDF_file='n-018_Ar_040.endf')\n",
+    "# Create list of material components\n",
+    "air_components = [Nitrogen,Oxygen,Argon]\n",
+    "air = nst.time_of_flight.LOS_material(density=1.293,length=20.0,components=air_components)\n",
+    "# Create list of line of sight materials, in this case just air\n",
+    "LOS_materials = [air]\n",
+    "# Compute LOS attenuation model for the various materials in the beam line\n",
+    "LOS_attenuation = nst.time_of_flight.get_LOS_attenuation(LOS_materials)\n",
+    "# Combine this attenuation with the scintillator sensitivity model\n",
+    "total_sensitivity = nst.time_of_flight.combine_detector_sensitivities([NLO_Verbinski_sensitivity,LOS_attenuation])\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# ToF demonstration\n",
+    "\n",
     "We can also demonstrate a forward in time of flight space using NeSST's functionality. Here we consider the case of fitting a DT peak which if affected by fluid velocity effects (giving a Doppler shift and broadening).\n",
     "\n",
     "First, we produce some synthetic data in the nToF space:"

src/NeSST/endf_interface.py

@@ -183,6 +183,16 @@ def retrieve_ENDF_data_from_manifest(manifest):
 
     return return_dict
 
+def retrieve_total_cross_section_from_ENDF_file(filename):
+
+    mat = endf.Material(ENDF_dir+filename)
+
+    total_dataset = mat.section_data[MF_XSECTION,MT_TOTXSECTION]
+    total_table = total_dataset['sigma']
+    total_E, total_sigma = total_table.x, total_table.y
+
+    return total_E, total_sigma 
+
 def convert_to_NeSST_legendre_format(ENDF_al):
     max_Nl = max([len(a) for a in ENDF_al])
     uniform_len_al = [np.pad(a,(0,max_Nl-len(a))) for a in ENDF_al]

src/NeSST/time_of_flight.py

@@ -2,10 +2,55 @@
 from NeSST.collisions import *
 from NeSST.spectral_model import mat_H
 from NeSST.utils import *
+from NeSST.endf_interface import retrieve_total_cross_section_from_ENDF_file
 import numpy as np
 from scipy.integrate import cumtrapz
 from scipy.special import erf
 
+from dataclasses import dataclass,field
+from typing import List
+
+@dataclass
+class LOS_material_component:
+    number_fraction : float
+    A : float
+    ENDF_file : str
+
+@dataclass
+class LOS_material:
+    density : float
+    ntot    : float = field(init=False)
+    mavg    : float = field(init=False)
+    length  : float
+    components : List[LOS_material_component]
+
+    def __post_init__(self):
+        self.mavg = 0.0
+        for comp in self.components:
+            self.mavg += comp.number_fraction*comp.A
+        self.mavg *= sc.atomic_mass
+        self.ntot = self.density/self.mavg
+
+def get_LOS_attenuation(LOS_materials : List[LOS_material]):
+    tau_interp_list = []
+    for LOS_material in LOS_materials:
+        ntot_barn   = 1e-28*LOS_material.ntot
+        L           = LOS_material.length
+        for LOS_component in LOS_material.components:
+            ncomp = ntot_barn*LOS_component.number_fraction
+            E,sigma_tot = retrieve_total_cross_section_from_ENDF_file(LOS_component.ENDF_file)
+            tau_interp  = interpolate_1d(E,L*ncomp*sigma_tot,method='linear')
+            tau_interp_list.append(tau_interp)
+
+    def LOS_attenuation(E):
+        total_tau = tau_interp_list[0](E)
+        for i in range(1,len(tau_interp_list)):
+            total_tau += tau_interp_list[i](E)
+        transmission = np.exp(-total_tau)
+        return transmission
+
+    return LOS_attenuation
+
 def get_power_law_NLO(p,Enorm=E0_DT):
     A = mat_H.sigma(Enorm)*Enorm**p
     def power_law_NLO(E):
@@ -32,6 +77,14 @@ def unity_sensitivity(En):
         return np.ones_like(En)
     return unity_sensitivity
 
+def combine_detector_sensitivities(model_list):
+    def total_sensitivity(E):
+        sensitivity = model_list[0](E)
+        for i in range(1,len(model_list)):
+            sensitivity *= model_list[i](E)
+        return sensitivity
+    return total_sensitivity
+
 def get_transit_time_tophat_IRF(scintillator_thickness):
     def transit_time_tophat_IRF(t_detected,En):
         vn = Ekin_2_beta(En,Mn)*c