FIX: take absolute value of systematic extrema (#234)

* DOC: make "factor" column gray * MAINT: apply `pre-commit` formatting to Julia * MAINT: implement `extract_polarizations()` * MAINT: prepend `SYST_` prefixes
ComPWA · Oct 21, 2022 · 0aee4eb · 0aee4eb
1 parent b3f23b7
commit 0aee4eb
Show file tree

Hide file tree

Showing 2 changed files with 47 additions and 41 deletions.
diff --git a/docs/zz.polarization-fit.ipynb b/docs/zz.polarization-fit.ipynb
@@ -349,7 +349,7 @@
     "    return optimizer.migrad()\n",
     "\n",
     "\n",
-    "FULL_FIT_RESULTS = [\n",
+    "SYST_FIT_RESULTS_FULL = [\n",
     "    perform_full_fit(PHSP, i)\n",
     "    for i in tqdm(range(17), desc=\"Performing fits\", disable=NO_TQDM)\n",
     "]"
@@ -371,7 +371,7 @@
    },
    "outputs": [],
    "source": [
-    "FULL_FIT_RESULTS[0]"
+    "SYST_FIT_RESULTS_FULL[0]"
    ]
   },
   {
@@ -388,24 +388,27 @@
    },
    "outputs": [],
    "source": [
+    "def extract_polarizations(fit_results: list[iminuit.Minuit]) -> np.ndarray:\n",
+    "    return np.array([extract_polarization(fit) for fit in fit_results])\n",
+    "\n",
+    "\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",
+    "    fit_results: list[iminuit.Minuit],\n",
+    "    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",
+    "    P_syst = 100 * extract_polarizations(fit_results)\n",
+    "    P_nominal = P_syst[0]\n",
+    "    P_max = (P_syst[1:] - P_nominal).max(axis=0)\n",
+    "    P_min = (P_syst[1:] - P_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",
+    "        np.testing.assert_array_almost_equal(P_nominal, 100 * np.array(P), decimal=2)\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",
+    "        return f\"{P_nominal[i]:+.2f}_{{{P_min[i]:+.2f}}}^{{{P_max[i]:+.2f}}}\"\n",
     "\n",
     "    src = Rf\"\"\"\n",
     "    \\begin{{array}}{{ccc}}\n",
@@ -421,7 +424,7 @@
     "The polarization $\\vec{{P}}$ is determined to be (in %):\n",
     "\n",
     "$$\n",
-    "{render_fit_results(FULL_FIT_RESULTS, compare=True)}\n",
+    "{render_fit_results(SYST_FIT_RESULTS_FULL, compare=True)}\n",
     "$$\n",
     "\n",
     "with the upper and lower sign being the systematic extrema uncertainties as determined by\n",
@@ -485,7 +488,7 @@
     "    return Latex(src)\n",
     "\n",
     "\n",
-    "render_all_polarizations(FULL_FIT_RESULTS)"
+    "render_all_polarizations(SYST_FIT_RESULTS_FULL)"
    ]
   },
   {
@@ -538,11 +541,12 @@
     "    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",
+    "SYST_AVERAGED_POLARIMETERS = get_averaged_polarimeters()\n",
+    "assert SYST_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,)"
+    "diff_rate = compute_differential_decay_rate(*P, SYST_AVERAGED_POLARIMETERS[0], PHSP)\n",
+    "assert diff_rate.shape == (N_EVENTS,)\n",
+    "del diff_rate"
    ]
   },
   {
@@ -566,7 +570,7 @@
    "outputs": [],
    "source": [
     "def perform_averaged_fit(phsp: DataSample, model_id: int) -> iminuit.Minuit:\n",
-    "    averaged_alpha = AVERAGED_POLARIMETERS[model_id]\n",
+    "    averaged_alpha = SYST_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",
@@ -579,7 +583,7 @@
     "    return optimizer.migrad()\n",
     "\n",
     "\n",
-    "AVERAGED_FIT_RESULTS = [\n",
+    "SYST_FIT_RESULTS_AVERAGED = [\n",
     "    perform_averaged_fit(PHSP, i)\n",
     "    for i in tqdm(range(17), desc=\"Performing fits\", disable=NO_TQDM)\n",
     "]"
@@ -601,7 +605,7 @@
    },
    "outputs": [],
    "source": [
-    "AVERAGED_FIT_RESULTS[0]"
+    "SYST_FIT_RESULTS_AVERAGED[0]"
    ]
   },
   {
@@ -622,7 +626,7 @@
     "polarization&nbsp;$\\vec{{P}}$ is determined to be (in %):\n",
     "\n",
     "$$\n",
-    "{render_fit_results(AVERAGED_FIT_RESULTS)}\\,.\n",
+    "{render_fit_results(SYST_FIT_RESULTS_AVERAGED)}\\,.\n",
     "$$\n",
     "\n",
     "\n",
@@ -644,7 +648,7 @@
    },
    "outputs": [],
    "source": [
-    "render_all_polarizations(AVERAGED_FIT_RESULTS)"
+    "render_all_polarizations(SYST_FIT_RESULTS_AVERAGED)"
    ]
   },
   {
@@ -678,18 +682,20 @@
    "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",
+    "    \\begin{array}{c|ccc}\n",
+    "      \\sigma_\\mathrm{{model}}\n",
+    "      & \\vec\\alpha(\\tau) & \\vec{\\overline{\\alpha}} & \\color{gray}{\\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",
+    "    syst_P_full = 100 * extract_polarizations(SYST_FIT_RESULTS_FULL)\n",
+    "    syst_P_aver = 100 * extract_polarizations(SYST_FIT_RESULTS_AVERAGED)\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 += f\"  P_{xyz}\"\n",
+    "        syst_sigma_full = np.abs(syst_P_full[:, i] - syst_P_full[0, i]).max()\n",
+    "        syst_sigma_aver = np.abs(syst_P_aver[:, i] - syst_P_aver[0, i]).max()\n",
+    "        src += Rf\" & {syst_sigma_full:.2f} & {syst_sigma_aver:.2f}\"\n",
+    "        src += Rf\" & \\color{{gray}}{{{syst_sigma_aver/syst_sigma_full:.1f}}} \\\\\" \"\\n\"\n",
     "    src += R\"\\end{array}\"\n",
     "    return Latex(src)\n",
     "\n",

diff --git a/julia/notebooks/eulerangles.jl b/julia/notebooks/eulerangles.jl
@@ -20,7 +20,7 @@ begin
     using Parameters
     #
     using ThreeBodyDecay
-	using Lc2ppiKModelLHCb 
+	using Lc2ppiKModelLHCb
 end
 
 # ╔═╡ 049886e5-0839-4dce-9621-32cce58132f5
@@ -62,10 +62,10 @@ begin
     Rz(θ) = x -> Rz(x, θ)
     Ry(θ) = x -> Ry(x, θ)
 	R(ϕ, θ, χ) = x -> R(x, ϕ, θ, χ)
-	
+
 	Ry(ps::NamedTuple{(:p1,:p2,:p3)}, θ) =
 		(; p1=Ry(ps.p1,θ), p2=Ry(ps.p2,θ), p3=Ry(ps.p3,θ))
-	# 
+	#
 	function R(ps::NamedTuple{(:p1,:p2,:p3)}, ϕ, θ, χ)
 		@unpack p1, p2, p3 = ps
 	    p1′ = p1 |> Rz(χ) |> Ry(θ) |> Rz(ϕ)
@@ -136,7 +136,7 @@ ps0 = constructinvariants(σs0) |> Ry(π)
 # ╔═╡ 36526f74-2e4c-44c1-8097-67553277ed83
 let
 	@unpack p4p, p4pi, p4k = pointvariables
-	isapprox(ps0.p1, p4p[1:3]; atol=1e-8), 
+	isapprox(ps0.p1, p4p[1:3]; atol=1e-8),
 	isapprox(ps0.p2, p4pi[1:3]; atol=1e-8),
 	isapprox(ps0.p3, p4k[1:3]; atol=1e-8)
 end |> x->x=>(prod(x) && "✔")
@@ -156,7 +156,7 @@ begin
 	    #
 	    (ϕ1 ≈ ϕp) * 100 + (θ1 ≈ θp) * 10 + (ϕ23 ≈ χK)
 	end
-		# 
+		#
 	σsv = flatDalitzPlotSample(ms; Nev=10)
 	(checkconsistency.(σsv; ϕ1=-0.3, θ1=1.3, ϕ23=-0.4) .== 111) |>
 		x->x=>(prod(x) && "✔")
@@ -167,7 +167,7 @@ md"""
 #### Check-3:
 The angle-mapping relations are checked:
 ```math
-\begin{align} 
+\begin{align}
     &&\bar{\theta}_K &= \theta_{23}\,,\\ \nonumber
     \theta_{\Lambda} &= \pi - \zeta^0_{1(2)}\,,&\theta_p^{\Lambda} &= -(\pi - \theta_{31})\,, \\ \nonumber
     \theta_{\Delta} &= -(\pi-\zeta^0_{3(1)})\,,&\theta_p^{\Delta} &= \theta_{12}\,.
@@ -180,7 +180,7 @@ begin # computed the DPD angles
 	θ12 = acos(cosθ12(σs0, ms²))
 	θ23 = acos(cosθ23(σs0, ms²))
 	θ31 = acos(cosθ31(σs0, ms²))
-	# 
+	#
 	ζ⁰12 = acos(cosζ(wr(1, 2, 0), σs0, ms²))
 	ζ⁰31 = acos(cosζ(wr(3, 1, 0), σs0, ms²))
 end ;
@@ -232,7 +232,7 @@ function eulerfromalgebra(_pD, _pB,
 		_pP = [1,0,0] # x axes
 		)
 	pP, pX, pD, pB = _pP[1:3], _pX[1:3] ,_pD[1:3], _pB[1:3]
-	# 
+	#
 	ϕ = atan(-(pP × pX) · pD, (-(pP × pX) × (-pX) ./ norm(pX)) · pD)
 	θ = acos((pX · pD) / norm(pX) / norm(pD))
 	χ = atan((pX × pD) · pB, ((pX × pD) × -pD ./ norm(pD)) · pB)
@@ -266,8 +266,8 @@ The following cell validates the transformation of the angles under parity flip
 \end{array}
 \right.
 \\
-\theta&\to \pi-\theta \\ 
-\chi &\to 
+\theta&\to \pi-\theta \\
+\chi &\to
 \left\{
 \begin{array}{}
 -\pi-\chi &\text{ for } \chi < 0,\\
@@ -319,7 +319,7 @@ let # 3 angles x 10 tries
 		pB′ = -pB |> Ry(π)
 		pX′ = -[0,0,-1] |> Ry(π)
 		pP′ = -[1,0,0] |> Ry(π)
-		# 
+		#
 		ϕ, θ, χ = eulerfromalgebra(pD′, pB′, pX′, pP′)
 		#
 		mapϕ′(angles.ϕ1) ≈ ϕ && mapθ′(angles.θ1) ≈ θ && mapχ′(angles.ϕ23) ≈ χ