Skip to content

2435 hi l2 updated systematic error calculation #2560

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
subagonsouth wants to merge 6 commits into
base: dev
Choose a base branch
from
Open

2435 hi l2 updated systematic error calculation #2560

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
6891852
Update equation for ena_intensity_sys_err
subagonsouth Nov 18, 2025
3e26124
Fix systematic error calculation - multiply instead of divide
subagonsouth Dec 4, 2025
0baf9d6
Fix systematic error calculation - use bg_rates, not bg_rates_unc
subagonsouth Dec 8, 2025
1507a68
Remove systematic error from weights used to combine calibration prod…
subagonsouth Dec 10, 2025
84537f0
Update test to remove use of statistical uncertainty when combinging …
subagonsouth Dec 10, 2025
1f9faf4
Get rid of check that variance is >= 1 in calculation of imporved var…
subagonsouth Dec 15, 2025
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
18 changes: 6 additions & 12 deletions imap_processing/hi/hi_l2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,11 @@ def calculate_ena_intensity(
map_ds["ena_intensity_stat_uncert"] = (
map_ds["ena_signal_rate_stat_unc"] / flux_conversion_divisor
)
map_ds["ena_intensity_sys_err"] = map_ds["bg_rates_unc"] / flux_conversion_divisor
map_ds["ena_intensity_sys_err"] = (
np.sqrt(map_ds["bg_rates"] * map_ds["exposure_factor"])
/ map_ds["exposure_factor"]
/ flux_conversion_divisor
)

# Combine calibration products using proper weighted averaging
# as described in Hi Algorithm Document Section 3.1.2
Expand Down Expand Up @@ -367,16 +371,11 @@ def combine_calibration_products(
map_ds, geometric_factors, esa_energies
)

# Calculate total variance
# Note that sys_err contains uncertainty, so it must be squared to get
# the systematic variance needed in this equation.
total_variance = improved_stat_variance + sys_err**2

# Perform inverse-variance weighted averaging
# Handle divide by zero and invalid values
with np.errstate(divide="ignore", invalid="ignore"):
# Use total variance weights for flux combination
flux_weights = 1.0 / total_variance
flux_weights = 1.0 / improved_stat_variance
weighted_flux_sum = (ena_flux * flux_weights).sum(dim="calibration_prod")
combined_flux = weighted_flux_sum / flux_weights.sum(dim="calibration_prod")

Expand Down Expand Up @@ -450,11 +449,6 @@ def _calculate_improved_stat_variance(
# Total count rates for Poisson uncertainty calculation
total_count_rates_for_uncertainty = map_ds["bg_rates"] + averaged_signal_rates

# Ensure non-negative values for sqrt and minimum of 1 for uncertainty calculation
total_count_rates_for_uncertainty = xr.where(
total_count_rates_for_uncertainty < 1, 1, total_count_rates_for_uncertainty
)

logger.debug("Computing improved flux uncertainties")
# Statistical variance:
with np.errstate(divide="ignore", invalid="ignore"):
Expand Down
2 changes: 1 addition & 1 deletion imap_processing/tests/hi/test_hi_l2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ def test_statistical_uncertainty_combination_correctness():
# Manual calculation of expected statistical uncertainty combination
# combined_stat_unc = sqrt(1/sum(1 / stat_unc**2))
expected_combined_stat_unc = np.sqrt(1 / np.sum(1 / stat_unc_values**2))
flux_weights = 1.0 / (np.array([101, 101]) + np.array([4, 16]))
flux_weights = 1.0 / np.array([101, 101])
expected_flux = np.sum(flux_values.squeeze() * flux_weights) / np.sum(flux_weights)

np.testing.assert_almost_equal(
Expand Down