Skip to content

Commit

Permalink
Merge pull request #17 from LIBRA-project/lsc-measurements
Browse files Browse the repository at this point in the history
Class to read LSC csv files
  • Loading branch information
RemDelaporteMathurin authored Nov 8, 2024
2 parents 95b0019 + 3b15684 commit 920f3e6
Show file tree
Hide file tree
Showing 12 changed files with 531 additions and 79 deletions.
47 changes: 24 additions & 23 deletions docs/examples/fit_tritium_release.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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"
]
},
Expand Down Expand Up @@ -582,7 +581,7 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 14,
"metadata": {},
"outputs": [
{
Expand All @@ -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",
Expand Down
77 changes: 26 additions & 51 deletions docs/examples/tritium_model.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@
"metadata": {},
"outputs": [],
"source": [
"from libra_toolbox.tritium import model"
"from libra_toolbox.tritium import model, ureg"
]
},
{
Expand All @@ -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",
Expand All @@ -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",
")"
]
},
{
Expand All @@ -134,7 +109,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 3,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -168,7 +143,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 4,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -219,7 +194,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 5,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -257,7 +232,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 6,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -297,7 +272,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 7,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -328,7 +303,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 8,
"metadata": {},
"outputs": [
{
Expand Down
6 changes: 6 additions & 0 deletions libra_toolbox/tritium/__init__.py
Original file line number Diff line number Diff line change
@@ -1 +1,7 @@
import pint

ureg = pint.UnitRegistry()
ureg.setup_matplotlib()
ureg.define("neutron = 1 * particle = n")

from . import model
144 changes: 144 additions & 0 deletions libra_toolbox/tritium/lsc_measurements.py
Original file line number Diff line number Diff line change
@@ -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]
4 changes: 1 addition & 3 deletions libra_toolbox/tritium/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 920f3e6

Please sign in to comment.