Open in Colab

Differential cross-section uncertainty propagation with KDUQ

KDUQ is an uncertainty-quantified global optical potential. This notebook shows how to sample that posterior information, run the solver repeatedly, and summarize the resulting uncertainty in differential cross sections.

1# import stuff for nice plotting
2import matplotlib.pyplot as plt
3import numpy as np
4from tqdm import tqdm
5
6import jitr

1from jitr.optical_potentials import kduq

 1#  elastic reaction
 2target = (54, 26)
 3proton = (1, 1)
 4neutron = (1, 0)
 5projectile = proton
 6
 7reaction = jitr.reactions.ElasticReaction(
 8    target=target,
 9    projectile=projectile,
10)
11
12# energy
13Elab = 35
14kinematics = reaction.kinematics(Elab)
15
16# for plotting differential xs
17angles = np.linspace(0.1, np.pi, 180)
18
19# Lagrange Mesh
20core_solver = jitr.rmatrix.Solver(40)

Set up the solver for differential cross sections

 1# get kinematics and parameters for this experiment
 2
 3a = jitr.utils.interaction_range(target[0]) * kinematics.k + np.pi * 2
 4N = jitr.utils.suggested_basis_size(a)
 5assert N < core_solver.kernel.quadrature.nbasis
 6channel_radius_fm = a / kinematics.k
 7
 8# build solver
 9solver = jitr.xs.elastic.DifferentialWorkspace.build_from_system(
10    reaction=reaction,
11    kinematics=kinematics,
12    channel_radius_fm=channel_radius_fm,
13    solver=core_solver,
14    lmax=50,
15    angles=angles,
16)

Sample the KDUQ posterior

1help(kduq.get_samples)

Help on function get_samples in module jitr.optical_potentials.kduq:

get_samples(projectile: tuple[int, int], posterior: str = 'federal') -> numpy.ndarray
    Get the posterior samples for the given projectile (neutron or
    proton) from the KDUQ Federal or Democratic posteriors.

    These samples are distinct from the original KD03 frequentist parameter
    set returned by :func:`get_kd03`.

    See [Pruitt, et al., 2023]
    (https://journals.aps.org/prc/pdf/10.1103/PhysRevC.107.014602) for
    details on the KDUQ posteriors.

    Args:
        projectile: tuple (Ap, Zp) of the projectile. Must be ``(1, 0)``
            for neutron or ``(1, 1)`` for proton.
        posterior: Which KDUQ posterior to return samples from. Must be
            either ``"federal"`` or ``"democratic"``. Defaults to ``"federal"``.

    Returns:
        An array of shape ``(NUM_POSTERIOR_SAMPLES, num_params)`` containing
        the posterior samples for the given projectile.

1kduq_samples = kduq.get_samples(projectile)

Set up the potential

Let’s take a look at the functionality built into jitr:

1omp = kduq.KDUQ(projectile)
2help(omp)

Help on KDUQ in module jitr.optical_potentials.kduq object:

class KDUQ(jitr.optical_potentials.omp.SingleChannelOpticalModel)
 |  KDUQ(projectile: tuple)
 |
 |  Koning-Delaroche Uncertainty Quantification (KDUQ) optical
 |  potential model.
 |
 |  Method resolution order:
 |      KDUQ
 |      jitr.optical_potentials.omp.SingleChannelOpticalModel
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(self, projectile: tuple)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  evaluate(self, rgrid: float | numpy.ndarray, reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]
 |      Evaluate the KDUQ central, spin-orbit, and Coulomb terms.
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'
 |      Call self as a function.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object

1help(jitr.optical_potentials.kduq.central)

Help on function central in module jitr.optical_potentials.kduq:

central(r: float | numpy.ndarray, Vv: float, Rv: float, av: float, Wv: float, Rwv: float, awv: float, Wd: float, Rd: float, ad: float) -> complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]]
    Koning-Delaroche central terms at a given energy.

    This matches Eq. (7) in Koning and Delaroche (2003).

    Args:
        r: The radius at which to evaluate the potential.
        Vv: The real central depth.
        Rv: The real central radius parameter.
        av: The real central diffuseness parameter.
        Wv: The imaginary volume depth.
        Rwv: The imaginary volume radius parameter.
        awv: The imaginary volume diffuseness parameter.
        Wd: The imaginary surface depth.
        Rd: The imaginary surface radius parameter.
        ad: The imaginary surface diffuseness parameter.

1help(jitr.optical_potentials.kduq.spin_orbit)

Help on function spin_orbit in module jitr.optical_potentials.kduq:

spin_orbit(r: float | numpy.ndarray, Vso: float, Rso: float, aso: float, Wso: float, Rwso: float, awso: float) -> complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]]
    Koning-Delaroche spin-orbit terms at a given energy.

    This matches Eq. (7) in Koning and Delaroche (2003).

    Args:
        r: The radius at which to evaluate the potential.
        Vso: The real spin-orbit depth.
        Rso: The real spin-orbit radius parameter.
        aso: The real spin-orbit diffuseness parameter.
        Wso: The imaginary spin-orbit depth.
        Rwso: The imaginary spin-orbit radius parameter.
        awso: The imaginary spin-orbit diffuseness parameter.

1help(jitr.optical_potentials.potential_forms.coulomb_charged_sphere)

Help on function coulomb_charged_sphere in module jitr.optical_potentials.potential_forms:

coulomb_charged_sphere(r: 'ArrayOrScalar', zz: 'float', r_c: 'float') -> 'ArrayOrScalar'
    Return the Coulomb potential of a uniformly charged sphere.

Run the calculation

This should demonstrate how the KDUQ class (omp) is built to interface with solver, which is an instance of jitr.xs.elastic.DifferentialWorkspace. This is how calculations are done in jitr.

1kduq_xs = np.zeros((len(angles), kduq.NUM_POSTERIOR_SAMPLES))
2rgrid = solver.radial_grid()
3
4for j, sample in enumerate(tqdm(kduq_samples)):
5    central_term, spin_orbit_term, coulomb_term = omp(
6        rgrid, reaction, kinematics, *sample
7    )
8    xs = solver.xs(central_term, spin_orbit_term, coulomb_term)
9    kduq_xs[:, j] = xs.dsdo / solver.rutherford

 89%|████████▉ | 372/416 [00:10<00:00, 216.23it/s]/home/kyle/umich/jitr/src/jitr/optical_potentials/kduq.py:470: RuntimeWarning: overflow encountered in exp
  d2 = d2_0 + d2_A / (1 + np.exp((A - d2_A3) / d2_A2))
100%|██████████| 416/416 [00:10<00:00, 40.80it/s]

Compute intervals and plot results

1kduq_pred_post = np.percentile(kduq_xs, [16, 84], axis=1)

 1fig, ax = plt.subplots(1, 1, figsize=(6, 4))
 2ax.fill_between(
 3    angles * 180 / np.pi,
 4    kduq_pred_post[0],  # / solver.rutherford,
 5    kduq_pred_post[1],  # / solver.rutherford,
 6    color="#ff4500",
 7    hatch="|-|-",
 8    alpha=0.2,
 9)
10plt.xlabel(r"$\theta$ [deg]")
11plt.ylabel(r"$\frac{d \sigma}{d\Omega} / \frac{d \sigma_{R}}{d\Omega} $")

Text(0, 0.5, '$\\frac{d \\sigma}{d\\Omega} / \\frac{d \\sigma_{R}}{d\\Omega} $')
../../_images/f149a3fc7f506d1ee07f610dd398d96ace3d4f62ad20cdb8496a9c087d3fd562.png 
1# NBVAL_CHECK_OUTPUT
2print(f"{kduq_pred_post[1][9]:1.6f}")

1.149638

In fact, KDUQ is just one definition of a global optical potential. jitr also has others built in:

1help(jitr.optical_potentials.wlh.WLH)

Help on class WLH in module jitr.optical_potentials.wlh:

class WLH(jitr.optical_potentials.omp.SingleChannelOpticalModel)
 |  WLH(projectile: tuple)
 |
 |  The Whitehead-Lim-Holt global optical potential for nucleon-nucleus
 |  scattering.
 |
 |  Method resolution order:
 |      WLH
 |      jitr.optical_potentials.omp.SingleChannelOpticalModel
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(self, projectile: tuple)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  evaluate(self, rgrid: float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]], reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]
 |      Evaluate the central, spin-orbit, and Coulomb terms of the WLH
 |      potential on the given radial grid for the specified reaction and
 |      kinematics, using the provided potential parameters.
 |
 |      Args:
 |          rgrid: Radial coordinate or grid in fm.
 |          reaction: Reaction for which to calculate the parameters.
 |          kinematics: Kinematics of the reaction channel.
 |          *params: Parameters of the WLH potential. See
 |              `Whitehead et al., 2021
 |              <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.127.182502>`_
 |              for details.
 |
 |      Returns:
 |          ``(U_central, U_spin_orbit, U_coulomb)`` evaluated on the radial grid.
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'
 |      Call self as a function.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object

1help(jitr.optical_potentials.chuq.CHUQ)

Help on class CHUQ in module jitr.optical_potentials.chuq:

class CHUQ(jitr.optical_potentials.omp.SingleChannelOpticalModel)
 |  Chapel-Hill Uncertainty Quantification (CHUQ) optical
 |  potential model.
 |
 |  Note that CH89 is Lane consistent, so the same parameters can be
 |  used for both neutron and proton projectiles.
 |
 |  Method resolution order:
 |      CHUQ
 |      jitr.optical_potentials.omp.SingleChannelOpticalModel
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(self)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  evaluate(self, rgrid: float | numpy.ndarray, reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]
 |      Evaluate the CHUQ central, spin-orbit, and Coulomb terms.
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'
 |      Call self as a function.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object

In fact, these are all derived from jitr.optical_potentials.SingleChannelOpticalModel, which sets the interface between the model for the effective interaction and the solver. You can come up with your own optical model and derive your own class from jitr.optical_potentials.SingleChannelOpticalModel to plug it into jitr!

1help(jitr.optical_potentials.SingleChannelOpticalModel)

Help on class SingleChannelOpticalModel in module jitr.optical_potentials.omp:

class SingleChannelOpticalModel(builtins.object)
 |  SingleChannelOpticalModel(params: 'list[str]') -> 'None'
 |
 |  Base class for local single-channel optical potentials.
 |
 |  Methods defined here:
 |
 |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'
 |      Call self as a function.
 |
 |  __init__(self, params: 'list[str]') -> 'None'
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  evaluate(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'
 |      Evaluate central, spin-orbit, and Coulomb terms on ``rgrid``.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object