Add bg_intensity_sys_err to Lo L2 maps

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

@@ -0,0 +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,Cntr_E_delta_minus,Cntr_E_delta_plus,GF_Trpl_H_unc
+# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV],[cm^2sr keV/keV]
+0,1,1,0.01633,0.00028,4.44705884E-05,4.02590505E-05,4.19612369E-05,0.00543,0.00543,3.11294119E-05
+0,2,2,0.03047,0.00043,5.01882355E-05,4.52880073E-05,4.72150361E-05,0.01002,0.01002,3.49411766E-05
+0,3,3,0.05576,0.00089,6.16235296E-05,5.58366786E-05,5.81934311E-05,0.01861,0.01861,4.32000002E-05
+0,4,4,0.10626,0.0019,7.11529414E-05,4.50746155E-05,4.88876658E-05,0.03331,0.03331,1.90588236E-05
+0,5,5,0.20004,0.0034,8.89411768E-05,5.85174388E-05,6.31197884E-05,0.06498,0.06498,2.85882354E-05
+0,6,6,0.40496,0.0073,1.12447059E-04,6.57915012E-05,7.22727021E-05,0.13164,0.13164,1.27058824E-05
+0,7,7,0.78729,0.022,1.42941177E-04,8.25394197E-05,9.08774549E-05,0.26235,0.26235,8.89411768E-06
+1,1,1,0.01719,0.00022,1.05E-04,8.82E-05,9.26E-05,0.00881,0.00881,
+1,2,2,0.03236,0.00036,1.22E-04,1.02E-04,1.07E-04,0.01604,0.01604,
+1,3,3,0.05948,0.00077,1.44E-04,1.21E-04,1.27E-04,0.02850,0.02850,
+1,4,4,0.11441,0.0013,1.73E-04,1.17E-04,1.26E-04,0.05313,0.05313,
+1,5,5,0.2137,0.003,2.15E-04,1.37E-04,1.49E-04,0.10560,0.10560,
+1,6,6,0.43736,0.0053,2.82E-04,1.65E-04,1.81E-04,0.21967,0.21967,
+1,7,7,0.83888,0.0117,3.61E-04,2.08E-04,2.29E-04,0.41360,0.41360,

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

@@ -0,0 +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,Cntr_E_delta_minus,Cntr_E_delta_plus
+# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV]
+0,1,1,0.01919,0.001372344105,1.98E-05,1.74E-05,1.81E-05,0.00582,0.00582
+0,2,2,0.03675,0.002537963082,3.18E-05,2.39E-05,2.53E-05,0.01110,0.01110
+0,3,3,0.07121,0.006591339559,5.34E-05,4.05E-05,4.29E-05,0.02178,0.02178
+0,4,4,0.14141,0.01743990907,8.64E-05,6.14E-05,6.56E-05,0.04147,0.04147
+0,5,5,0.274,0.03583900763,1.32E-04,9.07E-05,9.72E-05,0.08561,0.08561
+0,6,6,0.58503,0.09806681395,2.46E-04,1.53E-04,1.67E-04,0.18067,0.18067
+0,7,7,1.13506,0.2381076888,3.81E-04,2.23E-04,2.45E-04,0.36193,0.36193
+1,1,1,0.02043303552,0.001575924599,5.02E-05,4.41E-05,4.61E-05,0.00945,0.00945
+1,2,2,0.03972,0.002953020557,8.13E-05,6.09E-05,6.47E-05,0.01784,0.01784
+1,3,3,0.07648,0.007476494389,1.33E-04,1.01E-04,1.07E-04,0.03351,0.03351
+1,4,4,0.15353,0.01487953071,2.38E-04,1.69E-04,1.80E-04,0.06661,0.06661
+1,5,5,0.29846,0.04148950349,3.87E-04,2.67E-04,2.86E-04,0.13995,0.13995
+1,6,6,0.61524,0.08357629594,6.61E-04,4.11E-04,4.47E-04,0.30224,0.30224
+1,7,7,1.24282,0.2120806857,1.02E-03,6.10E-04,6.67E-04,0.56948,0.56948

imap_processing/lo/l2/lo_l2.py

@@ -69,6 +69,7 @@
     "bg_rate_stat_uncert",
     "bg_intensity",
     "bg_intensity_stat_uncert",
+    "bg_intensity_sys_err",
 )
 
 # The calibration ancillaries shipped with the package.
@@ -1115,6 +1116,10 @@ def _esa_calibration(species: str, esa_mode: int) -> xr.Dataset:
     therefore the downward excursion of the factor, and its ``_unc_minus`` the
     upward one.
 
+    Those bounds fold the recalibration of the factor in with its measurement
+    uncertainty. The measurement uncertainty alone, the symmetric ``_unc``
+    column, is what the background systematic error is derived from.
+
     Parameters
     ----------
     species : str
@@ -1142,6 +1147,9 @@ def _esa_calibration(species: str, esa_mode: int) -> xr.Dataset:
         - gf_data[f"{factor}_unc_plus"].to_numpy(),
         geometric_factor_high=geometric_factor
         + gf_data[f"{factor}_unc_minus"].to_numpy(),
+        geometric_factor_uncert=gf_data.get(
+            f"{factor}_unc", pd.Series(np.nan, index=gf_data.index)
+        ).to_numpy(),
     )
     data_vars = {
         name: xr.DataArray(arr, dims=[CoordNames.ENERGY_L2.value])
@@ -1776,6 +1784,9 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray:
     bg_rate_stat_uncert = np.sqrt(_divide(bg_rate, exposure))
     bg_intensity = _divide(bg_rate, geometric_factor * energy)
     bg_intensity_stat_uncert = _divide(bg_rate_stat_uncert, geometric_factor * energy)
+    gf_uncert = calibration["geometric_factor_uncert"]
+    gf_uncert_known = np.isfinite(gf_uncert)
+    bg_intensity_sys_err = bg_intensity * (gf_uncert / geometric_factor).fillna(0.0)
 
     # The Compton-Getting correction comes last, moving the intensities out of
     # the frame the instrument observed them in. The background is a spectrum
@@ -1811,9 +1822,12 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray:
             data=corrected_sys_err_minus
         )
 
-        bg_corrected, (bg_corrected_stat_uncert,) = _compton_getting_correct_intensity(
+        (
+            bg_corrected,
+            (bg_corrected_stat_uncert, bg_corrected_sys_err),
+        ) = _compton_getting_correct_intensity(
             bg_intensity.values,
-            (bg_intensity_stat_uncert.values,),
+            (bg_intensity_stat_uncert.values, bg_intensity_sys_err.values),
             cos_alpha.values,
             calibration,
             flux_corrector,
@@ -1822,6 +1836,7 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray:
         bg_intensity_stat_uncert = bg_intensity_stat_uncert.copy(
             data=bg_corrected_stat_uncert
         )
+        bg_intensity_sys_err = bg_intensity_sys_err.copy(data=bg_corrected_sys_err)
 
     variables = {
         "ena_count": counts,
@@ -1839,6 +1854,9 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray:
         "bg_rate_stat_uncert": bg_rate_stat_uncert,
         "bg_intensity": bg_intensity,
         "bg_intensity_stat_uncert": bg_intensity_stat_uncert,
+        "bg_intensity_sys_err": bg_intensity_sys_err.where(
+            gf_uncert_known, FILLVAL_FLOAT
+        ),
     }
 
     # A masked pixel has no ENA measurement to report.

imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv

@@ -1,19 +1,19 @@
-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,Cntr_E_delta_minus,Cntr_E_delta_plus
-# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV]
-0,1,1,0.01000,0.00010,1.00E-05,1.00E-06,1.00E-06,0.00100,0.00100
-0,2,1,0.01000,0.00010,9.99E-05,9.99E-06,9.99E-06,0.00100,0.00100
-0,2,2,0.02000,0.00020,2.00E-05,2.00E-06,2.00E-06,0.00200,0.00200
-0,3,2,0.02000,0.00020,9.99E-05,9.99E-06,9.99E-06,0.00200,0.00200
-0,3,3,0.04000,0.00040,3.00E-05,3.00E-06,3.00E-06,0.00400,0.00400
-0,4,4,0.08000,0.00080,4.00E-05,4.00E-06,4.00E-06,0.00800,0.00800
-0,5,5,0.16000,0.00160,5.00E-05,5.00E-06,5.00E-06,0.01600,0.01600
-0,6,6,0.32000,0.00320,6.00E-05,6.00E-06,6.00E-06,0.03200,0.03200
-0,7,7,0.64000,0.00640,7.00E-05,7.00E-06,7.00E-06,0.06400,0.06400
-1,1,1,0.01100,0.00011,1.10E-05,1.10E-06,1.10E-06,0.00110,0.00110
-1,2,1,0.01100,0.00011,9.99E-05,9.99E-06,9.99E-06,0.00110,0.00110
-1,2,2,0.02200,0.00022,2.20E-05,2.20E-06,2.20E-06,0.00220,0.00220
-1,3,3,0.04400,0.00044,3.30E-05,3.30E-06,3.30E-06,0.00440,0.00440
-1,4,4,0.08800,0.00088,4.40E-05,4.40E-06,4.40E-06,0.00880,0.00880
-1,5,5,0.17600,0.00176,5.50E-05,5.50E-06,5.50E-06,0.01760,0.01760
-1,6,6,0.35200,0.00352,6.60E-05,6.60E-06,6.60E-06,0.03520,0.03520
-1,7,7,0.70400,0.00704,7.70E-05,7.70E-06,7.70E-06,0.07040,0.07040
+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,Cntr_E_delta_minus,Cntr_E_delta_plus,GF_Trpl_H_unc
+# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV],[cm^2sr keV/keV]
+0,1,1,0.01000,0.00010,1.00E-05,1.00E-06,1.00E-06,0.00100,0.00100,5.00E-07
+0,2,1,0.01000,0.00010,9.99E-05,9.99E-06,9.99E-06,0.00100,0.00100,5.00E-06
+0,2,2,0.02000,0.00020,2.00E-05,2.00E-06,2.00E-06,0.00200,0.00200,1.00E-06
+0,3,2,0.02000,0.00020,9.99E-05,9.99E-06,9.99E-06,0.00200,0.00200,5.00E-06
+0,3,3,0.04000,0.00040,3.00E-05,3.00E-06,3.00E-06,0.00400,0.00400,1.50E-06
+0,4,4,0.08000,0.00080,4.00E-05,4.00E-06,4.00E-06,0.00800,0.00800,2.00E-06
+0,5,5,0.16000,0.00160,5.00E-05,5.00E-06,5.00E-06,0.01600,0.01600,2.50E-06
+0,6,6,0.32000,0.00320,6.00E-05,6.00E-06,6.00E-06,0.03200,0.03200,3.00E-06
+0,7,7,0.64000,0.00640,7.00E-05,7.00E-06,7.00E-06,0.06400,0.06400,3.50E-06
+1,1,1,0.01100,0.00011,1.10E-05,1.10E-06,1.10E-06,0.00110,0.00110,
+1,2,1,0.01100,0.00011,9.99E-05,9.99E-06,9.99E-06,0.00110,0.00110,
+1,2,2,0.02200,0.00022,2.20E-05,2.20E-06,2.20E-06,0.00220,0.00220,
+1,3,3,0.04400,0.00044,3.30E-05,3.30E-06,3.30E-06,0.00440,0.00440,
+1,4,4,0.08800,0.00088,4.40E-05,4.40E-06,4.40E-06,0.00880,0.00880,
+1,5,5,0.17600,0.00176,5.50E-05,5.50E-06,5.50E-06,0.01760,0.01760,
+1,6,6,0.35200,0.00352,6.60E-05,6.60E-06,6.60E-06,0.03520,0.03520,
+1,7,7,0.70400,0.00704,7.70E-05,7.70E-06,7.70E-06,0.07040,0.07040,

imap_processing/tests/lo/test_lo_l2.py

@@ -114,6 +114,10 @@
     1: np.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7]) * 1e-5,
 }
 
+# The symmetric geometric factor uncertainty of the same file, as a fraction of
+# each factor. Only the HiRes levels have one, as in the flight ancillary.
+GEO_FACTOR_UNCERT_FRACTION = 0.05
+
 # The ESA passband half-widths of the same file, a tenth of each center energy.
 ESA_ENERGY_DELTAS = {mode: energies / 10 for mode, energies in ESA_ENERGIES.items()}
 
@@ -605,6 +609,9 @@ def test_the_map_is_binned_in_its_esa_mode(
             )
 
         np.testing.assert_allclose(dataset["energy"].values, ESA_ENERGIES[esa_mode])
+        # Only HiRes has a geometric factor uncertainty to derive it from.
+        sys_err = dataset["bg_intensity_sys_err"].values
+        assert np.all(is_fill(sys_err)) == (esa_mode == 1)
 
 
 class TestMapStructure:
@@ -624,6 +631,7 @@ class TestMapStructure:
         "bg_rate_stat_uncert",
         "bg_intensity",
         "bg_intensity_stat_uncert",
+        "bg_intensity_sys_err",
     )
 
     def test_map_dims_and_source(self, full_map):
@@ -797,6 +805,16 @@ def test_background_rate(self, full_map):
             )
             assert np.all(is_fill(bg_rate[~exposed]))
 
+            # The systematic error is the fractional uncertainty of the factor.
+            bg_intensity = dataset["bg_intensity"].values[0, energy_index]
+            sys_err = dataset["bg_intensity_sys_err"].values[0, energy_index]
+            np.testing.assert_allclose(
+                sys_err[exposed],
+                bg_intensity[exposed] * GEO_FACTOR_UNCERT_FRACTION,
+                rtol=1e-5,
+            )
+            assert np.all(is_fill(sys_err[~exposed]))
+
     def test_systematic_error_bounds(self, full_map):
         """The systematic error is bracketed by the G-factor excursions."""
         dataset, _ = full_map
@@ -1171,6 +1189,14 @@ def test_intensities_are_corrected(self, cg_maps):
                 ~np.isclose(corrected[variable].values[lit], raw[variable].values[lit])
             ), f"{variable} was not corrected"
 
+        # The background systematic error moves with the background.
+        exposed = corrected["exposure_factor"].values > 0
+        np.testing.assert_allclose(
+            corrected["bg_intensity_sys_err"].values[exposed],
+            corrected["bg_intensity"].values[exposed] * GEO_FACTOR_UNCERT_FRACTION,
+            rtol=1e-5,
+        )
+
     def test_the_correction_leaks_no_nan(self, cg_maps):
         """A pixel the correction says nothing about comes out at zero, not NaN."""
         corrected, _ = cg_maps
@@ -1182,6 +1208,7 @@ def test_the_correction_leaks_no_nan(self, cg_maps):
             "ena_intensity_sys_err_minus",
             "bg_intensity",
             "bg_intensity_stat_uncert",
+            "bg_intensity_sys_err",
         ):
             values = corrected[variable].values
             assert np.isfinite(values).all(), f"{variable} holds NaN or inf"
@@ -1570,6 +1597,7 @@ def test_the_exposure_and_background_are_left_alone(self, masked_maps):
             "bg_rate_stat_uncert",
             "bg_intensity",
             "bg_intensity_stat_uncert",
+            "bg_intensity_sys_err",
         ):
             np.testing.assert_allclose(
                 masked[name].values, raw[name].values, rtol=1e-6, err_msg=name