Update ratio uncertainty calculation (#225)

* Update ratio uncertainty calculation

* Update changelog

changelog.md

@@ -2,6 +2,7 @@
 
 ### [Latest]
 
+- Update ratio uncertainty calculation [!225](https://github.com/umami-hep/puma/pull/225)
 - Add NaN filtering for HL plots. Bump atlas-ftag-tools version to 0.1.11. [!220](https://github.com/umami-hep/puma/pull/220)
 - Making Under- and Overflow bins default in histograms [!218](https://github.com/umami-hep/puma/pull/218)
 

puma/histogram.py

@@ -166,7 +166,6 @@ def divide(self, other):
             numerator=self.hist,
             denominator=other.hist,
             numerator_unc=self.unc,
-            denominator_unc=other.unc,
             step=False,
         )
         # To use the matplotlib.step() function later on, the first bin is duplicated
@@ -178,7 +177,6 @@ def divide(self, other):
     def divide_data_mc(
         self,
         ref_hist: np.ndarray,
-        ref_unc: np.ndarray,
     ) -> tuple:
         """
         Similar as divide, but the second item doesn't need to be a histogram object
@@ -201,7 +199,6 @@ def divide_data_mc(
             numerator=self.hist,
             denominator=ref_hist,
             numerator_unc=self.unc,
-            denominator_unc=ref_unc,
             step=False,
         )
         # To use the matplotlib.step() function later on, the first bin is duplicated
@@ -718,7 +715,6 @@ def plot_ratios(self):
 
                 ratio, ratio_unc = elem.divide_data_mc(
                     ref_hist=self.stacked_dict["band"],
-                    ref_unc=self.stacked_dict["unc"],
                 )
 
             else:

puma/tests/test_histogram.py

@@ -46,7 +46,7 @@ def test_divide_after_plotting_no_norm(self):
         # duplicated in the ratio calculation (the first one is so to say not plotted)
         # Therefore, we also use duplicated bins here
         expected_ratio = np.array([3, 3, 2 / 3])
-        expected_ratio_unc = np.array([3.46410162, 3.46410162, 0.60858062])
+        expected_ratio_unc = np.array([1.73205081, 1.73205081, 0.47140452])
 
         np.testing.assert_almost_equal(expected_ratio, hist_1.divide(hist_2)[0])
         np.testing.assert_almost_equal(expected_ratio_unc, hist_1.divide(hist_2)[1])
@@ -64,7 +64,7 @@ def test_divide_after_plotting_norm(self):
         # duplicated in the ratio calculation (the first one is so to say not plotted)
         # Therefore, we also use duplicated bins here
         expected_ratio = np.array([2.4, 2.4, 0.53333333])
-        expected_ratio_unc = np.array([2.77128129, 2.77128129, 0.4868645])
+        expected_ratio_unc = np.array([1.38564065, 1.38564065, 0.37712362])
 
         np.testing.assert_almost_equal(expected_ratio, hist_1.divide(hist_2)[0])
         np.testing.assert_almost_equal(expected_ratio_unc, hist_1.divide(hist_2)[1])
@@ -82,7 +82,7 @@ def test_ratio_same_histogram(self):
         # duplicated in the ratio calculation (the first one is so to say not plotted)
         # Therefore, we also use duplicated bins here
         expected_ratio = np.ones(3)
-        expected_ratio_unc = np.array([0.81649658, 0.81649658, 1])
+        expected_ratio_unc = np.array([0.57735027, 0.57735027, 0.70710678])
 
         np.testing.assert_almost_equal(expected_ratio, hist_1.divide(hist_2)[0])
         np.testing.assert_almost_equal(expected_ratio_unc, hist_1.divide(hist_2)[1])

puma/tests/utils/test_histogram_utils.py

@@ -244,15 +244,7 @@ def setUp(self):
         self.denominator_unc = np.array([1, 0.3, 2, 1, 5, 3])
         self.step = np.array([1.6666667, 1.6666667, 0.5, 1, 0.7142857, 0.6, 0.1666667])
         self.step_unc = np.array(
-            [
-                0.580017,
-                0.580017,
-                0.1685312,
-                1.0111874,
-                0.1059653,
-                0.3041381,
-                0.0485913,
-            ]
+            [0.16666667, 0.16666667, 0.16666667, 0.15, 0.02857143, 0.05, 0.025]
         )
 
     def test_hist_ratio(self):
@@ -261,7 +253,6 @@ def test_hist_ratio(self):
             numerator=self.numerator,
             denominator=self.denominator,
             numerator_unc=self.numerator_unc,
-            denominator_unc=self.denominator_unc,
         )
 
         np.testing.assert_almost_equal(step, self.step)
@@ -274,7 +265,6 @@ def test_hist_not_same_length_numerator_denominator(self):
                 numerator=np.ones(2),
                 denominator=np.ones(3),
                 numerator_unc=np.ones(3),
-                denominator_unc=np.ones(3),
             )
 
     def test_hist_not_same_length_numerator_and_unc(self):
@@ -284,15 +274,4 @@ def test_hist_not_same_length_numerator_and_unc(self):
                 numerator=np.ones(3),
                 denominator=np.ones(3),
                 numerator_unc=np.ones(2),
-                denominator_unc=np.ones(3),
-            )
-
-    def test_hist_not_same_length_denomiantor_and_unc(self):
-        """Test case where denominator and uncertainty have not the same length."""
-        with self.assertRaises(AssertionError):
-            _, _ = hist_ratio(
-                numerator=np.ones(3),
-                denominator=np.ones(3),
-                numerator_unc=np.ones(3),
-                denominator_unc=np.ones(2),
             )

puma/utils/histogram.py

@@ -161,7 +161,6 @@ def hist_ratio(
     numerator,
     denominator,
     numerator_unc,
-    denominator_unc,
     step: bool = True,
 ):
     """
@@ -176,8 +175,6 @@ def hist_ratio(
         Denominator in the ratio calculation.
     numerator_unc : array_like
         Uncertainty of the numerator.
-    denominator_unc : array_like
-        Uncertainty of the denominator.
     step : bool
         if True duplicates first bin to match with step plotting function,
         by default True
@@ -196,40 +193,24 @@ def hist_ratio(
         If inputs don't have the same shape.
 
     """
-    numerator, denominator, numerator_unc, denominator_unc = (
+    numerator, denominator, numerator_unc = (
         np.array(numerator),
         np.array(denominator),
         np.array(numerator_unc),
-        np.array(denominator_unc),
     )
     if numerator.shape != denominator.shape:
         raise AssertionError("Numerator and denominator don't have the same legth")
     if numerator.shape != numerator_unc.shape:
         raise AssertionError("Numerator and numerator_unc don't have the same legth")
-    if denominator.shape != denominator_unc.shape:
-        raise (
-            AssertionError("Denominator and denominator_unc don't have the same legth")
-        )
     step_ratio = save_divide(numerator, denominator, 1 if step else np.inf)
 
-    # Calculate rel uncertainties
-    numerator_rel_unc = save_divide(
-        numerator_unc, numerator, default=0 if step else np.inf
-    )
-    denominator_rel_unc = save_divide(
-        denominator_unc, denominator, default=0 if step else np.inf
-    )
-
-    # Calculate rel uncertainty
-    step_rel_unc = np.sqrt(numerator_rel_unc**2 + denominator_rel_unc**2)
-
-    # Calculate final uncertainty
-    step_unc = step_ratio * step_rel_unc
+    # Calculate ratio uncertainty
+    step_unc = save_divide(numerator_unc, denominator, default=0 if step else np.inf)
 
     if step:
         # Add an extra bin in the beginning to have the same binning as the input
         # Otherwise, the ratio will not be exactly above each other (due to step)
         step_ratio = np.append(np.array([step_ratio[0]]), step_ratio)
-        step_unc = np.append(np.array([step_rel_unc[0]]), step_rel_unc) * step_ratio
+        step_unc = np.append(np.array([step_unc[0]]), step_unc)
 
     return step_ratio, step_unc

puma/var_vs_var.py

@@ -121,7 +121,6 @@ def divide(self, other, inverse: bool = False):
             numerator=denom if inverse else nom,
             denominator=nom if inverse else denom,
             numerator_unc=denom_err if inverse else nom_err,
-            denominator_unc=nom_err if inverse else denom_err,
             step=False,
         )
         return (ratio, ratio_err)