BEHAVIOR: add systematic (detector) uncertainties (#184)

* FIX: swap uncertainties D(1600)/L(1600)

Co-authored-by: Remco de Boer <[email protected]>

data/model-definitions.yaml

@@ -43,10 +43,10 @@ Default amplitude model:
     AiL(1520)1: 0.044324 ± 0.022276 ± 0.12 ± 0.02
     ArL(1520)2: -0.160687 ± 0.075604 ± 0.69 ± 0.03
     AiL(1520)2: 1.498833 ± 0.032468 ± 1.3 ± 0.0
-    ArL(1600)1: 4.840649 ± 0.080342 ± 9.8 ± 0.2
-    AiL(1600)1: 3.082786 ± 0.313304 ± 1.4 ± 0.2
-    ArL(1600)2: -6.971233 ± 0.180278 ± 7.2 ± 0.2
-    AiL(1600)2: 0.842435 ± 0.300438 ± 2.1 ± 0.1
+    ArL(1600)1: 4.840649 ± 0.080342 ± 5.0 ± 0.1
+    AiL(1600)1: 3.082786 ± 0.313304 ± 3.7 ± 0.1
+    ArL(1600)2: -6.971233 ± 0.180278 ± 8.7 ± 0.1
+    AiL(1600)2: 0.842435 ± 0.300438 ± 2.0 ± 0.2
     ArL(1670)1: -0.339585 ± 0.013871 ± 0.35 ± 0.01
     AiL(1670)1: -0.144678 ± 0.022591 ± 0.22 ± 0.02
     ArL(1670)2: -0.570978 ± 0.059400 ± 0.46 ± 0.02
@@ -65,10 +65,10 @@ Default amplitude model:
     AiD(1232)1: 3.051805 ± 0.264359 ± 1.4 ± 0.1
     ArD(1232)2: -12.987193 ± 0.433180 ± 12 ± 0
     AiD(1232)2: 4.528336 ± 0.493430 ± 3.7 ± 0.3
-    ArD(1600)1: 11.401585 ± 0.437580 ± 5.0 ± 0.1
-    AiD(1600)1: -3.125511 ± 0.473191 ± 3.7 ± 0.1
-    ArD(1600)2: 6.729211 ± 0.283501 ± 8.7 ± 0.1
-    AiD(1600)2: -1.002383 ± 0.331137 ± 2.0 ± 0.2
+    ArD(1600)1: 11.401585 ± 0.437580 ± 9.8 ± 0.2
+    AiD(1600)1: -3.125511 ± 0.473191 ± 1.4 ± 0.2
+    ArD(1600)2: 6.729211 ± 0.283501 ± 7.2 ± 0.2
+    AiD(1600)2: -1.002383 ± 0.331137 ± 2.1 ± 0.1
     ArD(1700)1: 10.378280 ± 0.253080 ± 7.2 ± 0.3
     AiD(1700)1: 1.434872 ± 0.485532 ± 5.2 ± 0.2
     ArD(1700)2: 12.874102 ± 0.233827 ± 13 ± 0

docs/uncertainties.ipynb

@@ -401,16 +401,16 @@
     "np.testing.assert_almost_equal(\n",
     "    np.array(stat_intensities[bootstrap_selection, event_selection]),\n",
     "    [\n",
-    "        [2177.066401, 5459.644607, 6256.06364],\n",
-    "        [2263.337352, 5662.47171, 6339.399809],\n",
-    "        [2036.602322, 5124.725992, 6120.279915],\n",
-    "        [2060.930775, 5270.754879, 6131.571211],\n",
+    "        [2181.869539, 5457.018803, 6252.331785],\n",
+    "        [2272.483121, 5683.29221, 6329.830268],\n",
+    "        [2036.791629, 5087.462553, 6125.017938],\n",
+    "        [2053.192849, 5263.222729, 6142.861638],\n",
     "    ],\n",
     "    decimal=6,\n",
     ")\n",
     "np.testing.assert_almost_equal(\n",
     "    np.array(stat_polarimetry_norm[bootstrap_selection, 0]),\n",
-    "    [0.711835, 0.708065, 0.716237, 0.735825],\n",
+    "    [0.71179, 0.709137, 0.71447, 0.734856],\n",
     "    decimal=6,\n",
     ")\n",
     "del bootstrap_selection, event_selection"
@@ -1361,7 +1361,7 @@
    "source": [
     "np.testing.assert_almost_equal(\n",
     "    np.array(stat_decay_rates[\"L(1405)\"][:3]),\n",
-    "    [0.0796555, 0.0804976, 0.0789567],\n",
+    "    [0.0799202, 0.0809439, 0.0789926],\n",
     ")\n",
     "np.testing.assert_almost_equal(\n",
     "    np.array(syst_decay_rates[\"L(1405)\"]),\n",

src/polarimetry/lhcb/__init__.py

@@ -10,6 +10,7 @@
 import re
 import sys
 from copy import deepcopy
+from math import sqrt
 from pathlib import Path
 from typing import Generic, Iterable, TypeVar, overload
 
@@ -377,7 +378,7 @@ def _to_value_with_uncertainty(str_value: str) -> MeasuredParameter[float]:
     >>> par
     MeasuredParameter(value=0.94, hesse=0.042, model=0.35, systematic=0.04)
     >>> par.uncertainty
-    0.042
+    0.058
     """
     float_values = tuple(float(s) for s in str_value.split(" ± "))
     if len(float_values) == 2:
@@ -438,8 +439,14 @@ class MeasuredParameter(Generic[ParameterType]):
 
     @property
     def uncertainty(self) -> ParameterType:
-        # Will implement quadratic sum of uncertainties here
-        return self.hesse
+        if self.systematic is None:
+            return self.hesse
+        if isinstance(self.value, float):
+            return sqrt(self.hesse**2 + self.systematic**2)
+        return complex(
+            sqrt(self.hesse.real**2 + self.systematic.real**2),
+            sqrt(self.hesse.imag**2 + self.systematic.imag**2),
+        )
 
 
 def get_conversion_factor(