IMAP-Lo systematic error calculation fixes (#3315)

* fixes for issue 3314

* Derive Lo bg_rate_sys_err from the background, not the geometric factor

The algorithm document (Appendix A, eq 64) specifies the background rate
systematic as the exposure time weighted average of the per-pset systematic
error.

The systematic is a single symmetric value, so _plus and _minus are set equal
to it for now.

Note that the systematic errors are all 0s, since no constant has been
introduced yet to populate this (to be Imap-Lo team supplied) value.

Author: vineetbansal

imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv

@@ -1,16 +1,16 @@
 esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_H,GF_Trpl_H_unc_minus,GF_Trpl_H_unc_plus
 # [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV]
-0,1,1,0.01633,0.28,4.45E-05,4.03E-05,4.20E-05
-0,2,2,0.03047,0.43,5.02E-05,4.53E-05,4.72E-05
-0,3,3,0.05576,0.89,6.16E-05,5.58E-05,5.82E-05
-0,4,4,0.10626,1.9,7.12E-05,4.51E-05,4.89E-05
-0,5,5,0.20004,3.4,8.89E-05,5.85E-05,6.31E-05
-0,6,6,0.40496,7.3,1.12E-04,6.58E-05,7.23E-05
-0,7,7,0.78729,22,1.43E-04,8.25E-05,9.09E-05
-1,1,1,0.01719,0.22,1.05E-04,8.82E-05,9.26E-05
-1,2,2,0.03236,0.36,1.22E-04,1.02E-04,1.07E-04
-1,3,3,0.05948,0.77,1.44E-04,1.21E-04,1.27E-04
-1,4,4,0.11441,1.3,1.73E-04,1.17E-04,1.26E-04
-1,5,5,0.2137,3,2.15E-04,1.37E-04,1.49E-04
-1,6,6,0.43736,5.3,2.82E-04,1.65E-04,1.81E-04
-1,7,7,0.83888,11.7,3.61E-04,2.08E-04,2.29E-04
+0,1,1,0.01633,0.00028,4.45E-05,4.03E-05,4.20E-05
+0,2,2,0.03047,0.00043,5.02E-05,4.53E-05,4.72E-05
+0,3,3,0.05576,0.00089,6.16E-05,5.58E-05,5.82E-05
+0,4,4,0.10626,0.0019,7.12E-05,4.51E-05,4.89E-05
+0,5,5,0.20004,0.0034,8.89E-05,5.85E-05,6.31E-05
+0,6,6,0.40496,0.0073,1.12E-04,6.58E-05,7.23E-05
+0,7,7,0.78729,0.022,1.43E-04,8.25E-05,9.09E-05
+1,1,1,0.01719,0.00022,1.05E-04,8.82E-05,9.26E-05
+1,2,2,0.03236,0.00036,1.22E-04,1.02E-04,1.07E-04
+1,3,3,0.05948,0.00077,1.44E-04,1.21E-04,1.27E-04
+1,4,4,0.11441,0.0013,1.73E-04,1.17E-04,1.26E-04
+1,5,5,0.2137,0.003,2.15E-04,1.37E-04,1.49E-04
+1,6,6,0.43736,0.0053,2.82E-04,1.65E-04,1.81E-04
+1,7,7,0.83888,0.0117,3.61E-04,2.08E-04,2.29E-04

imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv

@@ -1,16 +1,16 @@
 esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_O,GF_Trpl_O_unc_minus,GF_Trpl_O_unc_plus
 # [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV]
-0,1,1,0.01919,1.372344105,1.98E-05,1.74E-05,1.81E-05
-0,2,2,0.03675,2.537963082,3.18E-05,2.39E-05,2.53E-05
-0,3,3,0.07121,6.591339559,5.34E-05,4.05E-05,4.29E-05
-0,4,4,0.14141,17.43990907,8.64E-05,6.14E-05,6.56E-05
-0,5,5,0.274,35.83900763,1.32E-04,9.07E-05,9.72E-05
-0,6,6,0.58503,98.06681395,2.46E-04,1.53E-04,1.67E-04
-0,7,7,1.13506,238.1076888,3.81E-04,2.23E-04,2.45E-04
-1,1,1,0.02043303552,1.575924599,5.02E-05,4.41E-05,4.61E-05
-1,2,2,0.03972,2.953020557,8.13E-05,6.09E-05,6.47E-05
-1,3,3,0.07648,7.476494389,1.33E-04,1.01E-04,1.07E-04
-1,4,4,0.15353,14.87953071,2.38E-04,1.69E-04,1.80E-04
-1,5,5,0.29846,41.48950349,3.87E-04,2.67E-04,2.86E-04
-1,6,6,0.61524,83.57629594,6.61E-04,4.11E-04,4.47E-04
-1,7,7,1.24282,212.0806857,1.02E-03,6.10E-04,6.67E-04
+0,1,1,0.01919,0.001372344105,1.98E-05,1.74E-05,1.81E-05
+0,2,2,0.03675,0.002537963082,3.18E-05,2.39E-05,2.53E-05
+0,3,3,0.07121,0.006591339559,5.34E-05,4.05E-05,4.29E-05
+0,4,4,0.14141,0.01743990907,8.64E-05,6.14E-05,6.56E-05
+0,5,5,0.274,0.03583900763,1.32E-04,9.07E-05,9.72E-05
+0,6,6,0.58503,0.09806681395,2.46E-04,1.53E-04,1.67E-04
+0,7,7,1.13506,0.2381076888,3.81E-04,2.23E-04,2.45E-04
+1,1,1,0.02043303552,0.001575924599,5.02E-05,4.41E-05,4.61E-05
+1,2,2,0.03972,0.002953020557,8.13E-05,6.09E-05,6.47E-05
+1,3,3,0.07648,0.007476494389,1.33E-04,1.01E-04,1.07E-04
+1,4,4,0.15353,0.01487953071,2.38E-04,1.69E-04,1.80E-04
+1,5,5,0.29846,0.04148950349,3.87E-04,2.67E-04,2.86E-04
+1,6,6,0.61524,0.08357629594,6.61E-04,4.11E-04,4.47E-04
+1,7,7,1.24282,0.2120806857,1.02E-03,6.10E-04,6.67E-04

imap_processing/lo/l2/lo_l2.py

@@ -499,6 +499,7 @@ def normalize_pset_coordinates(pset: xr.Dataset, species: str) -> xr.Dataset:
         f"{species}_counts": "counts",
         f"{species}_background_rates": "bg_rate",
         f"{species}_background_rates_stat_uncert": "bg_rate_stat_uncert",
+        f"{species}_background_rates_sys_err": "bg_rate_sys_err",
     }
     pset_renamed = pset_renamed.rename_vars(rename_map)
 
@@ -590,6 +591,10 @@ def calculate_efficiency_corrected_quantities(
     pset["bg_rate_stat_uncert_exposure_factor2"] = (
         pset["bg_rate_stat_uncert"] ** 2 * pset["exposure_factor"] ** 2
     )
+    # background systematic * exposure_factor for weighted average.
+    pset["bg_rate_sys_err_exposure_factor"] = (
+        pset["bg_rate_sys_err"] * pset["exposure_factor"]
+    )
 
     return pset
 
@@ -629,8 +634,10 @@ def project_pset_to_map(
         "counts_over_eff_squared",
         "bg_rate",
         "bg_rate_stat_uncert",
+        "bg_rate_sys_err",
         "bg_rate_exposure_factor",
         "bg_rate_stat_uncert_exposure_factor2",
+        "bg_rate_sys_err_exposure_factor",
     ]
     if cg_correct:
         value_keys.append("energy_sc_exposure_factor")
@@ -1009,16 +1016,24 @@ def calculate_intensities(dataset: xr.Dataset) -> xr.Dataset:
         dataset["counts_over_eff_squared"]
     ) / (dataset["geometric_factor"] * dataset["energy"] * dataset["exposure_factor"])
 
-    for suffix in ("minus", "plus"):
-        dataset[f"ena_intensity_sys_err_{suffix}"] = (
-            dataset["ena_intensity"]
-            * dataset[f"geometric_factor_stat_uncert_{suffix}"]
-            / dataset["geometric_factor"]
-        )
+    plus_multiplier = dataset["geometric_factor"] / (
+        dataset["geometric_factor"] - dataset["geometric_factor_stat_uncert_minus"]
+    )
+    minus_multiplier = dataset["geometric_factor"] / (
+        dataset["geometric_factor"] + dataset["geometric_factor_stat_uncert_plus"]
+    )
+
+    dataset["ena_intensity_sys_err_plus"] = (
+        dataset["ena_intensity"] * plus_multiplier
+    ) - dataset["ena_intensity"]
+
+    dataset["ena_intensity_sys_err_minus"] = dataset["ena_intensity"] - (
+        dataset["ena_intensity"] * minus_multiplier
+    )
 
-    # Symmetric systematic error (mean of the asymmetric minus/plus bounds)
-    dataset["ena_intensity_sys_err"] = 0.5 * (
-        dataset["ena_intensity_sys_err_minus"] + dataset["ena_intensity_sys_err_plus"]
+    # Symmetric systematic error
+    dataset["ena_intensity_sys_err"] = np.sqrt(
+        dataset["ena_intensity_sys_err_minus"] * dataset["ena_intensity_sys_err_plus"]
     )
 
     return dataset
@@ -1039,25 +1054,27 @@ def calculate_backgrounds(dataset: xr.Dataset) -> xr.Dataset:
         Dataset with calculated background rates and intensities for the
         specified species.
     """
-    # Equation 6 from mapping document (background rate)
+    # Equation 62 from mapping document (background rate)
     # exposure time weighted average of the background rates
     dataset["bg_rate"] = dataset["bg_rate_exposure_factor"] / dataset["exposure_factor"]
-    # Equation 7 from mapping document (background statistical uncertainty)
+    # Equation 63 from mapping document (background statistical uncertainty)
     dataset["bg_rate_stat_uncert"] = np.sqrt(
         dataset["bg_rate_stat_uncert_exposure_factor2"]
         / dataset["exposure_factor"] ** 2
     )
+    # Equation 64 from mapping document (background systematic error).
+    dataset["bg_rate_sys_err"] = (
+        dataset["bg_rate_sys_err_exposure_factor"] / dataset["exposure_factor"]
+    )
+    # The background rate systematic is a single symmetric value.
+    dataset["bg_rate_sys_err_plus"] = dataset["bg_rate_sys_err"].copy()
+    dataset["bg_rate_sys_err_minus"] = dataset["bg_rate_sys_err"].copy()
 
-    for suffix in ("minus", "plus"):
-        dataset[f"bg_rate_sys_err_{suffix}"] = (
-            dataset["bg_rate"]
-            * dataset[f"geometric_factor_stat_uncert_{suffix}"]
-            / dataset["geometric_factor"]
-        )
-
-    # Symmetric systematic error (mean of the asymmetric minus/plus bounds)
-    dataset["bg_rate_sys_err"] = 0.5 * (
-        dataset["bg_rate_sys_err_minus"] + dataset["bg_rate_sys_err_plus"]
+    plus_multiplier = dataset["geometric_factor"] / (
+        dataset["geometric_factor"] - dataset["geometric_factor_stat_uncert_minus"]
+    )
+    minus_multiplier = dataset["geometric_factor"] / (
+        dataset["geometric_factor"] + dataset["geometric_factor_stat_uncert_plus"]
     )
 
     # Background intensity
@@ -1068,16 +1085,17 @@ def calculate_backgrounds(dataset: xr.Dataset) -> xr.Dataset:
         dataset["geometric_factor"] * dataset["energy"]
     )
 
-    for suffix in ("minus", "plus"):
-        dataset[f"bg_intensity_sys_err_{suffix}"] = (
-            dataset["bg_intensity"]
-            * dataset[f"geometric_factor_stat_uncert_{suffix}"]
-            / dataset["geometric_factor"]
-        )
+    dataset["bg_intensity_sys_err_plus"] = (
+        dataset["bg_intensity"] * plus_multiplier
+    ) - dataset["bg_intensity"]
+
+    dataset["bg_intensity_sys_err_minus"] = dataset["bg_intensity"] - (
+        dataset["bg_intensity"] * minus_multiplier
+    )
 
-    # Symmetric systematic error (mean of the asymmetric minus/plus bounds)
-    dataset["bg_intensity_sys_err"] = 0.5 * (
-        dataset["bg_intensity_sys_err_minus"] + dataset["bg_intensity_sys_err_plus"]
+    # Symmetric systematic error
+    dataset["bg_intensity_sys_err"] = np.sqrt(
+        dataset["bg_intensity_sys_err_minus"] * dataset["bg_intensity_sys_err_plus"]
     )
 
     return dataset
@@ -1427,6 +1445,7 @@ def cleanup_intermediate_variables(dataset: xr.Dataset) -> xr.Dataset:
         "counts_over_eff_squared",
         "bg_rate_exposure_factor",
         "bg_rate_stat_uncert_exposure_factor2",
+        "bg_rate_sys_err_exposure_factor",
     ]
 
     for potential_var in potential_vars:

imap_processing/tests/lo/test_lo_l2.py

@@ -390,6 +390,13 @@ def sample_dataset_with_background_intermediates():
         bg_rate_stat_uncert_exposure_factor2,
     )
 
+    # Background systematic error times exposure time
+    bg_rate_sys_err_exposure_factor = np.ones((1, n_energy)) * 0.05
+    dataset["bg_rate_sys_err_exposure_factor"] = (
+        ("epoch", "energy"),
+        bg_rate_sys_err_exposure_factor,
+    )
+
     # Add exposure time (using current naming convention)
     exposure = np.ones((1, n_energy)) * 1.0  # 1 second
     dataset["exposure_factor"] = (("epoch", "energy"), exposure)
@@ -848,6 +855,10 @@ def test_normalize_coordinates_basic(self, species):
                     PSET_DIMS,
                     np.ones(PSET_SHAPE) * 0.01,
                 ),
+                f"{species}_background_rates_sys_err": (
+                    PSET_DIMS,
+                    np.ones(PSET_SHAPE) * 0.02,
+                ),
             },
             coords={
                 "epoch": [8.1794907049e17],
@@ -880,6 +891,7 @@ def test_normalize_coordinates_basic(self, species):
         assert "exposure_factor" in result.data_vars
         assert "bg_rate" in result.data_vars
         assert "bg_rate_stat_uncert" in result.data_vars
+        assert "bg_rate_sys_err" in result.data_vars
 
         # Check that old variable names are gone
         assert f"{species}_counts" not in result.data_vars
@@ -919,6 +931,10 @@ def test_normalize_coordinates_removes_old_coordinate(self):
                     PSET_DIMS,
                     np.ones(PSET_SHAPE) * 0.01,
                 ),
+                f"{species}_background_rates_sys_err": (
+                    PSET_DIMS,
+                    np.ones(PSET_SHAPE) * 0.02,
+                ),
                 "esa_energy_step_var": xr.DataArray([1, 2, 3, 4, 5, 6, 7]),  # Variable
             },
             coords={
@@ -1010,6 +1026,7 @@ def test_calculate_efficiency_corrected_quantities(self):
                 "exposure_factor": (("energy",), np.ones(7) * 1.0),  # 1 second
                 "bg_rate": (("energy",), np.ones(7) * 0.1),  # 0.1 counts/s
                 "bg_rate_stat_uncert": (("energy",), np.ones(7) * 0.01),  # uncertainty
+                "bg_rate_sys_err": (("energy",), np.ones(7) * 0.02),  # systematic
                 "efficiency": (
                     ("energy",),
                     np.array([0.8, 0.85, 0.9, 0.95, 0.88, 0.92, 0.87]),
@@ -1025,6 +1042,7 @@ def test_calculate_efficiency_corrected_quantities(self):
         assert "counts_over_eff_squared" in result.data_vars
         assert "bg_rate_exposure_factor" in result.data_vars
         assert "bg_rate_stat_uncert_exposure_factor2" in result.data_vars
+        assert "bg_rate_sys_err_exposure_factor" in result.data_vars
 
         # Check dimensions
         assert result["counts_over_eff"].dims == pset["counts"].dims
@@ -1052,6 +1070,11 @@ def test_calculate_efficiency_corrected_quantities(self):
             result["bg_rate_stat_uncert_exposure_factor2"], expected_bg_uncert_exposure
         )
 
+        expected_bg_sys_err_exposure = pset["bg_rate_sys_err"] * pset["exposure_factor"]
+        xr.testing.assert_allclose(
+            result["bg_rate_sys_err_exposure_factor"], expected_bg_sys_err_exposure
+        )
+
 
 class TestCalculateRates:
     """Tests for the calculate_rates function."""
@@ -1134,17 +1157,17 @@ def test_calculate_intensities_basic(self, sample_dataset_with_geometric_factors
             result["ena_intensity_stat_uncert"], expected_stat_uncert
         )
 
-        # Check systematic uncertainty calculation. The single `_sys_err` is the
-        # mean of the asymmetric minus/plus bounds.
-        mean_gf_stat_uncert = 0.5 * (
-            sample_dataset_with_geometric_factors["geometric_factor_stat_uncert_minus"]
-            + sample_dataset_with_geometric_factors["geometric_factor_stat_uncert_plus"]
-        )
-        expected_sys_err = (
-            result["ena_intensity"]
-            * mean_gf_stat_uncert
-            / sample_dataset_with_geometric_factors["geometric_factor"]
-        )
+        # Check systematic uncertainty calculation
+        gf = sample_dataset_with_geometric_factors["geometric_factor"]
+        dg_minus = sample_dataset_with_geometric_factors[
+            "geometric_factor_stat_uncert_minus"
+        ]
+        dg_plus = sample_dataset_with_geometric_factors[
+            "geometric_factor_stat_uncert_plus"
+        ]
+        expected_sys_err_plus = result["ena_intensity"] * dg_minus / (gf - dg_minus)
+        expected_sys_err_minus = result["ena_intensity"] * dg_plus / (gf + dg_plus)
+        expected_sys_err = np.sqrt(expected_sys_err_minus * expected_sys_err_plus)
         xr.testing.assert_allclose(result["ena_intensity_sys_err"], expected_sys_err)
 
     def test_calculate_intensities_missing_variables(self):
@@ -1196,14 +1219,9 @@ def test_calculate_backgrounds_basic(
         )
         xr.testing.assert_allclose(result["bg_rate_stat_uncert"], expected_stat_uncert)
 
-        # Check systematic uncertainty calculation
-        # (mean(geometric_factor_stat_uncert bounds) / geometric_factor) * bg_rate
-        mean_gf_stat_uncert = 0.5 * (
-            dataset["geometric_factor_stat_uncert_minus"]
-            + dataset["geometric_factor_stat_uncert_plus"]
-        )
+        # Check systematic error calculation
         expected_sys_err = (
-            result["bg_rate"] * mean_gf_stat_uncert / dataset["geometric_factor"]
+            dataset["bg_rate_sys_err_exposure_factor"] / dataset["exposure_factor"]
         )
         xr.testing.assert_allclose(result["bg_rate_sys_err"], expected_sys_err)
 
@@ -1219,6 +1237,10 @@ def test_calculate_backgrounds_zero_exposure(self):
                     ("epoch", "energy"),
                     np.ones((1, 7)) * 0.004,
                 ),
+                "bg_rate_sys_err_exposure_factor": (
+                    ("epoch", "energy"),
+                    np.ones((1, 7)) * 0.05,
+                ),
                 "exposure_factor": (("epoch", "energy"), np.zeros((1, 7))),
                 "geometric_factor": (("energy",), np.ones(7) * 1e-4),
                 "geometric_factor_stat_uncert_minus": (("energy",), np.ones(7) * 1e-5),
@@ -1839,6 +1861,7 @@ def test_cleanup_intermediate_variables(self):
                 "counts_over_eff_squared": (("energy",), np.ones(7)),
                 "bg_rate_exposure_factor": (("energy",), np.ones(7)),
                 "bg_rate_stat_uncert_exposure_factor2": (("energy",), np.ones(7)),
+                "bg_rate_sys_err_exposure_factor": (("energy",), np.ones(7)),
                 "ena_intensity": (("energy",), np.ones(7)),  # Should be kept
                 "exposure_factor": (("energy",), np.ones(7)),  # Should be kept
             }
@@ -1856,6 +1879,7 @@ def test_cleanup_intermediate_variables(self):
         assert "counts_over_eff_squared" not in result.data_vars
         assert "bg_rate_exposure_factor" not in result.data_vars
         assert "bg_rate_stat_uncert_exposure_factor2" not in result.data_vars
+        assert "bg_rate_sys_err_exposure_factor" not in result.data_vars
 
     def test_cleanup_partial_variables(self):
         """Test cleanup when only some intermediate variables exist."""
@@ -2327,6 +2351,10 @@ def test_calculate_all_rates_and_intensities_complete(self):
                     ("energy",),
                     np.ones(7) * 0.009,
                 ),
+                "bg_rate_sys_err_exposure_factor": (
+                    ("energy",),
+                    np.ones(7) * 0.06,
+                ),
             }
         )
 
@@ -2373,6 +2401,10 @@ def test_calculate_all_rates_with_cg_correction(
             ("epoch", "energy"),
             np.ones((1, 7)) * 0.009,
         )
+        dataset["bg_rate_sys_err_exposure_factor"] = (
+            ("epoch", "energy"),
+            np.ones((1, 7)) * 0.06,
+        )
         dataset["energy_sc_exposure_factor"] = xr.ones_like(dataset["ena_intensity"])
 
         # Mock the interpolation function
@@ -2425,6 +2457,10 @@ def test_calculate_all_rates_cg_with_other_corrections(
             ("epoch", "energy"),
             np.ones((1, 7)) * 0.009,
         )
+        dataset["bg_rate_sys_err_exposure_factor"] = (
+            ("epoch", "energy"),
+            np.ones((1, 7)) * 0.06,
+        )
         dataset["energy_sc_exposure_factor"] = xr.ones_like(dataset["ena_intensity"])
 
         with patch(
@@ -2869,8 +2905,10 @@ def test_project_pset_to_map_value_keys(self, minimal_pset_for_species):
             "counts_over_eff_squared",
             "bg_rate",
             "bg_rate_stat_uncert",
+            "bg_rate_sys_err",
             "bg_rate_exposure_factor",
             "bg_rate_stat_uncert_exposure_factor2",
+            "bg_rate_sys_err_exposure_factor",
         ]
 
         for key in expected_keys: