diff --git a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v006.csv b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v006.csv new file mode 100644 index 0000000000..3fd3fae67c --- /dev/null +++ b/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, diff --git a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v006.csv b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v006.csv new file mode 100644 index 0000000000..9f2c8cff3d --- /dev/null +++ b/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 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 334b186463..56b5459d9e 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/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. diff --git a/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv index 0a6733a9e1..16d4faff19 100644 --- a/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv +++ b/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, diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 85a0c1fbca..e8068e37ec 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/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