{ "cells": [ { "cell_type": "markdown", "id": "ba38fa7e", "metadata": {}, "source": [ "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/beykyle/jitr/blob/main/examples/notebooks/qepn_analyzing_power.ipynb)\n", "\n", "On Colab the first `import exfor_tools` downloads and unpacks the EXFOR database, which takes a few minutes." ] }, { "cell_type": "code", "execution_count": 1, "id": "135b5879", "metadata": {}, "outputs": [], "source": [ "# Google Colab setup: installs jitr and this notebook's extra dependencies.\n", "# Does nothing when the notebook is run anywhere else.\n", "import sys\n", "\n", "if \"google.colab\" in sys.modules:\n", " import subprocess\n", "\n", " subprocess.run(\n", " [sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"jitr\", \"exfor-tools\"],\n", " check=True,\n", " )" ] }, { "cell_type": "markdown", "id": "204d393d", "metadata": {}, "source": [ "# Analyzing power of the quasi-elastic $(\\vec{p},n)$ reaction\n", "\n", "In 1976 Gosset, Mayer and Escudié measured the differential analyzing power $A_y(\\theta)$ of\n", "the quasi-elastic $(p,n)$ reaction to the isobaric analog state (IAS) on nine targets, with a\n", "22.8 MeV polarized proton beam at Saclay\n", "([Phys. Rev. C **14**, 878](https://doi.org/10.1103/PhysRevC.14.878)).\n", "Those data are EXFOR entry [O0669](https://www-nds.iaea.org/exfor/servlet/X4sGetSubent?reqx=51740&subID=O0669).\n", "\n", "Their point was that $A_y$ is a far sharper probe of the Lane potential than the cross section is.\n", "In the macroscopic Lane model the nucleon-nucleus optical potential splits into isoscalar and\n", "isovector parts,\n", "\n", "$$ U = U_0 + \\frac{4}{A} U_1\\, \\vec{t} \\cdot \\vec{T}, $$\n", "\n", "which for the proton-target, neutron-analog and neutron-target channels gives\n", "\n", "$$ U_{pT} = U_0 - \\frac{2 T_0}{A} U_1, \\qquad\n", " U_{nA} = U_0 + \\frac{2 (T_0 - 1)}{A} U_1, \\qquad\n", " U_{nT} = U_0 + \\frac{2 T_0}{A} U_1, $$\n", "\n", "with $T_0 = (N - Z)/2$. The off-diagonal part of $U_1$ drives the charge-exchange transition\n", "between $|T_0 T_0\\rangle$ and its analog $|T_0, T_0 - 1\\rangle$. Because the transition carries\n", "$l = s = j = 0$, only two amplitudes survive — one conserving and one inverting the projectile\n", "spin projection — and the analyzing power is\n", "\n", "$$ A_y = \\frac{2\\,\\mathrm{Im}(Y X^*)}{|X|^2 + |Y|^2}. $$\n", "\n", "The spin-flip amplitude $Y$ is built from the *difference* between the $j = l + \\tfrac{1}{2}$ and\n", "$j = l - \\tfrac{1}{2}$ radial matrix elements. \n", "\n", "This notebook solves the two coupled Lane channels exactly with\n", "[`jitr.xs.lane_pn`](https://jitr.readthedocs.io/en/latest/autoapi/jitr/xs/lane_pn/index.html)\n", "rather than in DWBA, and propagates the posteriors of three uncertainty-quantified global optical\n", "potentials — KDUQ, WLH and CHUQ — plus a deterministic JLMB calculation, through to $A_y$ for all\n", "nine targets.\n", "\n", "Notw that, while $^{49}$Ti, $^{117}$Sn and $^{165}$Ho all have non-zero ground-state spin, the measured transitions transfer no angular momentum, so the target spin is a spectator and the spin-1/2 geometry used here is still reasonable." ] }, { "cell_type": "code", "execution_count": 2, "id": "d732cdc9", "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import periodictable\n", "\n", "import jitr\n", "from jitr.folding import ILDAFolder\n", "from jitr.folding.jlm import (\n", " lambda_v0,\n", " lambda_v1,\n", " lambda_vso,\n", " lambda_w0,\n", " lambda_w1,\n", " lambda_wso,\n", " potential_JLMB,\n", " spin_orbit_jlmb,\n", ")\n", "from jitr.optical_potentials import chuq, kduq, wlh\n", "from jitr.utils.density import density_table\n", "from jitr.xs import lane_pn\n", "\n", "# set JITR_QUICK=1 to thin the posterior draws; used by the CI notebook run\n", "QUICK = os.environ.get(\"JITR_QUICK\", \"0\") == \"1\"" ] }, { "cell_type": "markdown", "id": "ef5ca7f2", "metadata": {}, "source": [ "## Pulling the measurements from EXFOR\n", "\n", "Each of the nine targets is a subentry of O0669, tabulated as an analyzing power against the center-of-mass angle, together with the excitation energy of the analog state in the residual nucleus. That excitation energy is what fixes the outgoing neutron energy, so we read it from the entry to set the kinematics in the exit channels." ] }, { "cell_type": "code", "execution_count": 3, "id": "9de4a4b9-c983-4d82-89ca-4a1f1f9d9b23", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using database version X4-2024-12-31 located in: /home/kyle/db/exfor/unpack_exfor-2024/X4-2024-12-31\n", "1.3\n" ] } ], "source": [ "import exfor_tools\n", "\n", "print(exfor_tools.__version__)" ] }, { "cell_type": "code", "execution_count": 4, "id": "8d1eb7d8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "@article{JGosset1976_O0669001,\n", " author = {J.Gosset and B.Mayer and J.L.Escudie},\n", " title = {Quasielastic (P,N) Reactions Induced By Polarized Protons.},\n", " journal = {Physical Review, Part C, Nuclear Physics},\n", " year = {1976},\n", " volume = {14},\n", " pages = {878},\n", " note = {EXFOR reference: Physical Review, Part C, Nuclear Physics 14, 878 (1976); EXFOR subent: O0669001}\n", "}\n" ] } ], "source": [ "from exfor_tools import ExforEntry\n", "from exfor_tools.reaction import Reaction as ExforReaction\n", "\n", "# the nine targets of Table I, in the order the paper presents them\n", "TARGETS = [\n", " (49, 22),\n", " (56, 26),\n", " (64, 28),\n", " (70, 30),\n", " (90, 40),\n", " (96, 40),\n", " (117, 50),\n", " (165, 67),\n", " (208, 82),\n", "]\n", "PROTON, NEUTRON = (1, 1), (1, 0)\n", "\n", "\n", "def nuclide(A, Z):\n", " return rf\"$^{{{A}}}${periodictable.elements[Z].symbol}\"\n", "\n", "\n", "datasets = []\n", "for A, Z in TARGETS:\n", " entry = ExforEntry(\n", " \"O0669\",\n", " ExforReaction(\n", " target=(A, Z), projectile=PROTON, product=NEUTRON, residual=(A, Z + 1)\n", " ),\n", " quantity=\"Ay\",\n", " )\n", " (measurement,) = entry.measurements\n", " datasets.append(\n", " {\n", " \"A\": A,\n", " \"Z\": Z,\n", " \"name\": f\"{A}{periodictable.elements[Z].symbol}\",\n", " \"residual_name\": f\"{A}{periodictable.elements[Z + 1].symbol}\",\n", " \"label\": nuclide(A, Z),\n", " \"subentry\": measurement.subentry,\n", " \"E_lab\": measurement.Einc,\n", " \"E_IAS\": measurement.Ex,\n", " \"theta\": measurement.x,\n", " \"Ay\": measurement.y,\n", " \"Ay_err\": measurement.statistical_err,\n", " }\n", " )\n", "\n", "print(entry.bibtex())" ] }, { "cell_type": "code", "execution_count": 5, "id": "1d3f5bd4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " target residual subentry E_IAS [MeV] points angular range [deg]\n", "--------------------------------------------------------------------------\n", " 49Ti 49V O0669002 5.36 18 16.6 - 103.2\n", " 56Fe 56Co O0669003 3.51 23 15.7 - 125.1\n", " 64Ni 64Cu O0669004 6.70 19 16.7 - 106.4\n", " 70Zn 70Ga O0669005 8.12 20 16.0 - 114.3\n", " 90Zr 90Nb O0669006 5.03 22 15.4 - 124.3\n", " 96Zr 96Nb O0669007 11.07 18 14.4 - 98.7\n", " 117Sn 117Sb O0669008 11.18 17 15.0 - 94.9\n", " 165Ho 165Er O0669009 15.49 8 16.9 - 61.0\n", " 208Pb 208Bi O0669010 15.33 4 17.3 - 31.4\n", "\n", "149 points, all at E_p = 22.8 MeV\n" ] } ], "source": [ "header = f\"{'target':>8} {'residual':>10} {'subentry':>10} {'E_IAS [MeV]':>12} {'points':>7} {'angular range [deg]':>22}\"\n", "print(header)\n", "print(\"-\" * len(header))\n", "for d in datasets:\n", " angles = f\"{d['theta'][0]:.1f} - {d['theta'][-1]:.1f}\"\n", " print(\n", " f\"{d['name']:>8} {d['residual_name']:>10} {d['subentry']:>10} \"\n", " f\"{d['E_IAS']:>12.2f} {len(d['theta']):>7} {angles:>22}\"\n", " )\n", "print(\n", " f\"\\n{sum(len(d['theta']) for d in datasets)} points, all at E_p = {datasets[0]['E_lab']} MeV\"\n", ")" ] }, { "cell_type": "markdown", "id": "8b46aad9", "metadata": {}, "source": [ "## Compiling one coupled-channels workspace per target\n", "\n", "`jitr.xs.lane_pn.Workspace` solves the two-channel problem\n", "\n", "$$ \\left[T_l + U_{pp} - E_p\\right] u_p + U_1 u_n = 0, \\qquad\n", " \\left[T_l + U_{nn} - E_n\\right] u_n + U_1 u_p = 0 $$\n", "\n", "for every $(l, j)$, with an incoming wave in the proton channel only." ] }, { "cell_type": "code", "execution_count": 6, "id": "026ba475", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 49Ti E_n = 15.92 MeV a = 10.47 fm nbasis = 20\n", " 56Fe E_n = 13.78 MeV a = 10.66 fm nbasis = 20\n", " 64Ni E_n = 13.50 MeV a = 10.85 fm nbasis = 20\n", " 70Zn E_n = 13.11 MeV a = 10.99 fm nbasis = 20\n", " 90Zr E_n = 10.74 MeV a = 11.40 fm nbasis = 20\n", " 96Zr E_n = 10.99 MeV a = 11.52 fm nbasis = 20\n", " 117Sn E_n = 8.96 MeV a = 11.88 fm nbasis = 20\n", " 165Ho E_n = 6.05 MeV a = 12.58 fm nbasis = 25\n", " 208Pb E_n = 3.72 MeV a = 13.10 fm nbasis = 25\n" ] } ], "source": [ "E_LAB = 22.8\n", "LMAX = 25\n", "ANGLES = np.radians(np.linspace(1.0, 180.0, 180))\n", "\n", "solvers = {}\n", "for d in datasets:\n", " A, Z = d[\"A\"], d[\"Z\"]\n", " reaction = jitr.reactions.Reaction(\n", " target=(A, Z), projectile=PROTON, product=NEUTRON, residual=(A, Z + 1)\n", " )\n", " # \"El\" so that the neutron potential is evaluated for n + residual\n", " exit_reaction = jitr.reactions.Reaction(\n", " target=(A, Z + 1), projectile=NEUTRON, process=\"El\"\n", " )\n", " kinematics_entrance = reaction.kinematics(E_LAB)\n", " kinematics_exit = reaction.kinematics_exit(\n", " kinematics_entrance, residual_excitation_energy=d[\"E_IAS\"]\n", " )\n", "\n", " a = jitr.utils.interaction_range(A) * kinematics_entrance.k + 2 * np.pi\n", " nbasis = jitr.utils.suggested_basis_size(a)\n", " solver = solvers.setdefault(nbasis, jitr.rmatrix.Solver(nbasis))\n", "\n", " d.update(\n", " reaction=reaction,\n", " exit_reaction=exit_reaction,\n", " kinematics_entrance=kinematics_entrance,\n", " kinematics_exit=kinematics_exit,\n", " nbasis=nbasis,\n", " channel_radius_fm=a / kinematics_entrance.k,\n", " workspace=lane_pn.Workspace(\n", " reaction,\n", " kinematics_entrance,\n", " kinematics_exit,\n", " solver,\n", " ANGLES,\n", " LMAX,\n", " a / kinematics_entrance.k,\n", " ),\n", " )\n", "\n", "for d in datasets:\n", " print(\n", " f\"{d['name']:>8} E_n = {d['kinematics_exit'].Elab:5.2f} MeV \"\n", " f\"a = {d['channel_radius_fm']:5.2f} fm nbasis = {d['nbasis']:2d}\"\n", " )" ] }, { "cell_type": "markdown", "id": "76ec66a9", "metadata": {}, "source": [ "## Propagating the optical model posteriors\n", "\n", "Let's propagate the KDUQ, WLH and CHUQ posteriors through the $(p,n)$\n", "$A_y$.\n", "\n", "Nothing here is a fit to these data: the potentials were calibrated to elastic scattering, and the Lane transition potential is taken to be the isovector difference\n", "\n", "$$ U_1(E) = -\\left(U_n(E - Q - E_x) - U_p(E)\\right) \\frac{\\sqrt{|N - Z|}}{N - Z - 1}, $$\n", "\n", "evaluated separately for the central and spin-orbit terms — so the isovector spin-orbit form factor that the analyzing power is sensitive to comes entirely from the difference between the proton and neutron spin-orbit potentials of each global parameterization. Note that the neutron and proton potentials going into $U_1$ are evaluated at the energies of their respective channels.\n", "\n", "CHUQ deserves a word: CH89 is Lane-consistent by construction, so it is the one built-in whose proton and neutron potentials really are the $U_{pT}$ and $U_{nA}$ of the equations above, and whose difference therefore really is $U_1$. For KDUQ and WLH the same construction is an approximation, and the diagonal potentials may already absorb part of the coupling." ] }, { "cell_type": "code", "execution_count": 7, "id": "eb337795", "metadata": {}, "outputs": [], "source": [ "from tqdm import tqdm\n", "\n", "OMPS = {\n", " \"KDUQ\": (\n", " kduq.KDUQ(PROTON),\n", " kduq.KDUQ(NEUTRON),\n", " kduq.get_samples(PROTON),\n", " kduq.get_samples(NEUTRON),\n", " ),\n", " \"WLH\": (\n", " wlh.WLH(PROTON),\n", " wlh.WLH(NEUTRON),\n", " wlh.get_samples(PROTON),\n", " wlh.get_samples(NEUTRON),\n", " ),\n", " # CH89 is Lane consistent: one parameter set serves both projectiles\n", " \"CHUQ\": (chuq.CHUQ(), chuq.CHUQ(), chuq.get_samples(), chuq.get_samples()),\n", "}\n", "\n", "\n", "def propagate(omp_p, omp_n, samples_p, samples_n):\n", " n_draws = min(len(samples_p), len(samples_n))\n", " if QUICK:\n", " n_draws = min(n_draws, 50)\n", "\n", " bands = []\n", " for d in tqdm(datasets):\n", " workspace = d[\"workspace\"]\n", " rgrid = workspace.radial_grid()\n", " Ay = np.zeros((len(ANGLES), n_draws))\n", " dsdo = np.zeros((len(ANGLES), n_draws))\n", " for j in range(n_draws):\n", " U_p, U_p_so, U_coulomb = omp_p(\n", " rgrid, d[\"reaction\"], d[\"kinematics_entrance\"], *samples_p[j]\n", " )\n", " U_n, U_n_so, _ = omp_n(\n", " rgrid, d[\"exit_reaction\"], d[\"kinematics_exit\"], *samples_n[j]\n", " )\n", " _, S = workspace.rsmatrix(U_coulomb, U_p, U_p_so, U_n, U_n_so)\n", " observables = workspace.observables_from_smatrix(S)\n", " Ay[:, j] = observables.Ay\n", " dsdo[:, j] = observables.dsdo\n", " bands.append(\n", " {\n", " \"Ay\": np.percentile(Ay, [16, 84], axis=1),\n", " \"Ay_median\": np.median(Ay, axis=1),\n", " \"dsdo\": np.percentile(dsdo, [16, 84], axis=1),\n", " \"dsdo_median\": np.median(dsdo, axis=1),\n", " }\n", " )\n", " return n_draws, bands" ] }, { "cell_type": "code", "execution_count": 8, "id": "d2e1cc24-28f7-48f7-9c38-76b704f12b80", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " 0%| | 0/9 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "BANDS = {\"KDUQ\": \"#0072B2\", \"WLH\": \"#D55E00\", \"CHUQ\": \"#009E73\"}\n", "JLMB_INK = \"#333333\"\n", "angles_deg = np.degrees(ANGLES)\n", "\n", "\n", "def grid_figure(ylabel):\n", " fig, axes = plt.subplots(3, 3, figsize=(12, 10), sharex=True)\n", " fig.subplots_adjust(hspace=0.12, wspace=0.22)\n", " for ax in axes[-1]:\n", " ax.set_xlabel(r\"$\\theta_{\\rm c.m.}$ [deg]\")\n", " for ax in axes[:, 0]:\n", " ax.set_ylabel(ylabel)\n", " for ax in axes.flat:\n", " ax.set_xlim(0, 140)\n", " ax.tick_params(direction=\"in\", top=True, right=True)\n", " for side in ax.spines.values():\n", " side.set_linewidth(0.8)\n", " return fig, axes\n", "\n", "\n", "fig, axes = grid_figure(r\"$A_y$\")\n", "for i, (ax, d, curve) in enumerate(zip(axes.flat, datasets, jlmb, strict=True)):\n", " for name, color in BANDS.items():\n", " lo, hi = predictions[name][i][\"Ay\"]\n", " ax.fill_between(angles_deg, lo, hi, color=color, alpha=0.35, lw=0, label=name)\n", " # ax.plot(angles_deg, predictions[name][i][\"Ay_median\"], color=color, lw=1.2)\n", " ax.plot(angles_deg, curve.Ay, color=JLMB_INK, lw=1.2, ls=\"--\", label=\"JLMB\")\n", " ax.errorbar(\n", " d[\"theta\"],\n", " d[\"Ay\"],\n", " yerr=d[\"Ay_err\"],\n", " fmt=\"o\",\n", " ms=3.5,\n", " color=\"k\",\n", " lw=0.9,\n", " capsize=1.5,\n", " label=\"Gosset (1976)\",\n", " )\n", " ax.axhline(0.0, color=\"0.75\", lw=0.7, zorder=0)\n", " ax.set_ylim(-1.05, 1.05)\n", " ax.text(0.04, 0.06, d[\"label\"], transform=ax.transAxes, fontsize=13)\n", "\n", "axes[0, 0].legend(loc=\"upper right\", fontsize=8, frameon=False, ncols=2)\n", "fig.suptitle(r\"$(\\vec{p},n)$ analyzing power to the IAS at $E_p = 22.8$ MeV\", y=0.92)\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }