Mass-model effects on transmission-coefficient uncertainty¶
This notebook explores how mass-model choices propagate into transmission coefficients and related uncertainty estimates. It is a good next step once you are comfortable with the built-in uncertainty-quantified optical potentials.
1import numpy as np
2from IPython.display import Math, display
3from periodictable import elements
4
5from jitr.utils import mass
1def neutron_fermi_energy_KD(A, Z):
2 return -11.2814 + 0.02646 * A
3
4
5def proton_fermi_energy_KD(A, Z):
6 return -8.4075 + 0.01378 * A
Compare mass-model inputs to Koning-Delaroche¶
1A, Z = (154, 56)
1name_core = str(elements[Z].symbol)
2display(Math(f"^{{{A}}} \\rm{{{name_core}}}"))
\[\displaystyle ^{154} \rm{Ba}\]
1for m in mass.__MASS_MODELS__:
2 Ef = mass.neutron_fermi_energy(A, Z, model=m)
3 print(f"{m:7}: {Ef[0]:1.4f} +/- {Ef[1]}")
BMA : -3.6700 +/- 2.928163874630591
UNEDF2 : -3.5264 +/- None
BCPM : -3.6764 +/- None
dz31 : -3.4939 +/- None
hfb31 : -3.6514 +/- None
ws4rbf : -3.4174 +/- None
FRDM2012: -3.5514 +/- None
SLY4 : -3.4264 +/- None
SKP : -3.6264 +/- None
ame2020: nan +/- None
SKM : -4.1414 +/- None
UNEDF0 : -3.6114 +/- None
D1M : -3.6014 +/- None
SVMIN : -3.6514 +/- None
UNEDF1 : -3.4964 +/- None
HFB24 : -3.7864 +/- None
1neutron_fermi_energy_KD(A, Z)
-7.20656
1for m in mass.__MASS_MODELS__:
2 Ef = mass.proton_fermi_energy(A, Z, model=m)
3 print(f"{m:7}: {Ef[0]:1.4f} +/- {Ef[1]}")
BMA : -13.1027 +/- 2.8786616888136276
UNEDF2 : -13.6486 +/- None
BCPM : -13.2686 +/- None
dz31 : -13.5521 +/- None
hfb31 : -12.6936 +/- None
ws4rbf : -13.8941 +/- None
FRDM2012: -13.5936 +/- None
SLY4 : -13.8936 +/- None
SKP : -13.1386 +/- None
ame2020: nan +/- None
SKM : -12.5986 +/- None
UNEDF0 : -13.3486 +/- None
D1M : -13.6836 +/- None
SVMIN : -13.4186 +/- None
UNEDF1 : -13.7586 +/- None
HFB24 : -13.2636 +/- None
1proton_fermi_energy_KD(A, Z)
-6.28538
Compute transmission coefficients¶
Let’s use jitr to calculate UQ’ed transmission coefficients using KDUQ Fermi energies and those from a variety of mass models
1from tqdm import tqdm
2
3import jitr
1neutron = (1, 0)
2target = (A, Z)
1# com_energy_grid = np.logspace(-1, 1.3, 100)
2com_energy_grid = np.linspace(0.01, 5, 100)
3range_fm = jitr.utils.interaction_range(A)
4lmax = 2
5print(range_fm)
6.432130093158435
1reaction = jitr.reactions.Reaction(target=target, projectile=neutron, process="EL")
1def set_up_grid(core, com_energy_grid, mass_model="BMA"):
2 solvers = []
3 mn = jitr.utils.constants.MASS_N
4 mcore = jitr.utils.mass.mass(*core, model=mass_model)[0]
5 for _i, Ecm in enumerate(tqdm(com_energy_grid)):
6 kinematics = jitr.utils.kinematics.classical_kinematics_cm(mcore, mn, Ecm)
7 a = range_fm * kinematics.k + np.pi
8 N = 50
9 assert N > jitr.utils.suggested_basis_size(a)
10
11 solvers.append(
12 jitr.xs.elastic.IntegralWorkspace(
13 reaction=reaction,
14 kinematics=kinematics,
15 channel_radius_fm=a / kinematics.k,
16 solver=jitr.rmatrix.Solver(N),
17 lmax=lmax,
18 smatrix_abs_tol=0,
19 )
20 )
21 return solvers
1solvers = set_up_grid(target, com_energy_grid, mass_model="BMA")
100%|███████████████████████████████████████████████████████████████████████████████| 100/100 [00:16<00:00, 5.89it/s]
Run the uncertainty propagation¶
1samples = jitr.optical_potentials.kduq.get_samples(neutron)
KDUQ¶
1tcoeff_kduq = np.zeros((com_energy_grid.size, samples.shape[0], 2, lmax + 1))
2for j, sample in enumerate(tqdm(samples)):
3 for i, _Ecm in enumerate(com_energy_grid):
4 rgrid = solvers[i].radial_grid()
5 central_params, spin_orbit_params, coulomb_params = (
6 jitr.optical_potentials.kduq.calculate_params(
7 neutron,
8 target,
9 solvers[i].kinematics.Elab,
10 *sample,
11 )
12 )
13
14 tplus, tminus = solvers[i].transmission_coefficients(
15 jitr.optical_potentials.kduq.central(rgrid, *central_params),
16 jitr.optical_potentials.kduq.spin_orbit(rgrid, *spin_orbit_params),
17 jitr.optical_potentials.kduq.coulomb_charged_sphere(rgrid, *coulomb_params),
18 )
19 tcoeff_kduq[i, j, 0, :] = tplus
20 tcoeff_kduq[i, j, 1, :] = tminus
97%|████████████████████████████████████████████████████████████████████████████▌ | 403/416 [00:44<00:01, 11.49it/s]/home/kyle/umich/jitr/src/jitr/optical_potentials/kduq.py:467: RuntimeWarning: overflow encountered in exp
d2 = d2_0 + d2_A / (1 + np.exp((A - d2_A3) / d2_A2))
100%|███████████████████████████████████████████████████████████████████████████████| 416/416 [00:45<00:00, 9.14it/s]
KDUQ with BMA Fermi energies¶
1param_order = jitr.optical_potentials.kduq.get_param_names(neutron)
2idx_Ef0 = [i for i in range(len(param_order)) if param_order[i] == "Ef_0"][0]
3idx_EfA = [i for i in range(len(param_order)) if param_order[i] == "Ef_A"][0]
4
5idx_Ef0, idx_EfA
(0, 1)
1# set Fermi energy to BMA
2Ef = mass.neutron_fermi_energy(A, Z, model="BMA")
1tcoeff_kduq_bma = np.zeros((com_energy_grid.size, samples.shape[0], 2, lmax + 1))
2for j, sample in enumerate(tqdm(samples)):
3 # reset Ef
4 sample[idx_EfA] = 0
5 sample[idx_Ef0] = np.random.normal(loc=Ef[0], scale=Ef[1]) # sample from BMA
6
7 for i, _Ecm in enumerate(com_energy_grid):
8 rgrid = solvers[i].radial_grid()
9 central_params, spin_orbit_params, coulomb_params = (
10 jitr.optical_potentials.kduq.calculate_params(
11 neutron,
12 target,
13 solvers[i].kinematics.Elab,
14 *sample,
15 )
16 )
17
18 tplus, tminus = solvers[i].transmission_coefficients(
19 jitr.optical_potentials.kduq.central(rgrid, *central_params),
20 jitr.optical_potentials.kduq.spin_orbit(rgrid, *spin_orbit_params),
21 jitr.optical_potentials.kduq.coulomb_charged_sphere(rgrid, *coulomb_params),
22 )
23 tcoeff_kduq_bma[i, j, 0, :] = tplus
24 tcoeff_kduq_bma[i, j, 1, :] = tminus
100%|███████████████████████████████████████████████████████████████████████████████| 416/416 [00:34<00:00, 11.97it/s]
Extract confidence intervals¶
1def get_confidence_intervals(tcoeffs):
2 ci = []
3 for l in range(3):
4 ci.append([])
5 ci[l].append(
6 (
7 np.percentile(tcoeffs[:, :, 0, l], 16, axis=1),
8 np.percentile(tcoeffs[:, :, 0, l], 84, axis=1),
9 )
10 )
11 if l > 0:
12 ci[l].append(
13 (
14 np.percentile(tcoeffs[:, :, 1, l], 16, axis=1),
15 np.percentile(tcoeffs[:, :, 1, l], 84, axis=1),
16 )
17 )
18 return ci
1kduq_ci = get_confidence_intervals(tcoeff_kduq)
2kduq_bma_ci = get_confidence_intervals(tcoeff_kduq_bma)
Plot the results¶
1kduq_color = "tab:blue"
2bma_color = "tab:orange"
1from matplotlib import pyplot as plt
1fig, (ax1, ax2, ax3) = plt.subplots(
2 figsize=(4, 6), nrows=3, ncols=1, dpi=300, sharex=True
3)
4
5ax1.fill_between(
6 com_energy_grid,
7 kduq_bma_ci[0][0][0],
8 kduq_bma_ci[0][0][1],
9 color=bma_color,
10 alpha=0.3,
11 label="KDUQ + BMA",
12)
13ax1.fill_between(
14 com_energy_grid,
15 kduq_ci[0][0][0],
16 kduq_ci[0][0][1],
17 color=kduq_color,
18 alpha=0.3,
19 label="KDUQ",
20)
21
22ax2.fill_between(
23 com_energy_grid,
24 kduq_bma_ci[1][0][0],
25 kduq_bma_ci[1][0][1],
26 color=bma_color,
27 alpha=0.2,
28 hatch="|||",
29)
30ax2.fill_between(
31 com_energy_grid,
32 kduq_ci[1][0][0],
33 kduq_ci[1][0][1],
34 color=kduq_color,
35 alpha=0.2,
36 hatch="|||",
37)
38ax2.fill_between(
39 com_energy_grid,
40 kduq_bma_ci[1][1][0],
41 kduq_bma_ci[1][1][1],
42 color=bma_color,
43 alpha=0.2,
44 hatch="---",
45)
46ax2.fill_between(
47 com_energy_grid,
48 kduq_ci[1][1][0],
49 kduq_ci[1][1][1],
50 color=kduq_color,
51 alpha=0.2,
52 hatch="---",
53)
54
55ax3.fill_between(
56 com_energy_grid,
57 kduq_bma_ci[2][0][0],
58 kduq_bma_ci[2][0][1],
59 color=bma_color,
60 alpha=0.2,
61 hatch="|||",
62)
63ax3.fill_between(
64 com_energy_grid,
65 kduq_ci[2][0][0],
66 kduq_ci[2][0][1],
67 color=kduq_color,
68 alpha=0.2,
69 hatch="|||",
70)
71ax3.fill_between(
72 com_energy_grid,
73 kduq_bma_ci[2][1][0],
74 kduq_bma_ci[2][1][1],
75 color=bma_color,
76 alpha=0.2,
77 hatch="---",
78)
79ax3.fill_between(
80 com_energy_grid,
81 kduq_ci[2][1][0],
82 kduq_ci[2][1][1],
83 color=kduq_color,
84 alpha=0.2,
85 hatch="---",
86)
87
88# style legend
89swave_style = [
90 ax3.fill_between(
91 [-1, 0],
92 [-1, 0],
93 [-1, 0],
94 color="k",
95 alpha=0.2,
96 ),
97]
98lwave_style = [
99 ax3.fill_between(
100 [-1, 0],
101 [-1, 0],
102 [-1, 0],
103 color="k",
104 alpha=0.2,
105 hatch="---",
106 ),
107 ax3.fill_between(
108 [-1, 0],
109 [-1, 0],
110 [-1, 0],
111 color="k",
112 alpha=0.2,
113 hatch="|||",
114 ),
115]
116
117
118style_legend2 = ax2.legend(
119 lwave_style,
120 [r"$p_{1/2}$", r"$p_{3/2}$"],
121 ncols=2,
122 # loc="upper right",
123 fontsize=12,
124 framealpha=1,
125)
126style_legend3 = ax3.legend(
127 lwave_style,
128 [r"$d_{1/2}$", r"$d_{3/2}$"],
129 ncols=2,
130 loc="lower right",
131 fontsize=12,
132 framealpha=1,
133)
134fig.legend(loc="upper center", ncols=2, framealpha=0.5)
135
136ax1.set_ylim([0, 1])
137ax2.set_ylim([0, 1])
138ax3.set_ylim([0, 1])
139
140ax1.set_ylabel(r"$T_{jl}(E)$ for $s_{1/2}$")
141ax2.set_ylabel(r"$T_{jl}(E)$ for $p$-wave")
142ax3.set_ylabel(r"$T_{jl}(E)$ for $d$-wave")
143
144plt.xlabel(r"$E_{cm}$ [MeV]")
145plt.xlim([0, 5])
(0.0, 5.0)
