From 0aee4ebad16fb567c127dc550f73e31c8ccbe6a8 Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Fri, 21 Oct 2022 17:04:57 +0200 Subject: [PATCH] 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 --- docs/zz.polarization-fit.ipynb | 66 ++++++++++++++++++---------------- julia/notebooks/eulerangles.jl | 22 ++++++------ 2 files changed, 47 insertions(+), 41 deletions(-) diff --git a/docs/zz.polarization-fit.ipynb b/docs/zz.polarization-fit.ipynb index c68dfc7a..bc21a1ea 100644 --- 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 $\\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 index ee063982..249da8df 100644 --- 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) ≈ χ