From b3f23b77571fb39299c31e54888a291e22fbf484 Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Fri, 21 Oct 2022 14:35:55 +0200 Subject: [PATCH] FEAT: determine polarization with averaged polarimeter (#233) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * DOC: do not scroll interpolation code * FEAT: compute uncertainty increase on polarization * FIX: remove plus sign from |̅α| column --- .cspell.json | 1 + docs/uncertainties.ipynb | 2 +- docs/zz.polarization-fit.ipynb | 314 +++++++++++++++++++++++++++++---- 3 files changed, 278 insertions(+), 39 deletions(-) diff --git a/.cspell.json b/.cspell.json index 7252faf3..6cccc130 100644 --- a/.cspell.json +++ b/.cspell.json @@ -119,6 +119,7 @@ "mathbf", "mystnb", "nbsp", + "numref", "prereleased", "rgba", "shared_yaxes", diff --git a/docs/uncertainties.ipynb b/docs/uncertainties.ipynb index 5e8cfb65..6b64bdc0 100644 --- a/docs/uncertainties.ipynb +++ b/docs/uncertainties.ipynb @@ -1728,7 +1728,7 @@ " Δα = 1e3 * syst_weighted_alpha_diff[:, i]\n", " abs_Δα = 1e3 * syst_weighted_alpha_norm_diff[i]\n", " src += Rf\" \\textbf{{{i}}}\"\n", - " src += Rf\" & {α[0]:+.1f} & {α[1]:+.1f} & {α[2]:+.1f} & {abs_α:+.1f}\"\n", + " src += Rf\" & {α[0]:+.1f} & {α[1]:+.1f} & {α[2]:+.1f} & {abs_α:.1f}\"\n", " if i != 0:\n", " src += Rf\" & {Δα[0]:+.1f} & {Δα[1]:+.1f} & {Δα[2]:+.1f} & {abs_Δα:+.1f}\"\n", " src += R\" \\\\\" \"\\n\"\n", diff --git a/docs/zz.polarization-fit.ipynb b/docs/zz.polarization-fit.ipynb index e07dc8bb..c68dfc7a 100644 --- a/docs/zz.polarization-fit.ipynb +++ b/docs/zz.polarization-fit.ipynb @@ -10,7 +10,7 @@ "\n", "::::{only} html\n", ":::{margin}\n", - "This notebook has a `zz.` prefix because it has to be executed _after_ the polarimeter fields are exported in {doc}`/uncertainties`.\n", + "This notebook has a `zz.` prefix, because it has to be executed _after_ the polarimeter fields are exported in {doc}`/uncertainties`.\n", ":::\n", "::::\n", "\n", @@ -37,9 +37,11 @@ "source": [ "from __future__ import annotations\n", "\n", + "import json\n", "import logging\n", "import os\n", "from functools import lru_cache\n", + "from textwrap import dedent\n", "from warnings import filterwarnings\n", "\n", "import iminuit\n", @@ -81,7 +83,9 @@ "\n", "with $R$ the rotation matrix over the decay plane orientation, represented in Euler angles $\\left(\\phi, \\theta, \\chi\\right)$.\n", "\n", - "In this section, we show that it's possible to determine the polarization $\\vec{P}$ from a given intensity distribution $I$ of a $\\lambda_c$ decay if we the $\\vec\\alpha$ fields and the corresponding $I_0$ values of that $\\Lambda_c$ decay. We get $\\vec\\alpha$ and $I_0$ by interpolating the grid samples provided from {ref}`uncertainties:Exported distributions` using the method described in {ref}`appendix/serialization:Import and interpolate`.\n", + "In this section, we show that it's possible to determine the polarization $\\vec{P}$ from a given intensity distribution $I$ of a $\\lambda_c$ decay if we the $\\vec\\alpha$ fields and the corresponding $I_0$ values of that $\\Lambda_c$ decay. We get $\\vec\\alpha$ and $I_0$ by interpolating the grid samples provided from {ref}`uncertainties:Exported distributions` using the method described in {ref}`appendix/serialization:Import and interpolate`. We perform the same procedure with the averaged aligned polarimeter vector from {numref}`uncertainties:Average polarimetry values` in order to quantify the loss in precision when integrating over the Dalitz plane variables $\\tau$.\n", + "\n", + "## Polarized test distribution\n", "\n", "For this study, a phase space sample is uniformly generated over the Dalitz plane variables $\\tau$. The phase space sample is extended with uniform distributions over the decay plane angles $\\left(\\phi, \\theta, \\chi\\right)$, so that the phase space can be used to generate a hit-and-miss toy sample for a polarized intensity distribution." ] @@ -135,7 +139,6 @@ "code_prompt_show": "Code for interpolating α and I₀" }, "tags": [ - "scroll-input", "hide-input" ] }, @@ -290,6 +293,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "## Using the exported polarimeter grid\n", + "\n", "The generated distribution is now assumed to be a _measured distribution_ $I$ with unknown polarization $\\vec{P}$. It is shown below that the actual $\\vec{P}$ with which the distribution was generated can be found by performing a fit on Eq. {eq}`eq:master.intensity`. This is done with [`iminuit`](https://iminuit.rtfd.io), starting with a certain 'guessed' value for $\\vec{P}$ as initial parameters.\n", "\n", "To avoid having to generate a hit-and-miss intensity test distribution, the parameters $\\vec{P} = \\left(P_x, P_y, P_z\\right)$ are optimized with regard to a **weighted negative log likelihood estimator**:\n", @@ -321,7 +326,7 @@ "execution_count": null, "metadata": { "mystnb": { - "code_prompt_show": "Fit polarization to generated distribution" + "code_prompt_show": "Fit polarization with full polarimeter field" }, "tags": [ "hide-input" @@ -329,7 +334,7 @@ }, "outputs": [], "source": [ - "def perform_fit(phsp: DataSample, model_id: int) -> iminuit.Minuit:\n", + "def perform_full_fit(phsp: DataSample, model_id: int) -> iminuit.Minuit:\n", " I0 = interpolate_intensity(phsp, model_id)\n", " alpha = interpolate_polarimetry_field(phsp, model_id)\n", "\n", @@ -344,8 +349,8 @@ " return optimizer.migrad()\n", "\n", "\n", - "FIT_RESULTS = [\n", - " perform_fit(PHSP, i)\n", + "FULL_FIT_RESULTS = [\n", + " perform_full_fit(PHSP, i)\n", " for i in tqdm(range(17), desc=\"Performing fits\", disable=NO_TQDM)\n", "]" ] @@ -366,7 +371,7 @@ }, "outputs": [], "source": [ - "FIT_RESULTS[0]" + "FULL_FIT_RESULTS[0]" ] }, { @@ -383,26 +388,40 @@ }, "outputs": [], "source": [ - "P_fit_values = 100 * np.array([[p.value for p in fit.params] for fit in FIT_RESULTS])\n", - "P_fit_nominal = P_fit_values[0]\n", - "P_max = (P_fit_values[1:] - P_fit_nominal).max(axis=0)\n", - "P_min = (P_fit_values[1:] - P_fit_nominal).min(axis=0)\n", - "np.testing.assert_array_almost_equal(P_fit_nominal, 100 * np.array(P), decimal=2)\n", - "\n", - "\n", - "def render_p(i: int) -> str:\n", - " return f\"{P_fit_nominal[i]:+.2f}_{{{P_min[i]:+.2f}}}^{{{P_max[i]:+.2f}}}\"\n", + "def extract_polarization(fit_result: iminuit.Minuit) -> tuple[float, float, float]:\n", + " return tuple(p.value for p in fit_result.params)\n", + "\n", + "\n", + "def render_fit_results(\n", + " fit_results: list[iminuit.Minuit], compare: bool = False\n", + ") -> str:\n", + " P_fit_values = 100 * np.array([extract_polarization(fit) for fit in fit_results])\n", + " P_fit_nominal = P_fit_values[0]\n", + " P_max = (P_fit_values[1:] - P_fit_nominal).max(axis=0)\n", + " P_min = (P_fit_values[1:] - P_fit_nominal).min(axis=0)\n", + " if compare:\n", + " np.testing.assert_array_almost_equal(\n", + " P_fit_nominal, 100 * np.array(P), decimal=2\n", + " )\n", + "\n", + " def render_p(i: int) -> str:\n", + " return f\"{P_fit_nominal[i]:+.2f}_{{{P_min[i]:+.2f}}}^{{{P_max[i]:+.2f}}}\"\n", + "\n", + " src = Rf\"\"\"\n", + " \\begin{{array}}{{ccc}}\n", + " P_x &=& {render_p(0)} \\\\\n", + " P_y &=& {render_p(1)} \\\\\n", + " P_z &=& {render_p(2)} \\\\\n", + " \\end{{array}}\n", + " \"\"\"\n", + " return dedent(src).strip()\n", "\n", "\n", "src = Rf\"\"\"\n", "The polarization $\\vec{{P}}$ is determined to be (in %):\n", "\n", "$$\n", - "\\begin{{array}}{{ccc}}\n", - "P_x &=& {render_p(0)} \\\\\n", - "P_y &=& {render_p(1)} \\\\\n", - "P_z &=& {render_p(2)} \\\\\n", - "\\end{{array}}\n", + "{render_fit_results(FULL_FIT_RESULTS, compare=True)}\n", "$$\n", "\n", "with the upper and lower sign being the systematic extrema uncertainties as determined by\n", @@ -413,7 +432,12 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [] + }, "source": [ "This is to be compared with the model uncertainties reported by [^1]:\n", "\n", @@ -441,22 +465,236 @@ }, "outputs": [], "source": [ - "src = R\"\"\"\n", - "\\begin{array}{r|ccc|ccc}\n", - " \\textbf{Model} & \\mathbf{P_x} & \\mathbf{P_y} & \\mathbf{P_z}\n", - " & \\mathbf{\\Delta P_x} & \\mathbf{\\Delta P_y} & \\mathbf{\\Delta P_z} \\\\\n", - " \\hline\n", + "def render_all_polarizations(fit_results: list[iminuit.Minuit]) -> Latex:\n", + " src = R\"\"\"\n", + " \\begin{array}{r|ccc|ccc}\n", + " \\textbf{Model} & \\mathbf{P_x} & \\mathbf{P_y} & \\mathbf{P_z}\n", + " & \\mathbf{\\Delta P_x} & \\mathbf{\\Delta P_y} & \\mathbf{\\Delta P_z} \\\\\n", + " \\hline\n", + " \"\"\"\n", + " P_fit_values = np.array([extract_polarization(fit) for fit in fit_results])\n", + " P_fit_values *= 100\n", + " Px_nom, Py_nom, Pz_nom = P_fit_values[0]\n", + " for i, (Px, Py, Pz) in enumerate(P_fit_values):\n", + " src += Rf\" \\textbf{{{i}}} & {Px:+.2f} & {Py:+.2f} & {Pz:+.1f} & \"\n", + " if i != 0:\n", + " src += Rf\"{Px-Px_nom:+.2f} & {Py-Py_nom:+.2f} & {Pz-Pz_nom:+.2f}\"\n", + " src += R\" \\\\\" \"\\n\"\n", + " src += R\"\\end{array}\"\n", + " src = dedent(src).strip()\n", + " return Latex(src)\n", + "\n", + "\n", + "render_all_polarizations(FULL_FIT_RESULTS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using the averaged polarimeter vector\n", + "\n", + "Equation {eq}`eq:master.intensity` requires knowledge about the aligned polarimeter field $\\vec\\alpha(\\tau)$ and intensity distribution $I_0(\\tau)$ over all kinematic variables $\\tau$. It is, however, also possible to compute the differential decay rate from the averaged polarimeter vector $\\vec{\\overline{\\alpha}}$ (see {ref}`uncertainties:Average polarimetry values`). The equivalent formula to Eq. {eq}`eq:master.intensity` is:\n", + "\n", + "$$\n", + "\\frac{8\\pi^2}{\\Gamma}\\frac{\\mathrm{d}^3 \\Gamma}{\\mathrm{d}\\phi\\,\\mathrm{d}\\cos\\theta\\,\\mathrm{d}\\chi} = 1+\\sum_{i,j}P_i R_{ij}(\\phi, \\theta, \\chi) \\overline{\\alpha}_j\\,,\n", + "$$ (eq:I.alpha.averaged)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Code for computing differential decay rate" + }, + "tags": [ + "hide-input", + "scroll-input" + ] + }, + "outputs": [], + "source": [ + "def get_averaged_polarimeters(polar: bool = False) -> jnp.ndarray:\n", + " with open(\"_static/export/averaged-polarimeter-vectors.json\") as f:\n", + " json_data = json.load(f)\n", + " data = json_data[\"systematics\"]\n", + " typ = \"polar\" if polar else \"cartesian\"\n", + " return jnp.array([data[typ][i] for i in \"xyz\"]).T\n", + "\n", + "\n", + "def compute_differential_decay_rate(\n", + " Px: float,\n", + " Py: float,\n", + " Pz: float,\n", + " averaged_alpha: jnp.array,\n", + " phsp: DataSample,\n", + ") -> jnp.array:\n", + " P = jnp.array([Px, Py, Pz])\n", + " R = compute_rotation_matrix(phsp)\n", + " return 1 + jnp.einsum(\"i,ij...,j...->...\", P, R, averaged_alpha)\n", + "\n", + "\n", + "AVERAGED_POLARIMETERS = get_averaged_polarimeters()\n", + "assert AVERAGED_POLARIMETERS.shape == (17, 3)\n", + "\n", + "diff_rate = compute_differential_decay_rate(*P, AVERAGED_POLARIMETERS[0], PHSP)\n", + "assert diff_rate.shape == (N_EVENTS,)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use this equation along with Eq. {eq}`eq:weighted-nll` to determine $\\vec{P}$ with {class}`~iminuit.Minuit`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "mystnb": { + "code_prompt_show": "Fit polarization with averaged polarimeter" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "def perform_averaged_fit(phsp: DataSample, model_id: int) -> iminuit.Minuit:\n", + " averaged_alpha = AVERAGED_POLARIMETERS[model_id]\n", + "\n", + " def weighted_nll(Px: float, Py: float, Pz: float) -> float:\n", + " I_new = compute_differential_decay_rate(Px, Py, Pz, averaged_alpha, phsp)\n", + " I_new /= jnp.sum(I_new)\n", + " estimator_value = -jnp.sum(jnp.log(I_new) * I)\n", + " return estimator_value\n", + "\n", + " optimizer = iminuit.Minuit(weighted_nll, *P_guess)\n", + " optimizer.errordef = optimizer.LIKELIHOOD\n", + " return optimizer.migrad()\n", + "\n", + "\n", + "AVERAGED_FIT_RESULTS = [\n", + " perform_averaged_fit(PHSP, i)\n", + " for i in tqdm(range(17), desc=\"Performing fits\", disable=NO_TQDM)\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Show Minuit fit result for nominal model" + }, + "tags": [ + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "AVERAGED_FIT_RESULTS[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "src = Rf\"\"\"\n", + "Using the averaged polarimeter vector $\\vec{{\\overline{{\\alpha}}}}$, the\n", + "polarization $\\vec{{P}}$ is determined to be (in %):\n", + "\n", + "$$\n", + "{render_fit_results(AVERAGED_FIT_RESULTS)}\\,.\n", + "$$\n", + "\n", + "\n", + "The polarimeter values for each model are (in %):\n", "\"\"\"\n", - "Px_nom, Py_nom, Pz_nom = P_fit_nominal\n", - "for i, (Px, Py, Pz) in enumerate(P_fit_values):\n", - " src += Rf\" \\textbf{{{i}}} & {Px:+.2f} & {Py:+.2f} & {Pz:+.1f} & \"\n", - " if i != 0:\n", - " src += Rf\"{Px-Px_nom:+.2f} & {Py-Py_nom:+.2f} & {Pz-Pz_nom:+.2f}\"\n", - " src += R\" \\\\\" \"\\n\"\n", - " del Px, Py, Pz\n", - "del Px_nom, Py_nom, Pz_nom\n", - "src += R\"\\end{array}\"\n", - "Latex(src)" + "Markdown(src)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "render_all_polarizations(AVERAGED_FIT_RESULTS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Increase in uncertainties\n", + "\n", + "When the polarization is determined with the averaged aligned polarimeter vector $\\vec{\\overline{\\alpha}}$ instead of the aligned polarimeter vector field $\\vec\\alpha(\\tau)$ over all Dalitz variables $\\tau$, the uncertainty is expected to increase by a factor $S_0 / \\overline{S}_0 \\approx 3$, with:\n", + "\n", + "$$\n", + "S_0^2 = 3 \\int I_0 \\left|\\vec{\\alpha}\\right|^2 \\mathrm{d}^n \\tau \\,\\big /\\, \\int I_0\\,\\mathrm{d}^n \\tau \\\\\n", + "\\overline{S}_0^2 = 3(\\overline{\\alpha}_x^2+\\overline{\\alpha}_y^2+\\overline{\\alpha}_z^2)\\,.\n", + "$$ (eq:s0.integrals)\n", + "\n", + "The following table shows the maximal deviation (systematic uncertainty) of the determined polarization $\\vec{P}$ for each alternative model. The second and third column indicate the systematic uncertainty (in %) as determined with the full vector field and with the averaged vector, respectively." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "def render_uncertainty_increase() -> Latex:\n", + " src = R\"\"\"\n", + " \\begin{array}{c|cc|c}\n", + " \\sigma_\\mathrm{{model}} & \\vec\\alpha(\\tau) & \\vec{\\overline{\\alpha}} & \\text{factor} \\\\\n", + " \\hline\n", + " \"\"\"\n", + " src = dedent(src)\n", + " full_P = np.array([extract_polarization(fit) for fit in FULL_FIT_RESULTS])\n", + " aver_P = np.array([extract_polarization(fit) for fit in AVERAGED_FIT_RESULTS])\n", + " for i, xyz in enumerate(\"xyz\"):\n", + " full_sigma = 100 * (full_P[:, i] - full_P[0, i]).max()\n", + " aver_sigma = 100 * (aver_P[:, i] - aver_P[0, i]).max()\n", + " src += Rf\" P_{xyz} & {full_sigma:.2f} & {aver_sigma:.2f}\"\n", + " src += Rf\" & {aver_sigma/full_sigma:.1f} \\\\\" \"\\n\"\n", + " src += R\"\\end{array}\"\n", + " return Latex(src)\n", + "\n", + "\n", + "render_uncertainty_increase()" ] } ],