Skip to content

IMAP-Lo systematic error calculation fixes #3315

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
vineetbansal wants to merge 1 commit into
base: dev
Choose a base branch
from
Open

IMAP-Lo systematic error calculation fixes #3315

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
28 changes: 14 additions & 14 deletions imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
71 changes: 44 additions & 27 deletions imap_processing/lo/l2/lo_l2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1009,16 +1009,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"]

# 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"]
dataset["ena_intensity_sys_err_minus"] = dataset["ena_intensity"] - (
dataset["ena_intensity"] * minus_multiplier
)

# Symmetric systematic error
dataset["ena_intensity_sys_err"] = np.sqrt(
dataset["ena_intensity_sys_err_minus"] * dataset["ena_intensity_sys_err_plus"]
)

return dataset
Expand Down Expand Up @@ -1048,16 +1056,24 @@ def calculate_backgrounds(dataset: xr.Dataset) -> xr.Dataset:
/ dataset["exposure_factor"] ** 2
)

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"]
)
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["bg_rate_sys_err_plus"] = (dataset["bg_rate"] * plus_multiplier) - dataset[
"bg_rate"
]

dataset["bg_rate_sys_err_minus"] = dataset["bg_rate"] - (
dataset["bg_rate"] * minus_multiplier
)

# 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"]
# Symmetric systematic error
dataset["bg_rate_sys_err"] = np.sqrt(
dataset["bg_rate_sys_err_minus"] * dataset["bg_rate_sys_err_plus"]
)

# Background intensity
Expand All @@ -1068,16 +1084,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
Expand Down
36 changes: 17 additions & 19 deletions imap_processing/tests/lo/test_lo_l2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1134,17 +1134,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):
Expand Down Expand Up @@ -1197,14 +1197,12 @@ 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"]
)
expected_sys_err = (
result["bg_rate"] * mean_gf_stat_uncert / dataset["geometric_factor"]
)
gf = dataset["geometric_factor"]
dg_minus = dataset["geometric_factor_stat_uncert_minus"]
dg_plus = dataset["geometric_factor_stat_uncert_plus"]
expected_sys_err_plus = result["bg_rate"] * dg_minus / (gf - dg_minus)
expected_sys_err_minus = result["bg_rate"] * dg_plus / (gf + dg_plus)
expected_sys_err = np.sqrt(expected_sys_err_minus * expected_sys_err_plus)
xr.testing.assert_allclose(result["bg_rate_sys_err"], expected_sys_err)

def test_calculate_backgrounds_zero_exposure(self):
Expand Down
Loading