Class to read LSC csv files

docs/examples/fit_tritium_release.ipynb

@@ -24,7 +24,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from libra_toolbox.tritium.model import ureg\n",
+    "from libra_toolbox.tritium import ureg\n",
     "from libra_toolbox.tritium.helpers import substract_background_from_measurements\n",
     "\n",
     "background_activity = 0.334 * ureg.Bq\n",
@@ -207,35 +207,34 @@
     "from libra_toolbox.tritium.model import Model\n",
     "import numpy as np\n",
     "\n",
+    "P383_neutron_rate = 4.95e8 / 2 * ureg.neutron * ureg.s**-1\n",
+    "A325_neutron_rate = 2.13e8 / 2 * ureg.neutron * ureg.s**-1\n",
+    "\n",
+    "baby_diameter = 1.77 * ureg.inches - 2 * 0.06 * ureg.inches  # from CAD drawings\n",
+    "baby_radius = 0.5 * baby_diameter\n",
+    "baby_volume = 100 * ureg.mL\n",
+    "baby_cross_section = np.pi * baby_radius**2\n",
+    "baby_height = baby_volume / baby_cross_section\n",
+    "\n",
+    "exposure_time = 12 * ureg.hour\n",
+    "\n",
+    "irradiations = [\n",
+    "    [0 * ureg.hour, 0 + exposure_time],\n",
+    "    [24 * ureg.hour, 24 * ureg.hour + exposure_time],\n",
+    "]\n",
     "\n",
     "def make_model(k_top, k_wall):\n",
-    "    baby_diameter = 1.77 * ureg.inches - 2 * 0.06 * ureg.inches  # from CAD drawings\n",
-    "    baby_radius = 0.5 * baby_diameter\n",
-    "    baby_volume = 100 * ureg.mL\n",
-    "    baby_cross_section = np.pi * baby_radius**2\n",
-    "    baby_height = baby_volume / baby_cross_section\n",
     "\n",
     "    baby_model = Model(\n",
     "        radius=baby_radius,\n",
     "        height=baby_height,\n",
     "        TBR=5.4e-4 * ureg.particle * ureg.neutron**-1,\n",
+    "        k_top=k_top,\n",
+    "        k_wall=k_wall,\n",
+    "        irradiations=irradiations,\n",
+    "        neutron_rate=P383_neutron_rate + A325_neutron_rate,\n",
     "    )\n",
     "\n",
-    "    baby_model.k_top = k_top\n",
-    "    baby_model.k_wall = k_wall\n",
-    "\n",
-    "    exposure_time = 12 * ureg.hour\n",
-    "\n",
-    "    baby_model.irradiations = [\n",
-    "        [0 * ureg.hour, 0 + exposure_time],\n",
-    "        [24 * ureg.hour, 24 * ureg.hour + exposure_time],\n",
-    "    ]\n",
-    "\n",
-    "    P383_neutron_rate = 4.95e8 / 2 * ureg.neutron * ureg.s**-1\n",
-    "    A325_neutron_rate = 2.13e8 / 2 * ureg.neutron * ureg.s**-1\n",
-    "\n",
-    "    baby_model.neutron_rate = P383_neutron_rate + A325_neutron_rate\n",
-    "\n",
     "    return baby_model"
    ]
   },
@@ -582,7 +581,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
@@ -599,7 +598,9 @@
    "source": [
     "from libra_toolbox.tritium.plotting import plot_sample_activity_top\n",
     "\n",
-    "plot_sample_activity_top(model, replacement_times=sampling_times, color=\"black\", label=\"From model\")\n",
+    "plot_sample_activity_top(\n",
+    "    model, replacement_times=sampling_times, color=\"black\", label=\"From model\"\n",
+    ")\n",
     "plot_bars(measurements_after_background_sub, index=sampling_times)\n",
     "\n",
     "plt.legend(reverse=True)\n",

docs/examples/tritium_model.ipynb

@@ -47,7 +47,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from libra_toolbox.tritium import model"
+    "from libra_toolbox.tritium import model, ureg"
    ]
   },
   {
@@ -65,8 +65,7 @@
    "source": [
     "import numpy as np\n",
     "\n",
-    "ureg = model.ureg\n",
-    "\n",
+    "# geometry\n",
     "baby_diameter = 1.77 * ureg.inches - 2 * 0.06 * ureg.inches  # from CAD drawings\n",
     "baby_volume = 0.100 * ureg.L\n",
     "\n",
@@ -75,54 +74,30 @@
     "\n",
     "baby_height = baby_volume / baby_cross_section\n",
     "\n",
-    "my_model = model.Model(\n",
-    "    radius=baby_radius,\n",
-    "    height=baby_height,\n",
-    "    TBR=5.4e-4 * ureg.particle * ureg.neutron**-1,\n",
-    ")\n",
+    "# neutron rate\n",
+    "P383_neutron_rate = 4.95e8 / 2 * ureg.neutron * ureg.s**-1\n",
+    "A325_neutron_rate = 2.13e8 / 2 * ureg.neutron * ureg.s**-1\n",
+    "\n",
+    "neutron_rate = P383_neutron_rate + A325_neutron_rate\n",
+    "\n",
+    "# irradiation schedule\n",
     "\n",
-    "my_model.k_top = 1.6905e-6 * ureg.m * ureg.s**-1\n",
-    "my_model.k_wall = 5.0715e-8 * ureg.m * ureg.s**-1"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "It is possible to add a series of irradiations as a list of tuples:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {},
-   "outputs": [],
-   "source": [
     "exposure_time = 12 * ureg.hour\n",
     "\n",
-    "my_model.irradiations = [\n",
+    "irradiations = [\n",
     "    [0 * ureg.hour, 0 + exposure_time],\n",
     "    [24 * ureg.hour, 24 * ureg.hour + exposure_time],\n",
-    "]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We also provide the neutron rate:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "P383_neutron_rate = 4.95e8 / 2 * ureg.neutron * ureg.s**-1\n",
-    "A325_neutron_rate = 2.13e8 / 2 * ureg.neutron * ureg.s**-1\n",
+    "]\n",
     "\n",
-    "my_model.neutron_rate = P383_neutron_rate + A325_neutron_rate"
+    "my_model = model.Model(\n",
+    "    radius=baby_radius,\n",
+    "    height=baby_height,\n",
+    "    TBR=5.4e-4 * ureg.particle * ureg.neutron**-1,\n",
+    "    neutron_rate=neutron_rate,\n",
+    "    k_top=1.6905e-6 * ureg.m * ureg.s**-1,\n",
+    "    k_wall=5.0715e-8 * ureg.m * ureg.s**-1,\n",
+    "    irradiations=irradiations,\n",
+    ")"
    ]
   },
   {
@@ -134,7 +109,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -168,7 +143,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -219,7 +194,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
@@ -257,7 +232,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
@@ -297,7 +272,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
@@ -328,7 +303,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
     {

libra_toolbox/tritium/__init__.py

@@ -1 +1,7 @@
+import pint
+
+ureg = pint.UnitRegistry()
+ureg.setup_matplotlib()
+ureg.define("neutron = 1 * particle = n")
+
 from . import model

libra_toolbox/tritium/lsc_measurements.py

@@ -0,0 +1,144 @@
+import pandas as pd
+from typing import List
+import pint
+from libra_toolbox.tritium import ureg
+from datetime import datetime, timedelta
+import warnings
+
+
+class LSCFileReader:
+    def __init__(self, file_path, vial_labels=None):
+        self.file_path = file_path
+        self.vial_labels = vial_labels
+        self.data = None
+        self.header_content = None
+
+    def read_file(self):
+        # first read the file without dataframe to find the line starting with S#
+        header_lines = []
+        with open(self.file_path, "r") as file:
+            lines = file.readlines()
+            for i, line in enumerate(lines):
+                if line.startswith("S#"):
+                    start = i
+                    break
+                header_lines.append(line)
+        self.header_content = "".join(header_lines)
+
+        # read the file with dataframe starting from the line with S#
+        self.data = pd.read_csv(self.file_path, skiprows=start)
+
+        # check if last column is all NaN
+        if self.data[self.data.columns[-1]].isnull().all():
+            warnings.warn(
+                "There seem to be an issue with the last column. Is the format of the file correct?"
+            )
+
+    def get_bq1_values(self):
+        return self.data["Bq:1"].tolist()
+
+    def get_bq1_values_with_labels(self):
+        if self.vial_labels is None:
+            raise ValueError("Vial labels must be provided")
+
+        assert len(self.vial_labels) == len(
+            self.get_bq1_values()
+        ), "Vial labels and Bq:1 values are not equal in length, remember to give None as a label for missing vials"
+
+        values = self.get_bq1_values()
+        labelled_values = {label: val for label, val in zip(self.vial_labels, values)}
+
+        return labelled_values
+
+    def get_count_times(self):
+        return self.data["Count Time"].tolist()
+
+    def get_lum(self):
+        return self.data["LUM"].tolist()
+
+
+class LSCSample:
+    def __init__(self, activity: pint.Quantity, name: str):
+        self.activity = activity
+        self.name = name
+        # TODO add other attributes available in LSC file
+        self.background_substracted = False
+
+    def __str__(self):
+        return f"Sample {self.name}"
+
+    def substract_background(self, background_sample: "LSCSample"):
+        if self.background_substracted:
+            raise ValueError("Background already substracted")
+        self.activity -= background_sample.activity
+        self.background_substracted = True
+
+    @staticmethod
+    def from_file(file_reader: LSCFileReader, vial_name):
+        values = file_reader.get_bq1_values_with_labels()
+        if vial_name not in values:
+            raise ValueError(f"Vial {vial_name} not found in the file reader.")
+        activity = values[vial_name] * ureg.Bq
+        return LSCSample(activity, vial_name)
+
+
+class LIBRASample:
+    def __init__(self, samples: List[LSCSample], time: str):
+        self.samples = samples
+        self._time = time
+
+    def get_relative_time(self, start_time: str) -> timedelta:
+        start_time = datetime.strptime(start_time, "%m/%d/%Y %I:%M %p")
+        sample_time = datetime.strptime(self._time, "%m/%d/%Y %I:%M %p")
+        return sample_time - start_time
+
+    def substract_background(self, background_sample: LSCSample):
+        for sample in self.samples:
+            sample.substract_background(background_sample)
+
+    def get_soluble_activity(self):
+        act = 0
+        for sample in self.samples[:2]:
+            act += sample.activity
+
+        return act
+
+    def get_insoluble_activity(self):
+        act = 0
+        for sample in self.samples[2:]:
+            act += sample.activity
+
+        return act
+
+    def get_total_activity(self):
+        return self.get_soluble_activity() + self.get_insoluble_activity()
+
+
+class LIBRARun:
+    def __init__(self, samples: List[LIBRASample], start_time: str):
+        self.samples = samples
+        self.start_time = start_time
+
+    def get_cumulative_activity(self, form: str = "total"):
+        # check that background has been substracted
+        for sample in self.samples:
+            for lsc_sample in sample.samples:
+                if not lsc_sample.background_substracted:
+                    raise ValueError(
+                        "Background must be substracted before calculating cumulative activity"
+                    )
+        cumulative_activity = []
+        for sample in self.samples:
+            if form == "total":
+                cumulative_activity.append(sample.get_total_activity())
+            elif form == "soluble":
+                cumulative_activity.append(sample.get_soluble_activity())
+            elif form == "insoluble":
+                cumulative_activity.append(sample.get_insoluble_activity())
+        cumulative_activity = ureg.Quantity.from_list(cumulative_activity)
+        cumulative_activity = cumulative_activity.cumsum()
+        return cumulative_activity
+
+    @property
+    def relative_times(self):
+        return [sample.get_relative_time(self.start_time) for sample in self.samples]

libra_toolbox/tritium/model.py

@@ -3,9 +3,7 @@
 from scipy.integrate import cumulative_trapezoid
 from scipy.integrate import solve_ivp
 
-ureg = pint.UnitRegistry()
-ureg.setup_matplotlib()
-ureg.define("neutron = 1 * particle = n")
+from libra_toolbox.tritium import ureg
 
 SPECIFIC_ACT = 3.57e14 * ureg.Bq * ureg.g**-1
 MOLAR_MASS = 6.032 / 2 * ureg.g * ureg.mol**-1