# 1. About TDCRPy
`TDCRPy` is a Python code to estimate detection efficiencies of liquid scintillation counters (TDCR or CIEMAT/NIST).
The calculation is based on a photo-physical stochastic model allowing to adress complexe decay schemes and radionuclide mixtures.
The code is developped and maintained by the BIPM (MIT license).
Technical details can be found in
http://dx.doi.org/10.13140/RG.2.2.15682.80321
## 1.1 Installation
`TDCRPy` requires that the following packages are installed in your `Python` environement.
```shell
pip install importlib.resources configparser numpy tqdm setuptools scipy
```
or in `conda` environement:
```shell
conda install importlib.resources configparser numpy tqdm setuptools scipy
```
Then, `TDCRPy` can be installed.
```shell
pip install TDCRPy
```
To obtain the last version.
```shell
pip install TDCRPy --upgrade
```
The module can be imported in your Python code such as.
```python
import tdcrpy
```
## 1.2 Test
To run the unit tests of the package:
```shell
python -m unittest tdcrpy.test.test_tdcrpy
```
or (using the coverage package)
```shell
coverage run -m unittest tdcrpy.test.test_tdcrpy
coverage report -m
```
# 2. Quick start with the TDCRPy code
```python
import tdcrpy
```
## 2.1 Symmetric PMTs
### 2.1.1 Estimation of detection efficiencies given the free parameter
```python
mode = "eff" # ask for efficiency calculation
L = 1.2 # free parameter in keV-1
Rad="Co-60" # radionuclide
pmf_1="1" # relative fraction of the radionulide
N = 1000 # number of Monte Carlo trials
kB =1.0e-5 # Birks constant in cm keV-1
V = 10 # volume of scintillator in mL
```
```python
result = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V, mode)
```
```python
print(f"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}")
print(f"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}")
```
efficiency S = 0.9859 +/- 0.0031
efficiency D = 0.9738 +/- 0.0043
efficiency T = 0.9527 +/- 0.0057
efficiency D (C/N system) = 0.9702 +/- 0.0045
### 2.1.2 Estimation of detection efficiencies given the measured TDCR parameter
```python
TD = 0.977667386529166 # TDCR parameter
```
```python
result = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)
```
global free parameter = 1.1991314352332345 keV-1
```python
print(f"free parameter = {round(result[0],4)} keV-1")
print(f"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}")
print(f"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}")
```
free parameter = 1.1991 keV-1
efficiency S = 0.9858 +/- 0.001
efficiency D = 0.973 +/- 0.0014
efficiency T = 0.9512 +/- 0.0018
efficiency D (C/N system) = 0.9692 +/- 0.0014
## 2.2 Asymmetric PMTs
### 2.2.1 Estimation of detection efficiencies given the free parameters
```python
mode = "eff" # ask for efficiency calculation
L = (1.1, 1.3, 1.2) # free parameter in keV-1
Rad="Co-60" # radionuclide
pmf_1="1" # relative fraction of the radionulide
N = 1000 # number of Monte Carlo trials
kB =1.0e-5 # Birks constant in cm keV-1
V = 10 # volume of scintillator in mL
```
```python
result = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V)
```
```python
print(f"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}")
print(f"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency AB = {round(result[6],4)} +/- {round(result[7],4)}")
print(f"efficiency BC = {round(result[8],4)} +/- {round(result[9],4)}")
print(f"efficiency AC = {round(result[10],4)} +/- {round(result[11],4)}")
print(f"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}")
```
efficiency S = 0.9904 +/- 0.0024
efficiency D = 0.9803 +/- 0.0037
efficiency T = 0.9625 +/- 0.005
efficiency AB = 0.9684 +/- 0.0045
efficiency BC = 0.9696 +/- 0.0044
efficiency AC = 0.9674 +/- 0.0046
efficiency D (C/N system) = 0.9772 +/- 0.0039
### 2.2.2 Estimation of detection efficiencies given the measured TDCR parameters
```python
TD = [0.977667386529166, 0.992232838598821, 0.992343419459002, 0.99275350064608]
```
```python
result = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)
```
global free parameter = 1.2369501044018718 keV-1
free parameters = [1.2369501 1.2369501 1.2369501] keV-1
```python
print(f"Global free parameter = {round(result[0],4)} keV-1")
print(f"free parameter of PMT A = {round(result[1][0],4)} keV-1")
print(f"free parameter of PMT B = {round(result[1][1],4)} keV-1")
print(f"free parameter of PMT C = {round(result[1][2],4)} keV-1")
print(f"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}")
print(f"efficiency AB = {round(result[8],4)} +/- {round(result[9],4)}")
print(f"efficiency BC = {round(result[10],4)} +/- {round(result[11],4)}")
print(f"efficiency AC = {round(result[12],4)} +/- {round(result[13],4)}")
print(f"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}")
```
Global free parameter = 1.237 keV-1
free parameter of PMT A = 1.237 keV-1
free parameter of PMT B = 1.237 keV-1
free parameter of PMT C = 1.237 keV-1
efficiency S = 0.9872 +/- 0.0009
efficiency D = 0.9751 +/- 0.0013
efficiency T = 0.9533 +/- 0.0018
efficiency AB = 0.9606 +/- 0.0016
efficiency BC = 0.9606 +/- 0.0016
efficiency AC = 0.9606 +/- 0.0016
efficiency D (C/N system) = 0.9714 +/- 0.0014
## 2.3 Radionuclide mixture
### 2.3.1 Estimation of detection efficiencies given the free parameter
```python
mode = "eff" # ask for efficiency calculation
mode2 = "sym" # specify that symmetric PMTs is considered
L = 1.2 # free parameter in keV-1
Rad="Co-60, H-3" # radionuclides
pmf_1="0.8, 0.2" # relatives fractions of the radionulides
N = 1000 # number of Monte Carlo trials
kB =1.0e-5 # Birks constant in cm keV-1
V = 10 # volume of scintillator in mL
```
```python
result = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V)
```
```python
print(f"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}")
print(f"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}")
```
efficiency S = 0.9622 +/- 0.0045
efficiency D = 0.9163 +/- 0.007
efficiency T = 0.8482 +/- 0.0096
efficiency D (C/N system) = 0.9037 +/- 0.0074
### 2.3.2 Estimation of detection efficiencies given the measured TDCR parameter
```python
TD = 0.977667386529166 # TDCR parameter
```
```python
result = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)
```
global free parameter = 4.999573818598962 keV-1
```python
print(f"free parameter = {round(result[0],4)} keV-1")
print(f"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}")
print(f"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}")
print(f"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}")
print(f"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}")
```
free parameter = 4.9996 keV-1
efficiency S = 0.9835 +/- 0.001
efficiency D = 0.9677 +/- 0.0015
efficiency T = 0.9397 +/- 0.002
efficiency D (C/N system) = 0.9629 +/- 0.0016
```python
```
# 3. Advanced settings
```python
import tdcrpy as td
```
## 3.1 Read the parameters
```python
print("\nparameters in a list: ", td.TDCR_model_lib.readParameters(disp=True))
```
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
parameters in a list: (1000, 1000, 0.96, 5.2, 11.04, 5, 100.0, 1.5, 2.0, 0.1, 50, 10, 20, False)
## 3.2 Change energy binning the quenching function of electrons
```python
td.TDCR_model_lib.modifynE_electron(100)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifynE_electron(1000)
```
New configuration:
number of integration bins for electrons = 100
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.3 Change energy binning the quenching function of alpha particles
```python
td.TDCR_model_lib.modifynE_alpha(100)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifynE_alpha(1000)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 100
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.4 Change the density (in g/cm3) of the LS source
```python
td.TDCR_model_lib.modifyDensity(1.02)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyDensity(0.96)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 1.02 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.5 Change the mean charge number of the LS source
```python
td.TDCR_model_lib.modifyZ(5.7)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyZ(5.2)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.7
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.6 Change the mean atomic mass number of the LS source
```python
td.TDCR_model_lib.modifyA(12.04)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyA(11.04)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 12.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.7 Change the depht paramerter of the spline interpolation of the quenching function
```python
td.TDCR_model_lib.modifyDepthSpline(7)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyDepthSpline(5)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 7
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.8 Change the energy threshold (in keV) above which the interpolation is applied (for alpha particles)
```python
td.TDCR_model_lib.modifyEinterp_a(200)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyEinterp_a(100)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 200.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.9 Change the energy threshold (in keV) above which the interpolation is applied (for electrons)
```python
td.TDCR_model_lib.modifyEinterp_e(2.0)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyEinterp_e(1.5)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 2.0 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.10 Activate/desactivate the micelles correction
```python
td.TDCR_model_lib.modifyMicCorr(False)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyMicCorr(True)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.11 Change the diameter of reverse micelles (in nm)
```python
td.TDCR_model_lib.modifyDiam_micelle(4.0)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyDiam_micelle(2.0)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 4.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.12 Change the acqueous fraction of the LS source
```python
td.TDCR_model_lib.modifyfAq(0.2)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyfAq(0.1)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.2
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.13 Change the coincidence resolving time of the TDCR counter (in ns)
```python
td.TDCR_model_lib.modifyTau(100)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyTau(50)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 100 ns
extended dead time = 10 µs
measurement time = 20 min
## 3.14 Change the extended dead time of the TDCR counter (in µs)¶
```python
td.TDCR_model_lib.modifyDeadTime(100)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyDeadTime(10)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 100 µs
measurement time = 20 min
## 3.15 Change the measurement time (in min)¶
```python
td.TDCR_model_lib.modifyMeasTime(60)
print("New configuration:")
td.TDCR_model_lib.readParameters(disp=True)
# back to the default value
td.TDCR_model_lib.modifyMeasTime(20)
```
New configuration:
number of integration bins for electrons = 1000
number of integration bins for alpha = 1000
density = 0.96 g/cm3
Z = 5.2
A = 11.04
depth of spline interp. = 5
energy above which interp. in implemented (for alpha) = 100.0 keV
energy above which interp. in implemented (for electron) = 1.5 keV
activation of the micelle correction = False
diameter of micelle = 2.0 nm
acqueous fraction = 0.1
coincidence resolving time = 50 ns
extended dead time = 10 µs
measurement time = 60 min
# 4. Read Beta spectrum
```python
import tdcrpy as td
import matplotlib.pyplot as plt
import numpy as np
```
```python
radionuclide = "Sr-89"
mode = "beta-" # 'beta-' or 'beta+'
level = 'tot' # 0,1,2,3 .... or 'tot'
```
## 4.1 Get information about the BetaShape version
```python
print(td.TDCR_model_lib.readBetaShapeInfo(radionuclide,mode,level))
```
-Beta Spectrum of the main transition from BetaShape 2.4
(05/2024), DDEP 2004 evaluation and Q-value from AME2020-
## 4.2 Read the energy spectrum tabulated from BetaShape
```python
energy, probability = td.TDCR_model_lib.readBetaShape(radionuclide,mode,level)
plt.figure("Energy spectrum")
plt.clf()
plt.plot(energy[:-1], probability)
plt.xlabel(r'$E$ /keV', fontsize=14)
plt.ylabel(r'$p(E)$ /(keV$^{-1}$)', fontsize=14)
```
Text(0, 0.5, '$p(E)$ /(keV$^{-1}$)')
## 4.3 Read the deposited energy spectrum
This spectrum is built by TDCRPy (function `buildBetaSpectra`) to be used for the analytical TDCR model
```python
energy2, probability2 = td.TDCR_model_lib.readBetaSpectra(radionuclide)
plt.figure("Deposied energy spectrum")
plt.clf()
plt.plot(energy[:-1], probability,'-b', label="emitted energy spectrum")
plt.plot(energy2, probability2,'-r', label="deposited energy spectrum")
plt.legend(fontsize=14)
plt.xlabel(r'$E$ /keV', fontsize=14)
plt.ylabel(r'$p(E)$ /(keV$^{-1}$)', fontsize=14)
```
Text(0, 0.5, '$p(E)$ /(keV$^{-1}$)')
```python
```
# 5. Stopping power
```python
import tdcrpy as td
import matplotlib.pyplot as plt
import numpy as np
```
## 5.1 Stopping power at a given energy
```python
density = 0.96 # in g/cm3
energy = 100 # keV
result_alpha = td.TDCR_model_lib.stoppingpowerA(energy,rho=density)
result_electron = td.TDCR_model_lib.stoppingpower(energy*1e3,rho=density)
print(f"dE/dx = {result_alpha:.7g} keV/cm for alpha particles" )
print(f"dE/dx = {result_electron*1e3:.4g} keV/cm for electrons" )
```
dE/dx = 1461120 keV/cm for alpha particles
dE/dx = 3419 keV/cm for electrons
## 5.2 Plot stopping power curves
```python
energy_vec = np.logspace(-2,4,1000) # in keV
w_a, w_e = [], []
for e in energy_vec:
w_a.append(td.TDCR_model_lib.stoppingpowerA(e,rho=density))
w_e.append(1e3*td.TDCR_model_lib.stoppingpower(e*1e3,rho=density))
plt.figure("Stopping power")
plt.clf()
plt.plot(energy_vec,w_a,"-k",label="alpha particles")
plt.plot(energy_vec,w_e,"-r",label="electrons")
plt.xscale('log')
plt.xlabel(r'$E$ /keV', fontsize=14)
plt.ylabel(r'd$E$/d$x$ /(keV/cm)', fontsize=14)
plt.legend()
```
<matplotlib.legend.Legend at 0x25d905a4c90>
# 6. Scintillation quenching model
```python
import tdcrpy as td
import matplotlib.pyplot as plt
import numpy as np
```
## 6.1 Calculate the quenched energy from the Birks model
```python
kB = 0.01 # Birks constant in cm/MeV
ei = 1000 # initial energy in keV
ed = 1000 # deposited energy in keV
nE = 10000 # number of points of the energy linear space
eq_e=td.TDCR_model_lib.E_quench_e(ei*1e3,ed*1e3,kB,nE)*1e-3
print(f"For electrons => Initial energy = {ei} keV, deposited energy = {ed} keV, quenched energy = {eq_e:.5g} keV")
eq_a=td.TDCR_model_lib.E_quench_a(ei,kB*1e-3,nE)
print(f"For alpha particles => Initial energy = {ei} keV, quenched energy = {eq_a:.5g} keV")
```
For electrons => Initial energy = 1000 keV, deposited energy = 1000 keV, quenched energy = 975.24 keV
For alpha particles => Initial energy = 1000 keV, quenched energy = 48.934 keV
## 6.2 Correction of the micelle effect
The deposited energy ratios were evaluation by GEANT4-DNA [1].
[1] Nedjadi, Y., Laedermann, J.-P., Bochud, F., Bailat, C., 2017. On the reverse micelle effect in liquid scintillation counting. Applied Radiation and Isotopes 125, 94–107. https://doi.org/10.1016/j.apradiso.2017.04.020
```python
energy_vec = np.logspace(-2,4,100) # in keV
fAq=0.1
s05, s1, s2, s3, s4, s5, s6, s7, s8, s10 = [], [], [], [], [], [], [], [], [], []
for e in energy_vec:
s05.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=0.5))
s1.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=1.0))
s2.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=2.0))
s3.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=3.0))
s4.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=4.0))
s5.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=5.0))
s6.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=6.0))
s7.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=7.0))
s8.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=8.0))
s10.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=10.0))
plt.figure("Deposited energy ratio")
plt.clf()
plt.plot(energy_vec,s05,label="$\Phi$=0.5 nm")
plt.plot(energy_vec,s1,label="$\Phi$=1.0 nm")
plt.plot(energy_vec,s2,label="$\Phi$=2.0 nm")
plt.plot(energy_vec,s3,label="$\Phi$=3.0 nm")
plt.plot(energy_vec,s4,label="$\Phi$=4.0 nm")
plt.plot(energy_vec,s5,label="$\Phi$=5.0 nm")
plt.plot(energy_vec,s6,label="$\Phi$=6.0 nm")
plt.plot(energy_vec,s7,label="$\Phi$=7.0 nm")
plt.plot(energy_vec,s8,label="$\Phi$=8.0 nm")
plt.plot(energy_vec,s10,label="$\Phi$=10.0 nm")
plt.xscale('log')
plt.xlabel(r'$E$ /keV', fontsize=14)
plt.ylabel(r'$S(E)/E$', fontsize=14)
plt.legend()
```
<matplotlib.legend.Legend at 0x2dab792de10>
## 6.3 Scintillation quenching function from the Birks model (with micelle effect)
```python
energy_vec = np.logspace(-2,4,100) # in keV
fAq = 0.1 # acqueous fraction
diam_micelle = 4.0 # diameter of micelles (in nm)
eq_a, eq_e, eq_e_m = [], [], []
for e in energy_vec:
eq_e.append(td.TDCR_model_lib.E_quench_e(e*1e3,e*1e3,kB,nE)*1e-3)
eq_e_m.append(td.TDCR_model_lib.E_quench_e(e*1e3,e*1e3,kB,nE)*1e-3*td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=diam_micelle))
eq_a.append(td.TDCR_model_lib.E_quench_a(e,kB*1e-3,nE))
plt.figure("Quenched Energy")
plt.clf()
plt.plot(energy_vec,eq_a/energy_vec,"-k",label="alpha particles")
plt.plot(energy_vec,eq_e/energy_vec,"-r",label="electrons")
plt.plot(energy_vec,eq_e_m/energy_vec,"--r",label="electrons (micelle effect)")
plt.xscale('log')
plt.xlabel(r'$E$ /keV', fontsize=14)
plt.ylabel(r'$Q_t(E)/E$', fontsize=14)
plt.legend()
```
<matplotlib.legend.Legend at 0x2dab779f910>
```python
```
```python
# 7. Interaction of ionizing particles
```
```python
# pip install opencv-python
```
```python
import tdcrpy as td
import matplotlib.pyplot as plt
import numpy as np
import cv2
```
## 7.1 Distributions for electrons from 0 keV to 200 keV
```python
A = td.TDCR_model_lib.Matrice10_e_1
C = np.flipud(A[0:])
C = np.clip(C, a_min=7e-6, a_max=1e-1)
C = np.log(C)
C = cv2.GaussianBlur(C, (3, 5), 10)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(0, A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.savefig("electrons_0-200.png")
plt.show()
```
## 7.2 Distributions for electrons from 200 keV to 2 MeV
```python
A = td.TDCR_model_lib.Matrice10_e_2
C = np.flipud(A[0:])
C = np.clip(C, a_min=1e-5, a_max=1e-0)
C = np.log(C)
C = cv2.GaussianBlur(C, (3, 3), 20)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(A[0,0], A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x, rotation=20)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.tight_layout()
plt.savefig("electrons_200-2000.png")
plt.show()
```
## 7.3 Distributions for electrons from 2 MeV to 8 MeV
```python
A = td.TDCR_model_lib.Matrice10_e_3
C = np.flipud(A[0:])
C = np.clip(C, a_min=1e-5, a_max=1e-0)
C = np.log(C)
C = cv2.GaussianBlur(C, (3, 3), 20)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(A[0,0], A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x, rotation=20)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.tight_layout()
plt.savefig("electrons_2000-10000.png")
plt.show()
```
## 7.4 Distributions for photons from 0 keV to 200 keV
```python
A = td.TDCR_model_lib.Matrice10_p_1
C = np.flipud(A[0:])
C = np.clip(C, a_min=1e-7, a_max=1e0)
C = np.log(C)
C = cv2.GaussianBlur(C, (3, 3), 20)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(0, A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.savefig("photons_0-200.png")
plt.show()
x_1 = 1
# Find the column index for the given x_1
col_idx = np.searchsorted(A[0, :], x_1)
# Extract the corresponding y values from the column
y_1 = C[:, col_idx]
x_plot = 0.2*np.arange(0, len(y_1), 1)
print("escape probability = ",np.exp(y_1[-2])/sum(np.exp(y_1)))
# Plot the extracted values
plt.plot(x_plot,np.flipud(y_1))
plt.xlabel(r"$E_1$ /keV", fontsize=14)
plt.ylabel(r"log($P(E_1|E_0)$)", fontsize=14)
plt.title(f"Distribution for $x_1 = {x_1}$ keV")
plt.grid(True)
plt.xlim([0,x_1*1.5])
plt.savefig(f"distribution_x_{x_1}.png")
plt.show()
```
escape probability = 0.5623307322837948
## 7.5 Distributions for photons from 200 keV to 2 MeV
```python
A = td.TDCR_model_lib.Matrice10_p_2
C = np.flipud(A[0:])
C = np.clip(C, a_min=1e-7, a_max=1e0)
C = np.log(C)
C = cv2.GaussianBlur(C, (3, 3), 20)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(A[0,0], A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x, rotation=20)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.tight_layout()
plt.savefig("photons_200-2000.png")
plt.show()
x_1 = 2000
# Find the column index for the given x_1
col_idx = np.searchsorted(A[0, :], x_1)
# Extract the corresponding y values from the column
y_1 = np.log(A[:, col_idx])
x_plot = 2*np.arange(0, len(y_1), 1)
print("escape probability = ",np.exp(y_1[0]))
# Plot the extracted values
plt.plot(x_plot,y_1)
plt.xlabel(r"$E_1$ /keV", fontsize=14)
plt.ylabel(r"log($P(E_1|E_0)$)", fontsize=14)
plt.title(f"Distribution for $x_1 = {x_1}$ keV")
plt.grid(True)
plt.xlim([0,x_1*1.5])
plt.savefig(f"distribution_x_{x_1}.png")
plt.show()
```
escape probability = 1999.9999999999998
## 7.6 Distributions for photons from 2 MeV to 10 MeV
```python
A = td.TDCR_model_lib.Matrice10_p_3
C = np.flipud(A[1:])
C = np.clip(C, a_min=1e-6, a_max=1e-2)
C = np.log(C)
C = cv2.GaussianBlur(C, (5, 5), 20)
extent = [A[0,0], A[0,-1], 0, A[0,-1]]
x = np.arange(A[0,0], A[0,-1], A[0,-1]/10)
y = np.arange(0, A[0,-1], A[0,-1]/10)
plt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')
plt.colorbar()
plt.xticks(x, rotation=20)
plt.yticks(y)
plt.xlabel(r"$E_0$ /keV", fontsize=14)
plt.ylabel(r"$E_1$ /keV", fontsize=14)
plt.tight_layout()
plt.savefig("photons_2000-10000.png")
plt.show()
x_1 = 5000
# Find the column index for the given x_1
col_idx = np.searchsorted(A[0, :], x_1)
# Extract the corresponding y values from the column
y_1 = C[:, col_idx]
x_plot = 10*np.arange(0, len(y_1), 1)
print("escape probability = ",np.exp(y_1[-1])/sum(np.exp(y_1)))
# Plot the extracted values
plt.plot(x_plot,np.flipud(y_1))
plt.xlabel(r"$E_1$ /keV", fontsize=14)
plt.ylabel(r"log($P(E_1|E_0)$)", fontsize=14)
plt.title(f"Distribution for $x_1 = {x_1}$ keV")
plt.grid(True)
plt.xlim([0,x_1*1.5])
plt.savefig(f"distribution_x_{x_1}.png")
plt.show()
```
escape probability = 0.03360101175407434
## 7.7 Sample a deposited energy given an initial energy
```python
ei = 100 # initial energy in keV
v = 10 # volume of the scintillator in mL
ed_g=td.TDCR_model_lib.energie_dep_gamma2(ei,v)
ed_e=td.TDCR_model_lib.energie_dep_beta2(ei,v)
print(f"The gamma ray of initial energy = {ei} keV, has deposited = {ed_g} keV in the scintillant.")
print(f"The electron of initial energy = {ei} keV, has deposited = {ed_e} keV in the scintillant.")
```
The gamma ray of initial energy = 100 keV, has deposited = 0 keV in the scintillant.
The electron of initial energy = 100 keV, has deposited = 100 keV in the scintillant.
## 7.8 Deposited energy of photons as a function of the sintillant volume
```python
V=np.arange(8,21,0.1)
N=100000
E=15 # initial energy in keV
e_vec, ue_vec = [], []
for v in V:
x=[]
for i in range(N):
out=td.TDCR_model_lib.energie_dep_gamma2(15,v)
x.append(out)
e_vec.append(np.mean(x))
ue_vec.append(np.std(x)/np.sqrt(N))
```
```python
plt.figure(r"Ed vs volume")
plt.clf()
plt.errorbar(V, e_vec, yerr=ue_vec, fmt="-k", label=rf'$E_0$ = {E} keV')
plt.legend()
plt.xlabel(r"$V$ /mL", fontsize=14)
plt.ylabel(r"$\bar{y}$ /(keV)", fontsize=14)
```
```python
```
# 8. Efficiency curves
```python
import tdcrpy as td
import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt
```
```python
mode = "eff" # ask for efficiency calculation
Rad="Co-60" # radionuclides
pmf_1="1" # relatives fractions of the radionulides
N = 1000 # number of Monte Carlo trials
kB =1.0e-5 # Birks constant in cm keV-1
V = 10 # volume of scintillator in mL
L=np.logspace(-3,2,num=100) # free parameter in keV-1
# 8.1 Record decay histories in temporary files
td.TDCRPy.TDCRPy(L[0], Rad, pmf_1, N, kB, V, mode, barp=True, record=True)
effS, u_effS, effD, u_effD, effT, u_effT, effD2, u_effD2 = [], [],[], [],[], [], [], []
for l in tqdm(L, desc="free parameters ", unit=" iterations"):
out = td.TDCRPy.TDCRPy(l, Rad, pmf_1, N, kB, V, mode, readRecHist=True)
effS.append(out[2])
u_effS.append(out[3])
effD.append(out[2])
u_effD.append(out[3])
effT.append(out[4])
u_effT.append(out[5])
effD2.append(out[12])
u_effD2.append(out[13])
```
______ ______ ______ _______ ________
|__ __|| ___ \| ___|| ___ | | ____ |
| | | | | || | | | | | | |___| |___ ___
| | | | | || | | |__| | | _____|\ \ | |
| | | |__| || |____| __ \ | | \ \ | |
|_| |_____/ |_____||_| \__\|_| \ \_| |
+++++++++++++++++++++++++++++++++++++++++/ /
________________________________________/ /
|______________________________________________/
version 2.0.2
BIPM 2023 - license MIT
distribution: https://pypi.org/project/TDCRPy
developement: https://github.com/RomainCoulon/TDCRPy
start calculation...
Processing: 100%|█████████████████████████████████████████████████████████████| 1000/1000 [00:24<00:00, 40.79 decays/s]
free parameters : 100%|█████████████████████████████████████████████████████| 100/100 [00:02<00:00, 35.96 iterations/s]
```python
effS=np.asarray(effS)
effT=np.asarray(effT)
effD=np.asarray(effD)
effD2=np.asarray(effD2)
u_effS=np.asarray(u_effS)
u_effT=np.asarray(u_effT)
u_effD=np.asarray(u_effD)
u_effD2=np.asarray(u_effD2)
tdcr=effT/effD
u_tdcr=np.sqrt(u_effD**2*effT**2/effD**4+u_effT**2/effD**2)
plt.figure("efficiency vs free parameter")
plt.clf()
plt.errorbar(L,effD,yerr=u_effD,fmt="-k",label="double coincidences")
plt.errorbar(L,effT,yerr=u_effT,fmt="-r",label="triple coincidences")
plt.errorbar(L,effD2,yerr=u_effD2,fmt="-g",label="double coincidences (CIEMAT/NIST)")
plt.xscale('log')
plt.xlabel(r'$L$ /keV$^{-1}$', fontsize=14)
plt.ylabel(r'$\epsilon$', fontsize=14)
plt.legend()
plt.figure("efficiency vs TDCR")
plt.clf()
plt.errorbar(tdcr,effD,xerr=u_tdcr,yerr=u_effD,fmt="-k")
#plt.xscale('log')
plt.xlabel(r'$R_T/R_D$', fontsize=14)
plt.ylabel(r'$\epsilon_{D}$', fontsize=14)
```
# 9. Efficiency with radionuclide mixtures
```python
import tdcrpy as td
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
```
## 9.1 Fixed parameters
```python
mode = "eff" # ask for efficiency calculation
Rad = "Fe-55, Ni-63" # list of radionuclides
pmf_1="1" # relatives fractions of the radionulides
N = 1000 # number of Monte Carlo trials
kB =1.0e-5 # Birks constant in cm keV-1
V = 10 # volume of scintillator in mL
L=1 # free parameter in keV-1
```
## 9.2 Mixture parameters
```python
C = np.arange(0,1,0.05) # relative fraction of the second nuclide
```
## 9.3 Efficiency calculation
```python
effS, ueffS, effD, ueffD, effT, ueffT, effD2, ueffD2 = [], [], [], [], [], [], [], []
for i in tqdm(C, desc="Processing items"):
pmf_1 = f"{1-i}, {i}"
result = td.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V, mode)
effS.append(result[0]); ueffS.append(result[1])
effD.append(result[2]); ueffD.append(result[3])
effT.append(result[4]); ueffT.append(result[5])
effD2.append(result[12]); ueffD2.append(result[13])
effD = np.asarray(effD); ueffD = np.asarray(ueffD)
effD2 = np.asarray(effD2); ueffD2 = np.asarray(ueffD2)
effT = np.asarray(effT); ueffT = np.asarray(ueffT)
TDCR = effT/effD
```
Processing items: 100%|████████████████████████████████████████████████████████████████| 20/20 [03:52<00:00, 11.63s/it]
## 9.4 Plot the efficiency curves
```python
plt.figure("Stopping power")
plt.clf()
plt.errorbar(C,effD,yerr=2*ueffD,fmt="-k",label=r"$\epsilon_D$")
plt.errorbar(C,TDCR,yerr=2*ueffD,fmt="-r",label=r"$\epsilon_T/\epsilon_D$")
plt.errorbar(C,effD2,yerr=2*ueffD2,fmt="-g",label=r"$\epsilon_{D2}$ (CIEMAT/NIST)")
#plt.xscale('log')
plt.xlabel(f'relative fraction of {Rad.split(",")[1]}', fontsize=14)
plt.ylabel(r' ', fontsize=14)
plt.legend()
plt.savefig("stopping_power_plot.png")
```
# 10. Dynamic efficiency estimation
This notebook presents how to link the `TDCRPy` package with the `radioactivedecay` package [1] to simulate dynamic TDCR measurements.
[1] Alex Malins & Thom Lemoine, radioactivedecay: A Python package for radioactive decay calculations. Journal of Open Source Software, 7 (71), 3318 (2022). https://doi.org/10.21105/joss.03318.
```python
# pip install radioactivedecay
```
```python
import radioactivedecay as rd
import tdcrpy as td
import numpy as np
import matplotlib.pyplot as plt
```
## 10.1 Input parameters
```python
radionuclide = 'Mo-99' # parent nuclide decaying during the measurement
activity_unit = "Bq" # unit of the initial activity
time_unit = "h" # time unit of the decay process
A0 = 1 # initial activity (set to 1 in order to obtain relative activities)
coolingTime = 30.0 # the cooling time
mode = "eff"
L = 1 # free parameter (keV-1)
kB = 1e-5 # Birks constant cm/keV
V = 10 # volume of scintillator (mL)
N = 1000 # number of simulated decays
```
## 10.2 Run radioactivedecay
```python
rad_t0 = rd.Inventory({radionuclide: A0}, activity_unit)
rad_t1 = rad_t0.decay(coolingTime, time_unit)
A_t1 = rad_t1.activities(activity_unit)
As_t1 = sum(A_t1.values())
print(f"Activity at {coolingTime} {time_unit}")
for key, val in A_t1.items(): print(f"\t {key}: {val} {activity_unit}")
print("Total activity = ", As_t1, activity_unit)
print(f"Relative activity at {coolingTime} {time_unit}")
for key, val in A_t1.items(): print(f"\t {key}: {val/As_t1}")
```
Activity at 30.0 h
Mo-99: 0.7295308772422591 Bq
Ru-99: 0.0 Bq
Tc-99: 7.44742326547114e-09 Bq
Tc-99m: 0.6738301487178286 Bq
Total activity = 1.4033610334075108 Bq
Relative activity at 30.0 h
Mo-99: 0.5198454708913215
Ru-99: 0.0
Tc-99: 5.306847695056773e-09
Tc-99m: 0.4801545238018309
## 10.3 Display the decay graph
```python
nuc = rd.Nuclide(radionuclide)
nuc.plot()
```
## 10.4 Display the decay curve
```python
rad_t0 .plot(coolingTime, time_unit, yunits=activity_unit)
```
## 10.5 Efficiency calculation as a function of the time
```python
timeVec = np.arange(0,30,2) # time vector
effS, ueffS, effD, ueffD, effT, ueffT, effD2, ueffD2 = [], [], [], [], [], [], [],[]
for t in timeVec:
rad_t1 = rad_t0.decay(t, time_unit)
A_t1 = rad_t1.activities(activity_unit)
As_t1 = sum(A_t1.values())
rads = ""; pmf1 = ""
for key, val in A_t1.items():
if key != "Pb-208" and key != "Ru-99":
rads += ', '+key
pmf1 += ', '+str(val/As_t1)
rads = rads[2:]
pmf1 = pmf1[2:]
out=td.TDCRPy.TDCRPy(L, rads, pmf1, N, kB, V, mode)
effD.append(out[2]); ueffD.append(out[3]); effT.append(out[4]); ueffT.append(out[5])
effS.append(out[0]); ueffS.append(out[1]); effD2.append(out[12]); ueffD2.append(out[13])
effD = np.asarray(effD); effT = np.asarray(effT); ueffD = np.asarray(ueffD); ueffT = np.asarray(ueffT);
effD2 = np.asarray(effD2); effS = np.asarray(effS); ueffD2 = np.asarray(ueffD2); ueffS = np.asarray(ueffS);
tdcr = effT/effD
```
## 10.6 Plot efficiency curves
```python
# Create the plot
plt.figure(figsize=(10, 6))
plt.errorbar(timeVec, effS, yerr=ueffS, fmt='o', capsize=5, label=r'$\epsilon_S$')
plt.errorbar(timeVec, effD, yerr=ueffD, fmt='o', capsize=5, label=r'$\epsilon_D$')
plt.errorbar(timeVec, tdcr, yerr=ueffT , fmt='o', capsize=5, label=r'$\epsilon_T/\epsilon_D$')
plt.errorbar(timeVec, effD2, yerr=ueffD2, fmt='o', capsize=5, label=r'$\epsilon_D$ (CIEMAT/NIST)')
# Adding titles and labels
#plt.title('Efficiency (effD) as a Function of Time')
plt.xlabel(f'cooling time /{time_unit}',fontsize=16)
plt.ylabel(r'$\epsilon_D$ and $\epsilon_T/\epsilon_D$',fontsize=16)
plt.legend(fontsize=16)
plt.grid(True)
plt.savefig("plotDecay")
# Show the plot
plt.show()
```
```python
```
Raw data
{
"_id": null,
"home_page": "https://pypi.org/project/TDCRPy/",
"name": "TDCRPy",
"maintainer": null,
"docs_url": null,
"requires_python": ">=3.11",
"maintainer_email": null,
"keywords": "Python, TDCR, Monte-Carlo, radionuclide, scintillation, counting",
"author": "RomainCoulon (Romain Coulon)",
"author_email": "<romain.coulon@bipm.org>",
"download_url": "https://files.pythonhosted.org/packages/ec/e1/cf1e37671e2105af63e3c3e83cb1a51a2276b08f1459b938970f717a966a/tdcrpy-2.0.9.tar.gz",
"platform": null,
"description": "# 1. About TDCRPy\n\n`TDCRPy` is a Python code to estimate detection efficiencies of liquid scintillation counters (TDCR or CIEMAT/NIST).\nThe calculation is based on a photo-physical stochastic model allowing to adress complexe decay schemes and radionuclide mixtures.\n\nThe code is developped and maintained by the BIPM (MIT license).\n\nTechnical details can be found in\n\nhttp://dx.doi.org/10.13140/RG.2.2.15682.80321\n\n## 1.1 Installation\n\n`TDCRPy` requires that the following packages are installed in your `Python` environement.\n\n```shell\npip install importlib.resources configparser numpy tqdm setuptools scipy\n```\nor in `conda` environement:\n\n```shell\nconda install importlib.resources configparser numpy tqdm setuptools scipy\n```\n\nThen, `TDCRPy` can be installed.\n\n```shell\npip install TDCRPy\n```\n\nTo obtain the last version.\n\n```shell\npip install TDCRPy --upgrade\n```\n\nThe module can be imported in your Python code such as.\n\n```python\nimport tdcrpy\n```\n\n## 1.2 Test\n\nTo run the unit tests of the package:\n\n```shell\npython -m unittest tdcrpy.test.test_tdcrpy\n```\n\nor (using the coverage package)\n\n```shell\ncoverage run -m unittest tdcrpy.test.test_tdcrpy\ncoverage report -m\n```\n\n# 2. Quick start with the TDCRPy code\n\n```python\nimport tdcrpy\n```\n\n## 2.1 Symmetric PMTs\n\n### 2.1.1 Estimation of detection efficiencies given the free parameter\n\n\n```python\nmode = \"eff\" # ask for efficiency calculation\nL = 1.2 # free parameter in keV-1\nRad=\"Co-60\" # radionuclide\npmf_1=\"1\" # relative fraction of the radionulide\nN = 1000 # number of Monte Carlo trials\nkB =1.0e-5 # Birks constant in cm keV-1\nV = 10 # volume of scintillator in mL\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V, mode)\n```\n\n\n```python\nprint(f\"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}\")\nprint(f\"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}\")\n```\n\n efficiency S = 0.9859 +/- 0.0031\n efficiency D = 0.9738 +/- 0.0043\n efficiency T = 0.9527 +/- 0.0057\n efficiency D (C/N system) = 0.9702 +/- 0.0045\n \n\n### 2.1.2 Estimation of detection efficiencies given the measured TDCR parameter\n\n\n```python\nTD = 0.977667386529166 # TDCR parameter\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)\n```\n\n global free parameter = 1.1991314352332345 keV-1\n \n\n\n```python\nprint(f\"free parameter = {round(result[0],4)} keV-1\")\nprint(f\"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}\")\n```\n\n free parameter = 1.1991 keV-1\n efficiency S = 0.9858 +/- 0.001\n efficiency D = 0.973 +/- 0.0014\n efficiency T = 0.9512 +/- 0.0018\n efficiency D (C/N system) = 0.9692 +/- 0.0014\n \n\n## 2.2 Asymmetric PMTs\n\n### 2.2.1 Estimation of detection efficiencies given the free parameters\n\n\n```python\nmode = \"eff\" # ask for efficiency calculation\nL = (1.1, 1.3, 1.2) # free parameter in keV-1\nRad=\"Co-60\" # radionuclide\npmf_1=\"1\" # relative fraction of the radionulide\nN = 1000 # number of Monte Carlo trials\nkB =1.0e-5 # Birks constant in cm keV-1\nV = 10 # volume of scintillator in mL\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V)\n```\n\n\n```python\nprint(f\"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}\")\nprint(f\"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency AB = {round(result[6],4)} +/- {round(result[7],4)}\")\nprint(f\"efficiency BC = {round(result[8],4)} +/- {round(result[9],4)}\")\nprint(f\"efficiency AC = {round(result[10],4)} +/- {round(result[11],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}\")\n```\n\n efficiency S = 0.9904 +/- 0.0024\n efficiency D = 0.9803 +/- 0.0037\n efficiency T = 0.9625 +/- 0.005\n efficiency AB = 0.9684 +/- 0.0045\n efficiency BC = 0.9696 +/- 0.0044\n efficiency AC = 0.9674 +/- 0.0046\n efficiency D (C/N system) = 0.9772 +/- 0.0039\n \n\n### 2.2.2 Estimation of detection efficiencies given the measured TDCR parameters\n\n\n```python\nTD = [0.977667386529166, 0.992232838598821, 0.992343419459002, 0.99275350064608]\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)\n```\n\n global free parameter = 1.2369501044018718 keV-1\n free parameters = [1.2369501 1.2369501 1.2369501] keV-1\n \n\n\n```python\nprint(f\"Global free parameter = {round(result[0],4)} keV-1\")\nprint(f\"free parameter of PMT A = {round(result[1][0],4)} keV-1\")\nprint(f\"free parameter of PMT B = {round(result[1][1],4)} keV-1\")\nprint(f\"free parameter of PMT C = {round(result[1][2],4)} keV-1\")\nprint(f\"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}\")\nprint(f\"efficiency AB = {round(result[8],4)} +/- {round(result[9],4)}\")\nprint(f\"efficiency BC = {round(result[10],4)} +/- {round(result[11],4)}\")\nprint(f\"efficiency AC = {round(result[12],4)} +/- {round(result[13],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}\")\n```\n\n Global free parameter = 1.237 keV-1\n free parameter of PMT A = 1.237 keV-1\n free parameter of PMT B = 1.237 keV-1\n free parameter of PMT C = 1.237 keV-1\n efficiency S = 0.9872 +/- 0.0009\n efficiency D = 0.9751 +/- 0.0013\n efficiency T = 0.9533 +/- 0.0018\n efficiency AB = 0.9606 +/- 0.0016\n efficiency BC = 0.9606 +/- 0.0016\n efficiency AC = 0.9606 +/- 0.0016\n efficiency D (C/N system) = 0.9714 +/- 0.0014\n \n\n## 2.3 Radionuclide mixture\n\n### 2.3.1 Estimation of detection efficiencies given the free parameter\n\n\n```python\nmode = \"eff\" # ask for efficiency calculation\nmode2 = \"sym\" # specify that symmetric PMTs is considered\nL = 1.2 # free parameter in keV-1\nRad=\"Co-60, H-3\" # radionuclides\npmf_1=\"0.8, 0.2\" # relatives fractions of the radionulides\nN = 1000 # number of Monte Carlo trials\nkB =1.0e-5 # Birks constant in cm keV-1\nV = 10 # volume of scintillator in mL\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V)\n```\n\n\n```python\nprint(f\"efficiency S = {round(result[0],4)} +/- {round(result[1],4)}\")\nprint(f\"efficiency D = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency T = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[12],4)} +/- {round(result[13],4)}\")\n```\n\n efficiency S = 0.9622 +/- 0.0045\n efficiency D = 0.9163 +/- 0.007\n efficiency T = 0.8482 +/- 0.0096\n efficiency D (C/N system) = 0.9037 +/- 0.0074\n \n\n### 2.3.2 Estimation of detection efficiencies given the measured TDCR parameter\n\n\n```python\nTD = 0.977667386529166 # TDCR parameter\n```\n\n\n```python\nresult = tdcrpy.TDCRPy.eff(TD, Rad, pmf_1, kB, V)\n```\n\n global free parameter = 4.999573818598962 keV-1\n \n\n\n```python\nprint(f\"free parameter = {round(result[0],4)} keV-1\")\nprint(f\"efficiency S = {round(result[2],4)} +/- {round(result[3],4)}\")\nprint(f\"efficiency D = {round(result[4],4)} +/- {round(result[5],4)}\")\nprint(f\"efficiency T = {round(result[6],4)} +/- {round(result[7],4)}\")\nprint(f\"efficiency D (C/N system) = {round(result[14],4)} +/- {round(result[15],4)}\")\n```\n\n free parameter = 4.9996 keV-1\n efficiency S = 0.9835 +/- 0.001\n efficiency D = 0.9677 +/- 0.0015\n efficiency T = 0.9397 +/- 0.002\n efficiency D (C/N system) = 0.9629 +/- 0.0016\n \n\n\n```python\n\n```\n\n# 3. Advanced settings\n\n```python\nimport tdcrpy as td\n```\n\n## 3.1 Read the parameters\n\n\n```python\nprint(\"\\nparameters in a list: \", td.TDCR_model_lib.readParameters(disp=True))\n```\n\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n parameters in a list: (1000, 1000, 0.96, 5.2, 11.04, 5, 100.0, 1.5, 2.0, 0.1, 50, 10, 20, False)\n \n\n## 3.2 Change energy binning the quenching function of electrons \n\n\n```python\ntd.TDCR_model_lib.modifynE_electron(100)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifynE_electron(1000)\n```\n\n New configuration:\n number of integration bins for electrons = 100\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.3 Change energy binning the quenching function of alpha particles\n\n\n```python\ntd.TDCR_model_lib.modifynE_alpha(100)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifynE_alpha(1000)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 100\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.4 Change the density (in g/cm3) of the LS source\n\n\n```python\ntd.TDCR_model_lib.modifyDensity(1.02)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyDensity(0.96)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 1.02 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.5 Change the mean charge number of the LS source\n\n\n```python\ntd.TDCR_model_lib.modifyZ(5.7)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyZ(5.2)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.7\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.6 Change the mean atomic mass number of the LS source\n\n\n```python\ntd.TDCR_model_lib.modifyA(12.04)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyA(11.04)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 12.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.7 Change the depht paramerter of the spline interpolation of the quenching function\n\n\n```python\ntd.TDCR_model_lib.modifyDepthSpline(7)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyDepthSpline(5)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 7\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.8 Change the energy threshold (in keV) above which the interpolation is applied (for alpha particles)\n\n\n```python\ntd.TDCR_model_lib.modifyEinterp_a(200)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyEinterp_a(100)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 200.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.9 Change the energy threshold (in keV) above which the interpolation is applied (for electrons)\n\n\n```python\ntd.TDCR_model_lib.modifyEinterp_e(2.0)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyEinterp_e(1.5)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 2.0 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.10 Activate/desactivate the micelles correction\n\n\n```python\ntd.TDCR_model_lib.modifyMicCorr(False)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyMicCorr(True)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.11 Change the diameter of reverse micelles (in nm)\n\n\n```python\ntd.TDCR_model_lib.modifyDiam_micelle(4.0)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyDiam_micelle(2.0)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 4.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.12 Change the acqueous fraction of the LS source\n\n\n```python\ntd.TDCR_model_lib.modifyfAq(0.2)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyfAq(0.1)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.2\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.13 Change the coincidence resolving time of the TDCR counter (in ns)\n\n\n```python\ntd.TDCR_model_lib.modifyTau(100)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyTau(50)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 100 ns\n extended dead time = 10 \u00b5s\n measurement time = 20 min\n \n\n## 3.14 Change the extended dead time of the TDCR counter (in \u00b5s)\u00b6\n\n\n```python\ntd.TDCR_model_lib.modifyDeadTime(100)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyDeadTime(10)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 100 \u00b5s\n measurement time = 20 min\n \n\n## 3.15 Change the measurement time (in min)\u00b6\n\n\n```python\ntd.TDCR_model_lib.modifyMeasTime(60)\nprint(\"New configuration:\")\ntd.TDCR_model_lib.readParameters(disp=True)\n# back to the default value\ntd.TDCR_model_lib.modifyMeasTime(20)\n```\n\n New configuration:\n number of integration bins for electrons = 1000\n number of integration bins for alpha = 1000\n density = 0.96 g/cm3\n Z = 5.2\n A = 11.04\n depth of spline interp. = 5\n energy above which interp. in implemented (for alpha) = 100.0 keV\n energy above which interp. in implemented (for electron) = 1.5 keV\n activation of the micelle correction = False\n diameter of micelle = 2.0 nm\n acqueous fraction = 0.1\n coincidence resolving time = 50 ns\n extended dead time = 10 \u00b5s\n measurement time = 60 min\n \n# 4. Read Beta spectrum\n\n```python\nimport tdcrpy as td\nimport matplotlib.pyplot as plt\nimport numpy as np\n```\n\n\n```python\nradionuclide = \"Sr-89\"\nmode = \"beta-\" # 'beta-' or 'beta+'\nlevel = 'tot' # 0,1,2,3 .... or 'tot'\n```\n\n## 4.1 Get information about the BetaShape version\n\n\n```python\nprint(td.TDCR_model_lib.readBetaShapeInfo(radionuclide,mode,level))\n```\n\n -Beta Spectrum of the main transition from BetaShape 2.4\n (05/2024), DDEP 2004 evaluation and Q-value from AME2020-\n \n\n## 4.2 Read the energy spectrum tabulated from BetaShape\n\n\n```python\nenergy, probability = td.TDCR_model_lib.readBetaShape(radionuclide,mode,level)\n\nplt.figure(\"Energy spectrum\")\nplt.clf()\nplt.plot(energy[:-1], probability)\nplt.xlabel(r'$E$ /keV', fontsize=14)\nplt.ylabel(r'$p(E)$ /(keV$^{-1}$)', fontsize=14)\n```\n\n\n\n\n Text(0, 0.5, '$p(E)$ /(keV$^{-1}$)')\n\n\n\n\n \n\n \n\n\n## 4.3 Read the deposited energy spectrum\n\nThis spectrum is built by TDCRPy (function `buildBetaSpectra`) to be used for the analytical TDCR model\n\n\n```python\nenergy2, probability2 = td.TDCR_model_lib.readBetaSpectra(radionuclide)\n\nplt.figure(\"Deposied energy spectrum\")\nplt.clf()\nplt.plot(energy[:-1], probability,'-b', label=\"emitted energy spectrum\")\nplt.plot(energy2, probability2,'-r', label=\"deposited energy spectrum\")\nplt.legend(fontsize=14)\nplt.xlabel(r'$E$ /keV', fontsize=14)\nplt.ylabel(r'$p(E)$ /(keV$^{-1}$)', fontsize=14)\n```\n\n\n\n\n Text(0, 0.5, '$p(E)$ /(keV$^{-1}$)')\n\n\n\n\n \n\n \n\n\n\n```python\n\n```\n\n# 5. Stopping power\n\n```python\nimport tdcrpy as td\nimport matplotlib.pyplot as plt\nimport numpy as np\n```\n\n## 5.1 Stopping power at a given energy\n\n\n```python\ndensity = 0.96 # in g/cm3\nenergy = 100 # keV\nresult_alpha = td.TDCR_model_lib.stoppingpowerA(energy,rho=density)\nresult_electron = td.TDCR_model_lib.stoppingpower(energy*1e3,rho=density)\nprint(f\"dE/dx = {result_alpha:.7g} keV/cm for alpha particles\" )\nprint(f\"dE/dx = {result_electron*1e3:.4g} keV/cm for electrons\" )\n```\n\n dE/dx = 1461120 keV/cm for alpha particles\n dE/dx = 3419 keV/cm for electrons\n \n\n## 5.2 Plot stopping power curves\n\n\n```python\nenergy_vec = np.logspace(-2,4,1000) # in keV\nw_a, w_e = [], []\nfor e in energy_vec:\n w_a.append(td.TDCR_model_lib.stoppingpowerA(e,rho=density))\n w_e.append(1e3*td.TDCR_model_lib.stoppingpower(e*1e3,rho=density))\nplt.figure(\"Stopping power\")\nplt.clf()\nplt.plot(energy_vec,w_a,\"-k\",label=\"alpha particles\")\nplt.plot(energy_vec,w_e,\"-r\",label=\"electrons\")\nplt.xscale('log')\nplt.xlabel(r'$E$ /keV', fontsize=14)\nplt.ylabel(r'd$E$/d$x$ /(keV/cm)', fontsize=14)\nplt.legend()\n```\n\n\n\n\n <matplotlib.legend.Legend at 0x25d905a4c90>\n\n\n\n\n \n\n \n# 6. Scintillation quenching model\n\n```python\nimport tdcrpy as td\nimport matplotlib.pyplot as plt\nimport numpy as np\n```\n\n## 6.1 Calculate the quenched energy from the Birks model\n\n\n```python\nkB = 0.01 # Birks constant in cm/MeV\nei = 1000 # initial energy in keV\ned = 1000 # deposited energy in keV\nnE = 10000 # number of points of the energy linear space\neq_e=td.TDCR_model_lib.E_quench_e(ei*1e3,ed*1e3,kB,nE)*1e-3\nprint(f\"For electrons => Initial energy = {ei} keV, deposited energy = {ed} keV, quenched energy = {eq_e:.5g} keV\")\neq_a=td.TDCR_model_lib.E_quench_a(ei,kB*1e-3,nE)\nprint(f\"For alpha particles => Initial energy = {ei} keV, quenched energy = {eq_a:.5g} keV\")\n```\n\n For electrons => Initial energy = 1000 keV, deposited energy = 1000 keV, quenched energy = 975.24 keV\n For alpha particles => Initial energy = 1000 keV, quenched energy = 48.934 keV\n \n\n## 6.2 Correction of the micelle effect\n\nThe deposited energy ratios were evaluation by GEANT4-DNA [1].\n\n[1] Nedjadi, Y., Laedermann, J.-P., Bochud, F., Bailat, C., 2017. On the reverse micelle effect in liquid scintillation counting. Applied Radiation and Isotopes 125, 94\u2013107. https://doi.org/10.1016/j.apradiso.2017.04.020 \n\n\n```python\nenergy_vec = np.logspace(-2,4,100) # in keV\nfAq=0.1\ns05, s1, s2, s3, s4, s5, s6, s7, s8, s10 = [], [], [], [], [], [], [], [], [], []\nfor e in energy_vec:\n s05.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=0.5))\n s1.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=1.0))\n s2.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=2.0))\n s3.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=3.0))\n s4.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=4.0))\n s5.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=5.0))\n s6.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=6.0))\n s7.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=7.0))\n s8.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=8.0))\n s10.append(td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=10.0))\n\nplt.figure(\"Deposited energy ratio\")\nplt.clf()\nplt.plot(energy_vec,s05,label=\"$\\Phi$=0.5 nm\")\nplt.plot(energy_vec,s1,label=\"$\\Phi$=1.0 nm\")\nplt.plot(energy_vec,s2,label=\"$\\Phi$=2.0 nm\")\nplt.plot(energy_vec,s3,label=\"$\\Phi$=3.0 nm\")\nplt.plot(energy_vec,s4,label=\"$\\Phi$=4.0 nm\")\nplt.plot(energy_vec,s5,label=\"$\\Phi$=5.0 nm\")\nplt.plot(energy_vec,s6,label=\"$\\Phi$=6.0 nm\")\nplt.plot(energy_vec,s7,label=\"$\\Phi$=7.0 nm\")\nplt.plot(energy_vec,s8,label=\"$\\Phi$=8.0 nm\")\nplt.plot(energy_vec,s10,label=\"$\\Phi$=10.0 nm\")\nplt.xscale('log')\nplt.xlabel(r'$E$ /keV', fontsize=14)\nplt.ylabel(r'$S(E)/E$', fontsize=14)\nplt.legend()\n```\n\n\n\n\n <matplotlib.legend.Legend at 0x2dab792de10>\n\n\n\n\n \n\n \n\n\n## 6.3 Scintillation quenching function from the Birks model (with micelle effect)\n\n\n```python\nenergy_vec = np.logspace(-2,4,100) # in keV\nfAq = 0.1 # acqueous fraction\ndiam_micelle = 4.0 # diameter of micelles (in nm) \n\neq_a, eq_e, eq_e_m = [], [], []\nfor e in energy_vec:\n eq_e.append(td.TDCR_model_lib.E_quench_e(e*1e3,e*1e3,kB,nE)*1e-3)\n eq_e_m.append(td.TDCR_model_lib.E_quench_e(e*1e3,e*1e3,kB,nE)*1e-3*td.TDCR_model_lib.micelleLoss(e,fAq=fAq,diam_micelle=diam_micelle))\n eq_a.append(td.TDCR_model_lib.E_quench_a(e,kB*1e-3,nE))\n\nplt.figure(\"Quenched Energy\")\nplt.clf()\nplt.plot(energy_vec,eq_a/energy_vec,\"-k\",label=\"alpha particles\")\nplt.plot(energy_vec,eq_e/energy_vec,\"-r\",label=\"electrons\")\nplt.plot(energy_vec,eq_e_m/energy_vec,\"--r\",label=\"electrons (micelle effect)\")\nplt.xscale('log')\nplt.xlabel(r'$E$ /keV', fontsize=14)\nplt.ylabel(r'$Q_t(E)/E$', fontsize=14)\nplt.legend()\n```\n\n\n\n\n <matplotlib.legend.Legend at 0x2dab779f910>\n\n\n\n\n \n\n \n\n\n\n```python\n\n```\n\n```python\n# 7. Interaction of ionizing particles\n```\n\n```python\n# pip install opencv-python\n```\n\n\n```python\nimport tdcrpy as td\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport cv2\n```\n\n## 7.1 Distributions for electrons from 0 keV to 200 keV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_e_1\nC = np.flipud(A[0:])\nC = np.clip(C, a_min=7e-6, a_max=1e-1)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (3, 5), 10)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(0, A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.savefig(\"electrons_0-200.png\")\nplt.show()\n```\n\n\n \n\n \n\n\n## 7.2 Distributions for electrons from 200 keV to 2 MeV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_e_2\nC = np.flipud(A[0:])\nC = np.clip(C, a_min=1e-5, a_max=1e-0)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (3, 3), 20)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(A[0,0], A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x, rotation=20)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.tight_layout()\nplt.savefig(\"electrons_200-2000.png\")\nplt.show()\n```\n\n\n \n\n \n\n\n## 7.3 Distributions for electrons from 2 MeV to 8 MeV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_e_3\nC = np.flipud(A[0:])\nC = np.clip(C, a_min=1e-5, a_max=1e-0)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (3, 3), 20)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(A[0,0], A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x, rotation=20)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.tight_layout()\nplt.savefig(\"electrons_2000-10000.png\")\nplt.show()\n```\n\n\n \n\n \n\n\n## 7.4 Distributions for photons from 0 keV to 200 keV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_p_1\nC = np.flipud(A[0:])\nC = np.clip(C, a_min=1e-7, a_max=1e0)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (3, 3), 20)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(0, A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.savefig(\"photons_0-200.png\")\nplt.show()\n\nx_1 = 1\n\n# Find the column index for the given x_1\ncol_idx = np.searchsorted(A[0, :], x_1)\n\n# Extract the corresponding y values from the column\ny_1 = C[:, col_idx]\nx_plot = 0.2*np.arange(0, len(y_1), 1)\n\nprint(\"escape probability = \",np.exp(y_1[-2])/sum(np.exp(y_1)))\n# Plot the extracted values\nplt.plot(x_plot,np.flipud(y_1))\nplt.xlabel(r\"$E_1$ /keV\", fontsize=14)\nplt.ylabel(r\"log($P(E_1|E_0)$)\", fontsize=14)\nplt.title(f\"Distribution for $x_1 = {x_1}$ keV\")\nplt.grid(True)\nplt.xlim([0,x_1*1.5])\nplt.savefig(f\"distribution_x_{x_1}.png\")\nplt.show()\n```\n\n\n \n\n \n escape probability = 0.5623307322837948\n \n\n \n\n\n## 7.5 Distributions for photons from 200 keV to 2 MeV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_p_2\nC = np.flipud(A[0:])\nC = np.clip(C, a_min=1e-7, a_max=1e0)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (3, 3), 20)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(A[0,0], A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x, rotation=20)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.tight_layout()\nplt.savefig(\"photons_200-2000.png\")\nplt.show()\n\nx_1 = 2000\n\n# Find the column index for the given x_1\ncol_idx = np.searchsorted(A[0, :], x_1)\n\n# Extract the corresponding y values from the column\ny_1 = np.log(A[:, col_idx])\nx_plot = 2*np.arange(0, len(y_1), 1)\n\nprint(\"escape probability = \",np.exp(y_1[0]))\n# Plot the extracted values\nplt.plot(x_plot,y_1)\nplt.xlabel(r\"$E_1$ /keV\", fontsize=14)\nplt.ylabel(r\"log($P(E_1|E_0)$)\", fontsize=14)\nplt.title(f\"Distribution for $x_1 = {x_1}$ keV\")\nplt.grid(True)\nplt.xlim([0,x_1*1.5])\nplt.savefig(f\"distribution_x_{x_1}.png\")\nplt.show()\n```\n\n\n \n\n \n escape probability = 1999.9999999999998\n \n\n \n\n\n## 7.6 Distributions for photons from 2 MeV to 10 MeV\n\n\n```python\nA = td.TDCR_model_lib.Matrice10_p_3\nC = np.flipud(A[1:])\nC = np.clip(C, a_min=1e-6, a_max=1e-2)\nC = np.log(C)\nC = cv2.GaussianBlur(C, (5, 5), 20)\nextent = [A[0,0], A[0,-1], 0, A[0,-1]]\nx = np.arange(A[0,0], A[0,-1], A[0,-1]/10)\ny = np.arange(0, A[0,-1], A[0,-1]/10)\n\nplt.imshow(C, extent=extent, cmap='Greys', interpolation='nearest')\nplt.colorbar()\nplt.xticks(x, rotation=20)\nplt.yticks(y)\nplt.xlabel(r\"$E_0$ /keV\", fontsize=14)\nplt.ylabel(r\"$E_1$ /keV\", fontsize=14)\nplt.tight_layout()\nplt.savefig(\"photons_2000-10000.png\")\nplt.show()\n\nx_1 = 5000\n\n# Find the column index for the given x_1\ncol_idx = np.searchsorted(A[0, :], x_1)\n\n# Extract the corresponding y values from the column\ny_1 = C[:, col_idx]\nx_plot = 10*np.arange(0, len(y_1), 1)\n\nprint(\"escape probability = \",np.exp(y_1[-1])/sum(np.exp(y_1)))\n# Plot the extracted values\nplt.plot(x_plot,np.flipud(y_1))\nplt.xlabel(r\"$E_1$ /keV\", fontsize=14)\nplt.ylabel(r\"log($P(E_1|E_0)$)\", fontsize=14)\nplt.title(f\"Distribution for $x_1 = {x_1}$ keV\")\nplt.grid(True)\nplt.xlim([0,x_1*1.5])\nplt.savefig(f\"distribution_x_{x_1}.png\")\nplt.show()\n```\n\n\n \n\n \n escape probability = 0.03360101175407434\n \n\n \n\n\n## 7.7 Sample a deposited energy given an initial energy\n\n\n```python\nei = 100 # initial energy in keV\nv = 10 # volume of the scintillator in mL\ned_g=td.TDCR_model_lib.energie_dep_gamma2(ei,v)\ned_e=td.TDCR_model_lib.energie_dep_beta2(ei,v)\nprint(f\"The gamma ray of initial energy = {ei} keV, has deposited = {ed_g} keV in the scintillant.\")\nprint(f\"The electron of initial energy = {ei} keV, has deposited = {ed_e} keV in the scintillant.\")\n```\n\n The gamma ray of initial energy = 100 keV, has deposited = 0 keV in the scintillant.\n The electron of initial energy = 100 keV, has deposited = 100 keV in the scintillant.\n \n\n## 7.8 Deposited energy of photons as a function of the sintillant volume\n\n\n```python\nV=np.arange(8,21,0.1)\nN=100000\nE=15 # initial energy in keV\ne_vec, ue_vec = [], []\nfor v in V:\n x=[]\n for i in range(N):\n out=td.TDCR_model_lib.energie_dep_gamma2(15,v)\n x.append(out)\n e_vec.append(np.mean(x))\n ue_vec.append(np.std(x)/np.sqrt(N))\n\n```\n\n\n```python\nplt.figure(r\"Ed vs volume\")\nplt.clf()\nplt.errorbar(V, e_vec, yerr=ue_vec, fmt=\"-k\", label=rf'$E_0$ = {E} keV')\nplt.legend()\nplt.xlabel(r\"$V$ /mL\", fontsize=14)\nplt.ylabel(r\"$\\bar{y}$ /(keV)\", fontsize=14)\n```\n\n \n\n \n\n```python\n\n```\n\n# 8. Efficiency curves\n\n```python\nimport tdcrpy as td\nimport numpy as np\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\n```\n\n\n```python\nmode = \"eff\" # ask for efficiency calculation\nRad=\"Co-60\" # radionuclides\npmf_1=\"1\" # relatives fractions of the radionulides\nN = 1000 # number of Monte Carlo trials\nkB =1.0e-5 # Birks constant in cm keV-1\nV = 10 # volume of scintillator in mL\nL=np.logspace(-3,2,num=100) # free parameter in keV-1\n\n# 8.1 Record decay histories in temporary files\ntd.TDCRPy.TDCRPy(L[0], Rad, pmf_1, N, kB, V, mode, barp=True, record=True)\n\neffS, u_effS, effD, u_effD, effT, u_effT, effD2, u_effD2 = [], [],[], [],[], [], [], []\nfor l in tqdm(L, desc=\"free parameters \", unit=\" iterations\"):\n out = td.TDCRPy.TDCRPy(l, Rad, pmf_1, N, kB, V, mode, readRecHist=True)\n effS.append(out[2])\n u_effS.append(out[3])\n effD.append(out[2])\n u_effD.append(out[3])\n effT.append(out[4])\n u_effT.append(out[5])\n effD2.append(out[12])\n u_effD2.append(out[13])\n```\n\n \n ______ ______ ______ _______ ________\n |__ __|| ___ \\| ___|| ___ | | ____ |\n | | | | | || | | | | | | |___| |___ ___\n | | | | | || | | |__| | | _____|\\ \\ | |\n | | | |__| || |____| __ \\ | | \\ \\ | |\n |_| |_____/ |_____||_| \\__\\|_| \\ \\_| |\n +++++++++++++++++++++++++++++++++++++++++/ /\n ________________________________________/ /\n |______________________________________________/ \n \n \n version 2.0.2\n BIPM 2023 - license MIT \n distribution: https://pypi.org/project/TDCRPy \n developement: https://github.com/RomainCoulon/TDCRPy \n \n start calculation...\n \n\n Processing: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000/1000 [00:24<00:00, 40.79 decays/s]\n free parameters : 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 100/100 [00:02<00:00, 35.96 iterations/s]\n \n\n\n```python\neffS=np.asarray(effS)\neffT=np.asarray(effT)\neffD=np.asarray(effD)\neffD2=np.asarray(effD2)\nu_effS=np.asarray(u_effS)\nu_effT=np.asarray(u_effT)\nu_effD=np.asarray(u_effD)\nu_effD2=np.asarray(u_effD2)\n\ntdcr=effT/effD\nu_tdcr=np.sqrt(u_effD**2*effT**2/effD**4+u_effT**2/effD**2)\n\nplt.figure(\"efficiency vs free parameter\")\nplt.clf()\nplt.errorbar(L,effD,yerr=u_effD,fmt=\"-k\",label=\"double coincidences\")\nplt.errorbar(L,effT,yerr=u_effT,fmt=\"-r\",label=\"triple coincidences\")\nplt.errorbar(L,effD2,yerr=u_effD2,fmt=\"-g\",label=\"double coincidences (CIEMAT/NIST)\")\nplt.xscale('log')\nplt.xlabel(r'$L$ /keV$^{-1}$', fontsize=14)\nplt.ylabel(r'$\\epsilon$', fontsize=14)\nplt.legend()\n\nplt.figure(\"efficiency vs TDCR\")\nplt.clf()\nplt.errorbar(tdcr,effD,xerr=u_tdcr,yerr=u_effD,fmt=\"-k\")\n#plt.xscale('log')\nplt.xlabel(r'$R_T/R_D$', fontsize=14)\nplt.ylabel(r'$\\epsilon_{D}$', fontsize=14)\n```\n\n \n\n \n \n\n \n\n# 9. Efficiency with radionuclide mixtures\n\n```python\nimport tdcrpy as td\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom tqdm import tqdm\n```\n\n## 9.1 Fixed parameters\n\n```python\nmode = \"eff\" # ask for efficiency calculation\nRad = \"Fe-55, Ni-63\" # list of radionuclides\npmf_1=\"1\" # relatives fractions of the radionulides\nN = 1000 # number of Monte Carlo trials\nkB =1.0e-5 # Birks constant in cm keV-1\nV = 10 # volume of scintillator in mL\nL=1 # free parameter in keV-1\n```\n\n## 9.2 Mixture parameters\n\n```python\nC = np.arange(0,1,0.05) # relative fraction of the second nuclide\n```\n\n## 9.3 Efficiency calculation\n\n```python\neffS, ueffS, effD, ueffD, effT, ueffT, effD2, ueffD2 = [], [], [], [], [], [], [], []\nfor i in tqdm(C, desc=\"Processing items\"):\n pmf_1 = f\"{1-i}, {i}\"\n result = td.TDCRPy.TDCRPy(L, Rad, pmf_1, N, kB, V, mode)\n effS.append(result[0]); ueffS.append(result[1])\n effD.append(result[2]); ueffD.append(result[3])\n effT.append(result[4]); ueffT.append(result[5])\n effD2.append(result[12]); ueffD2.append(result[13])\neffD = np.asarray(effD); ueffD = np.asarray(ueffD)\neffD2 = np.asarray(effD2); ueffD2 = np.asarray(ueffD2)\neffT = np.asarray(effT); ueffT = np.asarray(ueffT)\nTDCR = effT/effD\n```\n\n Processing items: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 20/20 [03:52<00:00, 11.63s/it]\n \n\n## 9.4 Plot the efficiency curves\n\n\n```python\nplt.figure(\"Stopping power\")\nplt.clf()\nplt.errorbar(C,effD,yerr=2*ueffD,fmt=\"-k\",label=r\"$\\epsilon_D$\")\nplt.errorbar(C,TDCR,yerr=2*ueffD,fmt=\"-r\",label=r\"$\\epsilon_T/\\epsilon_D$\")\nplt.errorbar(C,effD2,yerr=2*ueffD2,fmt=\"-g\",label=r\"$\\epsilon_{D2}$ (CIEMAT/NIST)\")\n#plt.xscale('log')\nplt.xlabel(f'relative fraction of {Rad.split(\",\")[1]}', fontsize=14)\nplt.ylabel(r' ', fontsize=14)\nplt.legend()\nplt.savefig(\"stopping_power_plot.png\")\n```\n\n\n \n# 10. Dynamic efficiency estimation\n\nThis notebook presents how to link the `TDCRPy` package with the `radioactivedecay` package [1] to simulate dynamic TDCR measurements.\n\n[1] Alex Malins & Thom Lemoine, radioactivedecay: A Python package for radioactive decay calculations. Journal of Open Source Software, 7 (71), 3318 (2022). https://doi.org/10.21105/joss.03318.\n\n\n```python\n# pip install radioactivedecay\n```\n\n\n```python\nimport radioactivedecay as rd\nimport tdcrpy as td\nimport numpy as np\nimport matplotlib.pyplot as plt\n```\n\n## 10.1 Input parameters\n\n\n```python\nradionuclide = 'Mo-99' # parent nuclide decaying during the measurement\nactivity_unit = \"Bq\" # unit of the initial activity\ntime_unit = \"h\" # time unit of the decay process\nA0 = 1 # initial activity (set to 1 in order to obtain relative activities)\ncoolingTime = 30.0 # the cooling time\n\nmode = \"eff\"\nL = 1 # free parameter (keV-1)\nkB = 1e-5 # Birks constant cm/keV\nV = 10 # volume of scintillator (mL)\nN = 1000 # number of simulated decays\n```\n\n## 10.2 Run radioactivedecay\n\n\n```python\nrad_t0 = rd.Inventory({radionuclide: A0}, activity_unit)\nrad_t1 = rad_t0.decay(coolingTime, time_unit)\nA_t1 = rad_t1.activities(activity_unit)\nAs_t1 = sum(A_t1.values())\nprint(f\"Activity at {coolingTime} {time_unit}\") \nfor key, val in A_t1.items(): print(f\"\\t {key}: {val} {activity_unit}\")\nprint(\"Total activity = \", As_t1, activity_unit)\nprint(f\"Relative activity at {coolingTime} {time_unit}\")\nfor key, val in A_t1.items(): print(f\"\\t {key}: {val/As_t1}\")\n```\n\n Activity at 30.0 h\n \t Mo-99: 0.7295308772422591 Bq\n \t Ru-99: 0.0 Bq\n \t Tc-99: 7.44742326547114e-09 Bq\n \t Tc-99m: 0.6738301487178286 Bq\n Total activity = 1.4033610334075108 Bq\n Relative activity at 30.0 h\n \t Mo-99: 0.5198454708913215\n \t Ru-99: 0.0\n \t Tc-99: 5.306847695056773e-09\n \t Tc-99m: 0.4801545238018309\n \n\n## 10.3 Display the decay graph\n\n\n```python\nnuc = rd.Nuclide(radionuclide)\nnuc.plot()\n```\n \n\n \n## 10.4 Display the decay curve\n\n\n```python\nrad_t0 .plot(coolingTime, time_unit, yunits=activity_unit)\n```\n \n\n \n\n\n## 10.5 Efficiency calculation as a function of the time\n\n\n```python\ntimeVec = np.arange(0,30,2) # time vector\neffS, ueffS, effD, ueffD, effT, ueffT, effD2, ueffD2 = [], [], [], [], [], [], [],[]\nfor t in timeVec:\n rad_t1 = rad_t0.decay(t, time_unit)\n A_t1 = rad_t1.activities(activity_unit)\n As_t1 = sum(A_t1.values())\n \n rads = \"\"; pmf1 = \"\"\n for key, val in A_t1.items():\n if key != \"Pb-208\" and key != \"Ru-99\":\n rads += ', '+key\n pmf1 += ', '+str(val/As_t1)\n rads = rads[2:]\n pmf1 = pmf1[2:]\n \n out=td.TDCRPy.TDCRPy(L, rads, pmf1, N, kB, V, mode)\n effD.append(out[2]); ueffD.append(out[3]); effT.append(out[4]); ueffT.append(out[5])\n effS.append(out[0]); ueffS.append(out[1]); effD2.append(out[12]); ueffD2.append(out[13])\n\neffD = np.asarray(effD); effT = np.asarray(effT); ueffD = np.asarray(ueffD); ueffT = np.asarray(ueffT);\neffD2 = np.asarray(effD2); effS = np.asarray(effS); ueffD2 = np.asarray(ueffD2); ueffS = np.asarray(ueffS);\ntdcr = effT/effD\n```\n \n\n## 10.6 Plot efficiency curves\n\n\n```python\n# Create the plot\nplt.figure(figsize=(10, 6))\nplt.errorbar(timeVec, effS, yerr=ueffS, fmt='o', capsize=5, label=r'$\\epsilon_S$')\nplt.errorbar(timeVec, effD, yerr=ueffD, fmt='o', capsize=5, label=r'$\\epsilon_D$')\nplt.errorbar(timeVec, tdcr, yerr=ueffT , fmt='o', capsize=5, label=r'$\\epsilon_T/\\epsilon_D$')\nplt.errorbar(timeVec, effD2, yerr=ueffD2, fmt='o', capsize=5, label=r'$\\epsilon_D$ (CIEMAT/NIST)')\n\n# Adding titles and labels\n#plt.title('Efficiency (effD) as a Function of Time')\nplt.xlabel(f'cooling time /{time_unit}',fontsize=16)\nplt.ylabel(r'$\\epsilon_D$ and $\\epsilon_T/\\epsilon_D$',fontsize=16)\nplt.legend(fontsize=16)\nplt.grid(True)\n\nplt.savefig(\"plotDecay\")\n# Show the plot\nplt.show()\n```\n\n\n \n\n \n\n\n\n```python\n\n```\n\n\n\n",
"bugtrack_url": null,
"license": null,
"summary": "TDCR model",
"version": "2.0.9",
"project_urls": {
"Documentation": "https://github.com/RomainCoulon/TDCRPy/",
"Homepage": "https://pypi.org/project/TDCRPy/"
},
"split_keywords": [
"python",
" tdcr",
" monte-carlo",
" radionuclide",
" scintillation",
" counting"
],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "73aae63ffb5e4c1378f57e5c6629059eacab760ecaa3c3219a86400326d9f7d4",
"md5": "ad158fe67c047f4cf2f91f992e75ce62",
"sha256": "047a827c28be6501533b4d8a7c9e65334ff343179c3ee5e1106d6c9c2b8bf789"
},
"downloads": -1,
"filename": "TDCRPy-2.0.9-py3-none-any.whl",
"has_sig": false,
"md5_digest": "ad158fe67c047f4cf2f91f992e75ce62",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": ">=3.11",
"size": 22982340,
"upload_time": "2024-09-12T11:33:37",
"upload_time_iso_8601": "2024-09-12T11:33:37.737086Z",
"url": "https://files.pythonhosted.org/packages/73/aa/e63ffb5e4c1378f57e5c6629059eacab760ecaa3c3219a86400326d9f7d4/TDCRPy-2.0.9-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"blake2b_256": "ece1cf1e37671e2105af63e3c3e83cb1a51a2276b08f1459b938970f717a966a",
"md5": "391b73f787c684d9417916c27ebc1bee",
"sha256": "19654cb77ccbaca4ec31f587df2f32174cc258dc325bbf68f7159addd65f08f7"
},
"downloads": -1,
"filename": "tdcrpy-2.0.9.tar.gz",
"has_sig": false,
"md5_digest": "391b73f787c684d9417916c27ebc1bee",
"packagetype": "sdist",
"python_version": "source",
"requires_python": ">=3.11",
"size": 22037071,
"upload_time": "2024-09-12T11:33:40",
"upload_time_iso_8601": "2024-09-12T11:33:40.812369Z",
"url": "https://files.pythonhosted.org/packages/ec/e1/cf1e37671e2105af63e3c3e83cb1a51a2276b08f1459b938970f717a966a/tdcrpy-2.0.9.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2024-09-12 11:33:40",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "RomainCoulon",
"github_project": "TDCRPy",
"travis_ci": false,
"coveralls": false,
"github_actions": true,
"lcname": "tdcrpy"
}
JSON
