diff --git a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv index bcfc4cf1d..f110e64c2 100644 --- a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv +++ b/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 diff --git a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv index f494346c1..0b07d8199 100644 --- a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv +++ b/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 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 2160a3edd..50021a0bb 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/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: diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index f4868791b..c053ff8db 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/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: