Uncertainty quantification of partial wave transmission coefficients

This notebook extends the uncertainty-propagation workflow to angular reaction observables. It shows how the same model choices that affect transmission coefficients flow through to angle-dependent reaction outputs.

1import numpy as np
2from IPython.display import Math, display
3from periodictable import elements

1from matplotlib import pyplot as plt

Compare mass-model inputs to Koning-Delaroche

1A, Z = (24, 12)

1name_core = str(elements[Z].symbol)
2display(Math(f"^{{{A}}} \\rm{{{name_core}}}"))
\[\displaystyle ^{24} \rm{Mg}\]

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)
2proton = (1, 1)
3projectile = proton
4target = (A, Z)

1# we have 416 samples from the KDUQ posterior
2kduq_omp_samples = jitr.optical_potentials.kduq.get_samples(proton)

1# com_energy_grid = np.logspace(-1, 1.3, 100)
2lab_energy_grid = np.array([65, 200])
3range_fm = 15
4lmax = 20

1reaction = jitr.reactions.Reaction(target=target, projectile=projectile, process="EL")

 1def set_up_grid(core, lab_energy_grid):
 2    solvers = []
 3    for _i, Elab in enumerate(tqdm(lab_energy_grid)):
 4        kinematics = reaction.kinematics(Elab)
 5        a = range_fm * kinematics.k + np.pi / 2
 6        N = jitr.utils.suggested_basis_size(a)
 7        solvers.append(
 8            jitr.xs.elastic.IntegralWorkspace(
 9                reaction=reaction,
10                kinematics=kinematics,
11                channel_radius_fm=a / kinematics.k,
12                solver=jitr.rmatrix.Solver(N),
13                lmax=lmax,
14                smatrix_abs_tol=0,
15            )
16        )
17    return solvers

1solvers = set_up_grid(target, lab_energy_grid)

100%|███████████████████████████████████████████████████████████████████████████████████| 2/2 [00:14<00:00,  7.29s/it]

Run the uncertainty propagation

KDUQ

 1tcoeff_kduq = np.zeros((lab_energy_grid.size, kduq_omp_samples.shape[0], 2, lmax + 1))
 2for j, sample in enumerate(tqdm(kduq_omp_samples)):
 3    for i, _Elab in enumerate(lab_energy_grid):
 4        rgrid = solvers[i].radial_grid()
 5        central_params, spin_orbit_params, coulomb_params = (
 6            jitr.optical_potentials.kduq.calculate_params(
 7                projectile,
 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 [08:09<00:23,  1.81s/it]/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 [08:26<00:00,  1.22s/it]

 1data = [np.zeros(lmax + 1)] * 8
 2names = [
 3    "E=65 MeV, j = l + 1/2",
 4    "E=65 MeV, j = l - 1/2",
 5    "E=65 MeV, j = l + 1/2, err",
 6    "E=65 MeV, j = l - 1/2, err",
 7    "E=200 MeV, j = l + 1/2",
 8    "E=200 MeV, j = l - 1/2",
 9    "E=200 MeV, j = l + 1/2, err",
10    "E=200 MeV, j = l - 1/2, err",
11]

1from pandas import DataFrame as df

1data = df.from_dict(dict(zip(names, data, strict=False)))

1plus_color = "tab:blue"
2minus_color = "tab:orange"

1ci_plus = np.percentile(tcoeff_kduq[0, :, 0, :], 50, axis=0)
2ci_plus_errs = np.percentile(tcoeff_kduq[0, :, 0, :], 84, axis=0) - np.percentile(
3    tcoeff_kduq[0, :, 0, :], 16, axis=0
4)

1ci_minus = np.percentile(tcoeff_kduq[0, :, 1, :], 50, axis=0)
2ci_minus_errs = np.percentile(tcoeff_kduq[0, :, 1, :], 84, axis=0) - np.percentile(
3    tcoeff_kduq[0, :, 1, :], 16, axis=0
4)

1ci_minus[0] = None
2ci_minus_errs[0] = None

1data["E=65 MeV, j = l + 1/2"] = ci_plus
2data["E=65 MeV, j = l + 1/2, err"] = ci_plus_errs
3data["E=65 MeV, j = l - 1/2"] = ci_minus
4data["E=65 MeV, j = l - 1/2, err"] = ci_minus_errs

 1fig = plt.figure(figsize=(6, 4))
 2ls = np.arange(lmax + 1)
 3plt.errorbar(
 4    ls,
 5    ci_plus * (2 * ls + 1),
 6    ci_plus_errs * (2 * ls + 1),
 7    linestyle="none",
 8    marker="^",
 9    label="$j = l + 1/2$",
10)
11
12plt.errorbar(
13    ls[1:],
14    ci_minus[1:] * (2 * ls[1:] + 1),
15    ci_minus_errs[1:] * (2 * ls[1:] + 1),
16    linestyle="none",
17    marker="v",
18    label="$j = l - 1/2$",
19)
20
21# plt.yscale("log")
22plt.xticks([0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
23plt.xlabel(r"$\ell$")
24plt.ylabel(r"$(2  \ell +1) \mathcal{T}_{lj}$")
25plt.title(r"$^{24}$Mg(p,p) 65 MeV", fontsize=14)
26plt.legend()
27plt.tight_layout()
../../_images/90e92a34a642b08328e96a0dd60186c551a5d81510f2933a103a40ed9c68cfe4.png 
1ci_plus = np.percentile(tcoeff_kduq[1, :, 0, :], 50, axis=0)
2ci_plus_errs = np.percentile(tcoeff_kduq[1, :, 0, :], 84, axis=0) - np.percentile(
3    tcoeff_kduq[1, :, 0, :], 16, axis=0
4)

1ci_minus = np.percentile(tcoeff_kduq[1, :, 1, :], 50, axis=0)
2ci_minus_errs = np.percentile(tcoeff_kduq[1, :, 1, :], 84, axis=0) - np.percentile(
3    tcoeff_kduq[1, :, 1, :], 16, axis=0
4)
5ci_minus[0] = None
6ci_minus_errs[0] = None

1data["E=200 MeV, j = l + 1/2"] = ci_plus
2data["E=200 MeV, j = l + 1/2, err"] = ci_plus_errs
3data["E=200 MeV, j = l - 1/2"] = ci_minus
4data["E=200 MeV, j = l - 1/2, err"] = ci_minus_errs

 1fig = plt.figure(figsize=(6, 4))
 2plt.errorbar(
 3    ls,
 4    ci_plus * (2 * ls + 1),
 5    ci_plus_errs * (2 * ls + 1),
 6    linestyle="none",
 7    marker="^",
 8    label="$j = l + 1/2$",
 9)
10
11plt.errorbar(
12    ls[1:],
13    ci_minus[1:] * (2 * ls[1:] + 1),
14    ci_minus_errs[1:] * (2 * ls[1:] + 1),
15    linestyle="none",
16    marker="v",
17    label="$j = l - 1/2$",
18)
19
20# plt.yscale("log")
21plt.xticks([0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
22plt.xlabel(r"$\ell$")
23plt.ylabel(r"$(2  \ell +1) \mathcal{T}_{lj}$")
24plt.title(r"$^{24}$Mg(p,p) 200 MeV", fontsize=14)
25plt.legend()
26plt.tight_layout()
../../_images/7263f46b7ee6522a59eca580f017a8fe99eecb327ebd7357e9dba775ed3b43e4.png 
1# NBVAL_CHECK_OUTPUT
2data
E=65 MeV, j = l + 1/2 E=65 MeV, j = l - 1/2 E=65 MeV, j = l + 1/2, err E=65 MeV, j = l - 1/2, err E=200 MeV, j = l + 1/2 E=200 MeV, j = l - 1/2 E=200 MeV, j = l + 1/2, err E=200 MeV, j = l - 1/2, err
0 0.732306 NaN 0.103491 NaN 0.793631 NaN 0.083665 NaN
1 0.714397 0.733912 0.106325 0.101622 0.781721 0.806950 0.086771 0.084428
2 0.707490 0.740207 0.101053 0.093985 0.765326 0.809628 0.089057 0.087703
3 0.678968 0.730161 0.100371 0.088427 0.741524 0.808922 0.090451 0.094959
4 0.656383 0.725234 0.095936 0.075356 0.709459 0.804422 0.101630 0.100298
5 0.681355 0.689816 0.074087 0.082002 0.664013 0.797320 0.111231 0.107043
6 0.589953 0.538968 0.083585 0.095931 0.605166 0.781423 0.127492 0.114318
7 0.358919 0.329351 0.067499 0.065663 0.537147 0.751475 0.154021 0.130638
8 0.169756 0.165448 0.035843 0.034335 0.465712 0.698140 0.178390 0.146219
9 0.072202 0.073235 0.017226 0.016429 0.393139 0.621629 0.178045 0.154005
10 0.029473 0.030809 0.007992 0.007903 0.319778 0.527126 0.162889 0.147182
11 0.012078 0.012668 0.003880 0.003637 0.249148 0.422926 0.125635 0.130248
12 0.004944 0.005216 0.001846 0.001790 0.187534 0.319209 0.089697 0.102278
13 0.002036 0.002163 0.000890 0.000856 0.136752 0.229855 0.061229 0.077375
14 0.000841 0.000896 0.000424 0.000425 0.096277 0.158128 0.041666 0.055007
15 0.000356 0.000377 0.000211 0.000207 0.065568 0.105852 0.027758 0.039233
16 0.000152 0.000159 0.000100 0.000095 0.043796 0.069592 0.018947 0.025934
17 0.000065 0.000068 0.000047 0.000046 0.029016 0.044560 0.013011 0.017157
18 0.000028 0.000029 0.000023 0.000022 0.019050 0.028528 0.008826 0.011763
19 0.000012 0.000012 0.000011 0.000011 0.012400 0.018159 0.006208 0.007886
20 0.000005 0.000005 0.000005 0.000005 0.008033 0.011565 0.004265 0.005293