diff --git a/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml b/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml
index ec6ef36a2..22835efc8 100644
--- a/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml
+++ b/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml
@@ -183,6 +183,16 @@ ena_intensity_sys_err_plus:
   UNITS: cm -2 s -1 sr -1 keV -1
   VAR_TYPE: support_data
 
+ena_intensity_calibration_sys_err:
+  <<: *default_float32
+  DEPEND_0: epoch
+  DICT_KEY: SPASE>Particle>ParticleType:Atom,ParticleQuantity:NumberFlux,Qualifier:Uncertainty,CoordinateSystemName:HAE,CoordinateRepresentation:Spherical
+  DISPLAY_TYPE: map_image
+  FIELDNAM: Intensity calibration systematic error
+  LABLAXIS: Calibration Systematic Error
+  UNITS: cm -2 s -1 sr -1 keV -1
+  VAR_TYPE: support_data
+
 ena_count_rate:
   <<: *default_float32
   DEPEND_0: epoch
diff --git a/imap_processing/cdf/config/imap_enamaps_l2-rectangular_variable_attrs.yaml b/imap_processing/cdf/config/imap_enamaps_l2-rectangular_variable_attrs.yaml
index 64367758a..ec3b639fe 100644
--- a/imap_processing/cdf/config/imap_enamaps_l2-rectangular_variable_attrs.yaml
+++ b/imap_processing/cdf/config/imap_enamaps_l2-rectangular_variable_attrs.yaml
@@ -115,6 +115,14 @@ ena_intensity_sys_err_plus:
   LABL_PTR_2: longitude_label
   LABL_PTR_3: latitude_label
 
+ena_intensity_calibration_sys_err:
+  DEPEND_1: energy
+  DEPEND_2: longitude
+  DEPEND_3: latitude
+  LABL_PTR_1: energy_label
+  LABL_PTR_2: longitude_label
+  LABL_PTR_3: latitude_label
+
 ena_count_rate:
   DEPEND_1: energy
   DEPEND_2: longitude
diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py
index 73e72429e..ee8270695 100644
--- a/imap_processing/ena_maps/ena_maps.py
+++ b/imap_processing/ena_maps/ena_maps.py
@@ -756,6 +756,7 @@ class HiPointingSet(LoHiBasePointingSet):
         "exposure_times": "exposure_factor",
         "background_rates": "bg_rate",
         "background_rates_uncertainty": "bg_rate_sys_err",
+        "counts": "ena_count",
     }
 
     def __init__(self, dataset: xr.Dataset | str | Path):
diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py
index 191239d2c..95c078543 100644
--- a/imap_processing/ena_maps/utils/naming.py
+++ b/imap_processing/ena_maps/utils/naming.py
@@ -195,6 +195,7 @@ def build_map_var_catdesc(self, support_var_name: str) -> str | None:
             "ena_intensity_sys_err": "Inten Sys. Err.",
             "ena_intensity_sys_err_minus": "Inten Sys. Err. Lower",
             "ena_intensity_sys_err_plus": "Inten Sys. Err. Upper",
+            "ena_intensity_calibration_sys_err": "Inten Calibration Sys. Err.",
             "ena_spectral_index": "Spectral Index",
             "ena_spectral_index_stat_uncert": "Spectral Stat. Unc.",
             "ena_spectral_scalar": "Spectral Scalar",
diff --git a/imap_processing/hi/hi_l2.py b/imap_processing/hi/hi_l2.py
index c0a31208a..5b3bcc62c 100644
--- a/imap_processing/hi/hi_l2.py
+++ b/imap_processing/hi/hi_l2.py
@@ -26,7 +26,7 @@
 logger = logging.getLogger(__name__)
 
 SC_FRAME_VARS_TO_PROJECT = {
-    "counts",
+    "ena_count",
     "exposure_factor",
     "bg_rate",
     "bg_rate_sys_err",
@@ -322,7 +322,7 @@ def calculate_all_rates_and_intensities(
     Parameters
     ----------
     map_ds : xarray.Dataset
-        Map dataset with projected PSET data (counts, exposure_factor, bg_rate,
+        Map dataset with projected PSET data (ena_count, exposure_factor, bg_rate,
         energy_delta_minus, energy_delta_plus, etc.) and an `energy` coordinate
         containing the ESA nominal central energies in keV.
     l2_ancillary_path_dict : dict[str, pathlib.Path]
@@ -344,7 +344,21 @@ def calculate_all_rates_and_intensities(
     logger.debug("Calculating ENA intensities")
     map_ds = calculate_ena_intensity(map_ds, l2_ancillary_path_dict, descriptor)
 
-    # Step 3: Handle obs_date variable type conversion
+    # Step 3: Add calibration systematic uncertainty in quadrature with the
+    # background-associated systematic. This is a percentage of the intensity.
+    # calib_sys_err is computed from ena_intensity after flux correction has
+    # already been applied (see calculate_ena_intensity), so it is automatically
+    # consistent with the corrected intensity and needs no separate scaling.
+    # This must happen before the CG interpolation step below, which requires
+    # ena_intensity_sys_err to already exist (even though it deliberately
+    # leaves it unmodified -- see update_sys_err=False).
+    logger.debug("Adding calibration systematic uncertainty")
+    bg_sys_err = map_ds["bg_intensity_sys_err"]
+    calib_sys_err = CALIBRATION_UNCERTAINTY_FRACTION * map_ds["ena_intensity"]
+    map_ds["ena_intensity_calibration_sys_err"] = calib_sys_err
+    map_ds["ena_intensity_sys_err"] = np.sqrt(bg_sys_err**2 + calib_sys_err**2)
+
+    # Step 4: Handle obs_date variable type conversion
     # TODO: Handle variable types correctly in RectangularSkyMap.build_cdf_dataset
     obs_date = map_ds["obs_date"]
     # Replace non-finite values with the int64 sentinel before casting
@@ -355,13 +369,13 @@ def calculate_all_rates_and_intensities(
     )
     map_ds["obs_date_range"] = xr.zeros_like(map_ds["obs_date"])
 
-    # Step 4: Swap esa_energy_step dimension for energy coordinate
+    # Step 5: Swap esa_energy_step dimension for energy coordinate
     map_ds = map_ds.swap_dims({"esa_energy_step": "energy"})
     map_ds = map_ds.drop_vars(
         ["esa_energy_step", "esa_energy_step_label"], errors="ignore"
     )
 
-    # Step 5: Apply Compton-Getting interpolation for heliocentric frame maps
+    # Step 6: Apply Compton-Getting interpolation for heliocentric frame maps
     if descriptor.frame_descriptor == "hf":
         logger.debug("Applying Compton-Getting interpolation for heliocentric frame")
         # Convert energy coordinate from keV to eV for interpolation
@@ -378,13 +392,6 @@ def calculate_all_rates_and_intensities(
         # Drop any esa_energy_step_label that may have been re-added
         map_ds = map_ds.drop_vars(["esa_energy_step_label"], errors="ignore")
 
-    # Step 6: Add calibration systematic uncertainty in quadrature with the
-    # background-associated systematic. This is a percentage of the intensity.
-    logger.debug("Adding calibration systematic uncertainty")
-    bg_sys_err = map_ds["ena_intensity_sys_err"]
-    calib_sys_err = CALIBRATION_UNCERTAINTY_FRACTION * map_ds["ena_intensity"]
-    map_ds["ena_intensity_sys_err"] = np.sqrt(bg_sys_err**2 + calib_sys_err**2)
-
     # Step 7: Clean up intermediate variables
     map_ds = cleanup_intermediate_variables(map_ds)
 
@@ -398,7 +405,7 @@ def calculate_ena_signal_rates(map_ds: xr.Dataset) -> xr.Dataset:
     Parameters
     ----------
     map_ds : xarray.Dataset
-        Map dataset that has counts, exposure_factor, and bg_rate calculated.
+        Map dataset that has ena_count, exposure_factor, and bg_rate calculated.
 
     Returns
     -------
@@ -409,7 +416,7 @@ def calculate_ena_signal_rates(map_ds: xr.Dataset) -> xr.Dataset:
     with np.errstate(divide="ignore"):
         # Calculate the ENA Signal Rate
         map_ds["ena_signal_rates"] = (
-            map_ds["counts"] / map_ds["exposure_factor"] - map_ds["bg_rate"]
+            map_ds["ena_count"] / map_ds["exposure_factor"] - map_ds["bg_rate"]
         )
         # Calculate the ENA Signal Rate Uncertainties
         # The minimum count uncertainty is 1 for any pixel that has non-zero
@@ -417,7 +424,7 @@ def calculate_ena_signal_rates(map_ds: xr.Dataset) -> xr.Dataset:
         # we can ignore the non-zero exposure time condition when setting the
         # minimum count uncertainty because division by zero exposure time results
         # in the correct NaN value.
-        min_counts_unc = xr.ufuncs.maximum(map_ds["counts"], 1)
+        min_counts_unc = xr.ufuncs.maximum(map_ds["ena_count"], 1)
         map_ds["ena_signal_rate_stat_unc"] = (
             np.sqrt(min_counts_unc) / map_ds["exposure_factor"]
         )
@@ -454,7 +461,7 @@ def calculate_ena_intensity(
     -------
     map_ds : xarray.Dataset
         Map dataset with new variables: ena_intensity, ena_intensity_stat_uncert,
-        ena_intensity_sys_err.
+        bg_intensity_sys_err.
     """
     # read calibration product configuration file
     cal_prod_df = CalibrationProductConfig.from_csv(l2_ancillary_path_dict["cal-prod"])
@@ -478,9 +485,7 @@ def calculate_ena_intensity(
 
     # Convert the exposure-time weighted average of background rate systematic
     # uncertainty from rate to intensity units.
-    map_ds["ena_intensity_sys_err"] = (
-        map_ds["bg_rate_sys_err"] / flux_conversion_divisor
-    )
+    map_ds["bg_intensity_sys_err"] = map_ds["bg_rate_sys_err"] / flux_conversion_divisor
 
     # Combine calibration products using proper weighted averaging
     # as described in Hi Algorithm Document Section 3.1.2
@@ -494,6 +499,7 @@ def calculate_ena_intensity(
         # Flux correction
         corrector = PowerLawFluxCorrector(l2_ancillary_path_dict["esa-eta-fit-factors"])
         # Apply flux correction with xarray inputs
+        pre_correction_intensity = map_ds["ena_intensity"]
         map_ds["ena_intensity"], map_ds["ena_intensity_stat_uncert"] = (
             corrector.apply_flux_correction(
                 map_ds["ena_intensity"],
@@ -501,6 +507,13 @@ def calculate_ena_intensity(
                 esa_energy,
             )
         )
+        # Scale the background systematic error by the same flux correction
+        # ratio that was just applied to ena_intensity.
+        with np.errstate(divide="ignore", invalid="ignore"):
+            flux_correction_ratio = map_ds["ena_intensity"] / pre_correction_intensity
+        map_ds["bg_intensity_sys_err"] = (
+            map_ds["bg_intensity_sys_err"] * flux_correction_ratio
+        )
 
     return map_ds
 
@@ -530,11 +543,12 @@ def combine_calibration_products(
     -------
     map_ds : xarray.Dataset
         Map dataset with updated variables: ena_intensity, ena_intensity_stat_uncert,
-        ena_intensity_sys_err now combined across calibration products at each
+        bg_intensity_sys_err, ena_count, bg_rate,
+        bg_rate_sys_err now combined across calibration products at each
         energy level.
     """
     ena_flux = map_ds["ena_intensity"]
-    sys_err = map_ds["ena_intensity_sys_err"]
+    sys_err = map_ds["bg_intensity_sys_err"]
 
     # Calculate improved statistical variance estimates using geometric factor
     # ratios to reduce bias from Poisson uncertainty estimation
@@ -560,10 +574,36 @@ def combine_calibration_products(
     )
     # For systematic error, just do quadrature sum over the systematic error for
     # each calibration product.
-    map_ds["ena_intensity_sys_err"] = np.sqrt(
+    map_ds["bg_intensity_sys_err"] = np.sqrt(
         (sys_err**2).sum(dim="calibration_prod", skipna=True, min_count=1)
     )
 
+    # ena_count is a diagnostic tally of binned direct events, so it is simply
+    # summed across calibration products.
+    map_ds["ena_count"] = map_ds["ena_count"].sum(dim="calibration_prod")
+
+    # bg_rate is combined the same way as ena_intensity_stat_uncert above:
+    # inverse-variance weighted using each calibration product's own
+    # bg_rate_sys_err as its uncertainty.
+    with np.errstate(divide="ignore", invalid="ignore"):
+        bg_rate_weights = xr.where(
+            map_ds["bg_rate_sys_err"] > 0,
+            1.0 / (map_ds["bg_rate_sys_err"] ** 2),
+            0.0,
+        )
+        weight_sum = bg_rate_weights.sum(
+            dim="calibration_prod", skipna=True, min_count=1
+        )
+        weighted_bg_sum = (map_ds["bg_rate"] * bg_rate_weights).sum(
+            dim="calibration_prod", skipna=True, min_count=1
+        )
+        map_ds["bg_rate"] = xr.where(
+            weight_sum > 0, weighted_bg_sum / weight_sum, np.nan
+        )
+        map_ds["bg_rate_sys_err"] = xr.where(
+            weight_sum > 0, np.sqrt(1 / weight_sum), np.nan
+        )
+
     return map_ds
 
 
@@ -607,8 +647,8 @@ def combine_maps(sky_maps: dict[str, RectangularSkyMap]) -> RectangularSkyMap:
     combined_map = sky_maps["ram"]
     combined = ram_ds.copy()
 
-    # Additive variables: counts and exposure_factor
-    combined["counts"] = ram_ds["counts"] + anti_ds["counts"]
+    # Additive variables: ena_count and exposure_factor
+    combined["ena_count"] = ram_ds["ena_count"] + anti_ds["ena_count"]
     combined["exposure_factor"] = ram_ds["exposure_factor"] + anti_ds["exposure_factor"]
 
     # Compute weights for ram and anti-ram based on the inverse of the variance
@@ -642,16 +682,24 @@ def combine_maps(sky_maps: dict[str, RectangularSkyMap]) -> RectangularSkyMap:
         # ena_intensity_stat_uncertainty is combined using inverse quadrature sum
         combined["ena_intensity_stat_uncert"] = np.sqrt(1 / total_weight)
 
-    # Exposure-weighted average for systematic error
-    # NaNs in the systematic error should occur only where the exposure_factor
+    # Exposure-weighted average for systematic error and background/calibration
+    # rate and systematic error variables.
+    # NaNs in these variables should occur only where the exposure_factor
     # is zero. This means the correct NaN handling is to just replace NaNs in
-    # the systematic error with zeros so that the sum is not affected.
+    # these variables with zeros so that the sum is not affected.
     with np.errstate(divide="ignore", invalid="ignore"):
         total_exp = combined["exposure_factor"]
-        combined["ena_intensity_sys_err"] = (
-            ram_ds["ena_intensity_sys_err"].fillna(0) * ram_ds["exposure_factor"]
-            + anti_ds["ena_intensity_sys_err"].fillna(0) * anti_ds["exposure_factor"]
-        ) / total_exp
+        for var in (
+            "ena_intensity_sys_err",
+            "bg_rate",
+            "bg_rate_sys_err",
+            "bg_intensity_sys_err",
+            "ena_intensity_calibration_sys_err",
+        ):
+            combined[var] = (
+                ram_ds[var].fillna(0) * ram_ds["exposure_factor"]
+                + anti_ds[var].fillna(0) * anti_ds["exposure_factor"]
+            ) / total_exp
 
     # Exposure-weighted average for obs_date
     with np.errstate(divide="ignore", invalid="ignore"):
@@ -729,7 +777,7 @@ def _calculate_improved_stat_variance(
 
     logger.debug("Computing geometric factor normalized signal rates")
 
-    # signal_rates = counts / exposure_factor - bg_rate
+    # signal_rates = ena_count / exposure_factor - bg_rate
     # signal_rates shape is: (n_epoch, n_energy, n_cal_prod, n_spatial_pixels)
     signal_rates = map_ds["ena_signal_rates"]
 
@@ -783,8 +831,6 @@ def cleanup_intermediate_variables(dataset: xr.Dataset) -> xr.Dataset:
     """
     # Remove the intermediate variables from the map
     potential_vars = [
-        "bg_rate",
-        "bg_rate_sys_err",
         "energy_sc",
         "ena_signal_rates",
         "ena_signal_rate_stat_unc",
diff --git a/imap_processing/tests/ena_maps/test_ena_maps.py b/imap_processing/tests/ena_maps/test_ena_maps.py
index 76113f466..7f9ad1917 100644
--- a/imap_processing/tests/ena_maps/test_ena_maps.py
+++ b/imap_processing/tests/ena_maps/test_ena_maps.py
@@ -192,8 +192,9 @@ def test_init(self, hi_pset_cdf_path):
         np.testing.assert_array_equal(hi_pset.az_el_points.shape, (3600, 2))
         # check that the midpoint_j2000_et property is equal to the expected value
         assert hi_pset.midpoint_j2000_et == ttj2000ns_to_et(hi_pset.epoch + delta / 2)
-        for var_name in ["exposure_factor", "bg_rate", "bg_rate_sys_err"]:
+        for var_name in ["exposure_factor", "bg_rate", "bg_rate_sys_err", "ena_count"]:
             assert var_name in hi_pset.data
+        assert "counts" not in hi_pset.data
 
     def test_from_cdf(self, hi_pset_cdf_path):
         """Test coverage for instantiating HiPointingSet from cdf."""
@@ -208,8 +209,8 @@ def test_plays_nice_with_rectangular_sky_map(self, hi_pset_cdf_path):
         rect_map = ena_maps.RectangularSkyMap(
             spacing_deg=2, spice_frame=geometry.SpiceFrame.IMAP_HAE
         )
-        rect_map.project_pset_values_to_map(hi_pset, ["counts", "exposure_factor"])
-        assert rect_map.data_1d["counts"].max() > 0
+        rect_map.project_pset_values_to_map(hi_pset, ["ena_count", "exposure_factor"])
+        assert rect_map.data_1d["ena_count"].max() > 0
 
 
 @pytest.fixture
diff --git a/imap_processing/tests/hi/test_hi_l2.py b/imap_processing/tests/hi/test_hi_l2.py
index 42a274511..755c6bc3d 100644
--- a/imap_processing/tests/hi/test_hi_l2.py
+++ b/imap_processing/tests/hi/test_hi_l2.py
@@ -128,7 +128,7 @@ def sample_map_dataset():
             "ena_intensity_stat_uncert": xr.DataArray(
                 np.random.rand(*shape) * 10 + 5, dims=list(coords.keys())
             ),
-            "ena_intensity_sys_err": xr.DataArray(
+            "bg_intensity_sys_err": xr.DataArray(
                 np.random.rand(*shape) * 5 + 1, dims=list(coords.keys())
             ),
             "bg_rate": xr.DataArray(
@@ -140,6 +140,9 @@ def sample_map_dataset():
             "exposure_factor": xr.DataArray(
                 np.random.rand(*shape) * 5 + 1, dims=list(coords.keys())
             ),
+            "ena_count": xr.DataArray(
+                np.random.randint(0, 100, size=shape), dims=list(coords.keys())
+            ),
         },
         coords=coords,
     )
@@ -177,10 +180,23 @@ def test_hi_l2(
     assert l2_dataset.attrs["Logical_source"] == f"imap_hi_l2_{descriptor_str}"
     assert "Hi90" in l2_dataset.attrs["Logical_source_description"]
 
-    assert len(l2_dataset.data_vars) == 16
+    assert len(l2_dataset.data_vars) == 21
     np.testing.assert_array_equal(
         l2_dataset["ena_intensity"].dims, ["epoch", "energy", "longitude", "latitude"]
     )
+    # ena_count, bg_rate, bg_rate_sys_err, and the two split-out systematic
+    # error components must survive to the final CDF output (previously
+    # "counts" was silently dropped for lacking a CDF attribute definition,
+    # and bg_rate/bg_rate_sys_err were dropped as "intermediate" variables).
+    for var_name in [
+        "ena_count",
+        "bg_rate",
+        "bg_rate_sys_err",
+        "bg_intensity_sys_err",
+        "ena_intensity_calibration_sys_err",
+    ]:
+        assert var_name in l2_dataset.data_vars
+    assert "counts" not in l2_dataset.data_vars
     # Test ISTP compliance by writing the CDF
     write_cdf(l2_dataset, istp=True)
 
@@ -268,7 +284,7 @@ def test_create_sky_map_from_psets(
         assert "energy" in sky_map.data_1d.coords
 
         # Test that we got some non-zero values
-        for var_name in ["counts", "exposure_factor", "obs_date"]:
+        for var_name in ["ena_count", "exposure_factor", "obs_date"]:
             assert var_name in sky_map.data_1d.data_vars
             assert np.nanmax(sky_map.data_1d[var_name].data) > 0
 
@@ -320,10 +336,10 @@ def test_calculate_ena_signal_rates(empty_rectangular_map_dataset):
     # we ensure that each unique combination is encountered in a PSET bin.
     map_ds.update(
         {
-            "counts": xr.DataArray(
+            "ena_count": xr.DataArray(
                 np.arange(np.prod(tuple(map_ds.sizes.values()))).reshape(counts_shape)
                 % 5,
-                name="counts",
+                name="ena_count",
                 dims=list(map_ds.sizes.keys()),
             ),
             "exposure_factor": xr.DataArray(
@@ -349,7 +365,7 @@ def test_calculate_ena_signal_rates(empty_rectangular_map_dataset):
         assert var_name in result_ds
         assert result_ds[var_name].shape == counts_shape
     # Verify that there are no negative signal rates. The synthetic data combination
-    # where counts = 0, exposure_factor = 1, and bg_rate = 1 would result in
+    # where ena_count = 0, exposure_factor = 1, and bg_rate = 1 would result in
     # an ena_signal_rate of (0 / 1) - 1 = -1
     assert np.nanmin(result_ds["ena_signal_rates"].values) >= 0
     # Verify that the minimum finite uncertainty is sqrt(1) / exposure_factor.
@@ -358,7 +374,7 @@ def test_calculate_ena_signal_rates(empty_rectangular_map_dataset):
     assert np.nanmin(result_ds["ena_signal_rate_stat_unc"].values) == 1 / 2
 
 
-@pytest.fixture(scope="module")
+@pytest.fixture
 def ena_intensity_map_ds(empty_rectangular_map_dataset):
     """Fixture that produces a dataset to use in testing ena_intensity."""
     # Start with an empty (coords only) dataset
@@ -379,10 +395,16 @@ def ena_intensity_map_ds(empty_rectangular_map_dataset):
                 dims=list(map_ds.sizes.keys()),
             ),
             "bg_rate_sys_err": xr.DataArray(
-                np.arange(np.prod(tuple(map_ds.sizes.values()))).reshape(var_shape) % 3,
+                np.arange(np.prod(tuple(map_ds.sizes.values()))).reshape(var_shape) % 3
+                + 1,
                 name="bg_rate_sys_err",
                 dims=list(map_ds.sizes.keys()),
             ),
+            "ena_count": xr.DataArray(
+                np.arange(np.prod(tuple(map_ds.sizes.values()))).reshape(var_shape),
+                name="ena_count",
+                dims=list(map_ds.sizes.keys()),
+            ),
         }
     )
 
@@ -417,7 +439,7 @@ def test_calculate_ena_intensity(ena_intensity_map_ds, anc_path_dict):
     for var_name in [
         "ena_intensity",
         "ena_intensity_stat_uncert",
-        "ena_intensity_sys_err",
+        "bg_intensity_sys_err",
     ]:
         assert var_name in result_ds
         # Check that calibration_prod dimension has been removed
@@ -457,6 +479,47 @@ def test_calculate_ena_intensity_flux_correction_logic(
         mock_instance.apply_flux_correction.assert_not_called()
 
 
+@mock.patch("imap_processing.hi.hi_l2.PowerLawFluxCorrector", autospec=True)
+def test_calculate_ena_intensity_scales_background_systematic_by_flux_ratio(
+    mock_flux_corrector_class, ena_intensity_map_ds, anc_path_dict
+):
+    """Test that bg_intensity_sys_err is flux-corrected.
+
+    The background-derived systematic error is computed before the flux
+    correction is applied to ena_intensity, so it must be rescaled by the same
+    ratio the correction applies to ena_intensity, or it will become
+    inconsistent with the corrected intensity.
+    """
+    # Correction doubles the intensity (and leaves stat_unc unchanged) so the
+    # resulting flux correction ratio is exactly 2.0 everywhere.
+    mock_instance = mock_flux_corrector_class.return_value
+    mock_instance.apply_flux_correction.side_effect = (
+        lambda intensity, stat_unc, energy: (intensity * 2.0, stat_unc)
+    )
+
+    # "raw" descriptor skips flux correction entirely, giving the uncorrected
+    # baseline value of bg_intensity_sys_err to compare against.
+    raw_descriptor = MapDescriptor.from_string("h90-enaraw-h-sf-nsp-full-gcs-6deg-3mo")
+    uncorrected_ds = ena_intensity_map_ds.copy(deep=True)
+    uncorrected_result = calculate_ena_intensity(
+        uncorrected_ds, anc_path_dict, raw_descriptor
+    )
+    expected_uncorrected_bg_sys_err = uncorrected_result[
+        "bg_intensity_sys_err"
+    ].values.copy()
+
+    map_descriptor = MapDescriptor.from_string("h90-ena-h-sf-nsp-full-gcs-6deg-3mo")
+    corrected_ds = ena_intensity_map_ds.copy(deep=True)
+    corrected_result = calculate_ena_intensity(
+        corrected_ds, anc_path_dict, map_descriptor
+    )
+
+    np.testing.assert_allclose(
+        corrected_result["bg_intensity_sys_err"].values,
+        expected_uncorrected_bg_sys_err * 2.0,
+    )
+
+
 def test_combine_calibration_products(sample_map_dataset):
     """Test coverage for combine_calibration_products"""
     test_ds, geometric_factors, esa_energies = sample_map_dataset
@@ -472,7 +535,7 @@ def test_combine_calibration_products(sample_map_dataset):
     expected_vars = [
         "ena_intensity",
         "ena_intensity_stat_uncert",
-        "ena_intensity_sys_err",
+        "bg_intensity_sys_err",
     ]
     for var_name in expected_vars:
         assert var_name in result_ds
@@ -497,9 +560,9 @@ def test_combine_calibration_products(sample_map_dataset):
     )
 
     # Check systematic error combination (root sum of squares)
-    input_sys_err = test_ds["ena_intensity_sys_err"]
+    input_sys_err = test_ds["bg_intensity_sys_err"]
     expected_sys_err = np.sqrt((input_sys_err**2).sum(dim="calibration_prod"))
-    combined_sys_err = result_ds["ena_intensity_sys_err"]
+    combined_sys_err = result_ds["bg_intensity_sys_err"]
 
     np.testing.assert_array_almost_equal(
         combined_sys_err.values, expected_sys_err.values, decimal=10
@@ -592,7 +655,7 @@ def test_weighted_average_mathematical_correctness():
             "ena_intensity_stat_uncert": xr.DataArray(
                 stat_unc_values, dims=list(coords.keys())
             ),
-            "ena_intensity_sys_err": xr.DataArray(
+            "bg_intensity_sys_err": xr.DataArray(
                 sys_err_values, dims=list(coords.keys())
             ),
             "ena_signal_rates": xr.DataArray(
@@ -607,6 +670,12 @@ def test_weighted_average_mathematical_correctness():
                 np.array([2.0]).reshape(1, 1, 1, 1),
                 dims=[d for d in coords.keys() if d != "calibration_prod"],
             ),
+            "bg_rate_sys_err": xr.DataArray(
+                np.array([1.0, 2.0]).reshape(1, 1, 2, 1, 1), dims=list(coords.keys())
+            ),
+            "ena_count": xr.DataArray(
+                np.array([50.0, 100.0]).reshape(1, 1, 2, 1, 1), dims=list(coords.keys())
+            ),
         }
     )
 
@@ -621,12 +690,12 @@ def test_weighted_average_mathematical_correctness():
     # Check that results are finite and reasonable
     assert np.isfinite(result_ds["ena_intensity"].values[0, 0, 0, 0])
     assert result_ds["ena_intensity_stat_uncert"].values[0, 0, 0, 0] > 0
-    assert result_ds["ena_intensity_sys_err"].values[0, 0, 0, 0] > 0
+    assert result_ds["bg_intensity_sys_err"].values[0, 0, 0, 0] > 0
 
     # Systematic error should be root sum of squares
     expected_sys_err = np.sqrt(5.0**2 + 10.0**2)
     np.testing.assert_almost_equal(
-        result_ds["ena_intensity_sys_err"].values[0, 0, 0, 0],
+        result_ds["bg_intensity_sys_err"].values[0, 0, 0, 0],
         expected_sys_err,
         decimal=10,
     )
@@ -653,7 +722,7 @@ def test_statistical_uncertainty_combination_correctness():
             "ena_intensity_stat_uncert": xr.DataArray(
                 stat_unc_values, dims=list(coords.keys())
             ),
-            "ena_intensity_sys_err": xr.DataArray(
+            "bg_intensity_sys_err": xr.DataArray(
                 sys_err_values, dims=list(coords.keys())
             ),
             "ena_signal_rates": xr.DataArray(flux_values, dims=list(coords.keys())),
@@ -663,6 +732,12 @@ def test_statistical_uncertainty_combination_correctness():
             "exposure_factor": xr.DataArray(
                 np.array([1.0, 1.0]).reshape(1, 1, 2, 1, 1), dims=list(coords.keys())
             ),
+            "bg_rate_sys_err": xr.DataArray(
+                np.array([0.5, 1.0]).reshape(1, 1, 2, 1, 1), dims=list(coords.keys())
+            ),
+            "ena_count": xr.DataArray(
+                np.array([90.0, 210.0]).reshape(1, 1, 2, 1, 1), dims=list(coords.keys())
+            ),
         }
     )
 
@@ -711,9 +786,18 @@ def test_combine_calibration_products_edge_cases():
             "ena_intensity_stat_uncert": xr.DataArray(
                 np.array([10.0]).reshape(1, 1, 1, 1, 1), dims=list(coords.keys())
             ),
-            "ena_intensity_sys_err": xr.DataArray(
+            "bg_intensity_sys_err": xr.DataArray(
                 np.array([5.0]).reshape(1, 1, 1, 1, 1), dims=list(coords.keys())
             ),
+            "bg_rate": xr.DataArray(
+                np.array([1.0]).reshape(1, 1, 1, 1, 1), dims=list(coords.keys())
+            ),
+            "bg_rate_sys_err": xr.DataArray(
+                np.array([0.5]).reshape(1, 1, 1, 1, 1), dims=list(coords.keys())
+            ),
+            "ena_count": xr.DataArray(
+                np.array([100.0]).reshape(1, 1, 1, 1, 1), dims=list(coords.keys())
+            ),
         }
     )
 
@@ -730,7 +814,14 @@ def test_combine_calibration_products_edge_cases():
     np.testing.assert_almost_equal(result_ds["ena_intensity"].values[0, 0, 0, 0], 100.0)
 
     # Check that calibration_prod dimension was removed
-    for var in ["ena_intensity", "ena_intensity_stat_uncert", "ena_intensity_sys_err"]:
+    for var in [
+        "ena_intensity",
+        "ena_intensity_stat_uncert",
+        "bg_intensity_sys_err",
+        "bg_rate",
+        "bg_rate_sys_err",
+        "ena_count",
+    ]:
         assert "calibration_prod" not in result_ds[var].dims
 
 
@@ -761,6 +852,7 @@ def test_combine_calibration_products_nan_handling():
     bg_rate = np.full(shape, 20.0)
     bg_rate_sys_err = np.full(shape, 2.0)
     exposure_factor = np.full(shape, 1.0)
+    ena_count = np.full(shape, 250.0)
 
     # Set NaN in one calibration product's uncertainty at position [0,0,0,0,0]
     # The other calibration product is valid, so result should be finite
@@ -783,11 +875,12 @@ def test_combine_calibration_products_nan_handling():
         {
             "ena_intensity": xr.DataArray(intensity, dims=dim_names),
             "ena_intensity_stat_uncert": xr.DataArray(stat_uncert, dims=dim_names),
-            "ena_intensity_sys_err": xr.DataArray(sys_err, dims=dim_names),
+            "bg_intensity_sys_err": xr.DataArray(sys_err, dims=dim_names),
             "ena_signal_rates": xr.DataArray(signal_rates, dims=dim_names),
             "bg_rate": xr.DataArray(bg_rate, dims=dim_names),
             "bg_rate_sys_err": xr.DataArray(bg_rate_sys_err, dims=dim_names),
             "exposure_factor": xr.DataArray(exposure_factor, dims=dim_names),
+            "ena_count": xr.DataArray(ena_count, dims=dim_names),
         }
     )
 
@@ -808,10 +901,10 @@ def test_combine_calibration_products_nan_handling():
     assert np.isnan(result_ds["ena_intensity_stat_uncert"].values[0, 1, 0, 0])
 
     # Test 3: When one product has NaN sys_err, result should be finite
-    assert np.isfinite(result_ds["ena_intensity_sys_err"].values[0, 0, 1, 0])
+    assert np.isfinite(result_ds["bg_intensity_sys_err"].values[0, 0, 1, 0])
 
     # Test 4: When ALL products have NaN sys_err, result should be NaN
-    assert np.isnan(result_ds["ena_intensity_sys_err"].values[0, 1, 1, 0])
+    assert np.isnan(result_ds["bg_intensity_sys_err"].values[0, 1, 1, 0])
 
 
 # =============================================================================
@@ -898,10 +991,12 @@ def test_process_single_pset_renames_variables(
     assert "exposure_factor" in result
     assert "bg_rate" in result
     assert "bg_rate_sys_err" in result
+    assert "ena_count" in result
     # Original names should not exist
     assert "exposure_times" not in result
     assert "background_rates" not in result
     assert "background_rates_uncertainty" not in result
+    assert "counts" not in result
 
 
 @mock.patch("imap_processing.hi.hi_l2.calculate_ram_mask")
@@ -1079,7 +1174,7 @@ def mock_map_dataset_for_rates():
 
     map_ds = xr.Dataset(
         {
-            "counts": xr.DataArray(
+            "ena_count": xr.DataArray(
                 np.ones(shape) * 100.0, dims=list(coords.keys())[:5]
             ),
             "exposure_factor": xr.DataArray(
@@ -1249,15 +1344,16 @@ def test_cleanup_intermediate_variables():
     result = cleanup_intermediate_variables(ds)
 
     # Intermediate variables should be removed
-    assert "bg_rate" not in result
-    assert "bg_rate_sys_err" not in result
     assert "energy_sc" not in result
     assert "ena_signal_rates" not in result
     assert "ena_signal_rate_stat_unc" not in result
 
-    # Non-intermediate variables should remain
+    # Non-intermediate variables should remain. bg_rate/bg_rate_sys_err are now
+    # final output variables, not intermediates.
     assert "ena_intensity" in result
     assert "exposure_factor" in result
+    assert "bg_rate" in result
+    assert "bg_rate_sys_err" in result
 
 
 def test_cleanup_intermediate_variables_missing_vars():
@@ -1265,7 +1361,7 @@ def test_cleanup_intermediate_variables_missing_vars():
     # Create a dataset without all intermediate variables
     ds = xr.Dataset(
         {
-            "bg_rate": xr.DataArray([1, 2, 3], dims=["x"]),
+            "energy_sc": xr.DataArray([1, 2, 3], dims=["x"]),
             "ena_intensity": xr.DataArray([10, 20, 30], dims=["x"]),
         }
     )
@@ -1273,7 +1369,7 @@ def test_cleanup_intermediate_variables_missing_vars():
     # Should not raise an error
     result = cleanup_intermediate_variables(ds)
 
-    assert "bg_rate" not in result
+    assert "energy_sc" not in result
     assert "ena_intensity" in result
 
 
@@ -1295,7 +1391,7 @@ def _create_map(intensity_offset=0, exposure_offset=0):
         shape = (1, 3, 4, 2)  # epoch, energy, lon, lat
         sky_map.data_1d = xr.Dataset(
             {
-                "counts": xr.DataArray(
+                "ena_count": xr.DataArray(
                     np.ones(shape) * (100 + intensity_offset),
                     dims=["epoch", "energy", "longitude", "latitude"],
                 ),
@@ -1323,6 +1419,22 @@ def _create_map(intensity_offset=0, exposure_offset=0):
                     np.ones(shape) * 2.0,
                     dims=["epoch", "energy", "longitude", "latitude"],
                 ),
+                "bg_rate": xr.DataArray(
+                    np.ones(shape) * (3.0 + 0.1 * intensity_offset),
+                    dims=["epoch", "energy", "longitude", "latitude"],
+                ),
+                "bg_rate_sys_err": xr.DataArray(
+                    np.ones(shape) * (0.5 + 0.01 * intensity_offset),
+                    dims=["epoch", "energy", "longitude", "latitude"],
+                ),
+                "bg_intensity_sys_err": xr.DataArray(
+                    np.ones(shape) * (1.5 + 0.05 * intensity_offset),
+                    dims=["epoch", "energy", "longitude", "latitude"],
+                ),
+                "ena_intensity_calibration_sys_err": xr.DataArray(
+                    np.ones(shape) * (1.0 + 0.02 * intensity_offset),
+                    dims=["epoch", "energy", "longitude", "latitude"],
+                ),
             },
             coords={
                 "epoch": [0],
@@ -1361,8 +1473,8 @@ def test_combine_maps_two_maps(mock_sky_map_for_combine):
     # Check additive variables
     expected_counts = 100 + 120  # 100 + (100 + 20)
     np.testing.assert_array_almost_equal(
-        result.data_1d["counts"].values,
-        np.ones_like(result.data_1d["counts"].values) * expected_counts,
+        result.data_1d["ena_count"].values,
+        np.ones_like(result.data_1d["ena_count"].values) * expected_counts,
     )
 
     expected_exposure = 10 + 15  # 10 + (10 + 5)
@@ -1690,13 +1802,13 @@ def test_calculate_ena_intensity_uses_bg_rate_sys_err(
         ena_intensity_map_ds, anc_path_dict, map_descriptor
     )
 
-    # Verify ena_intensity_sys_err was calculated
-    assert "ena_intensity_sys_err" in result_ds
+    # Verify bg_intensity_sys_err was calculated
+    assert "bg_intensity_sys_err" in result_ds
 
     # The sys_err should be based on bg_rate_sys_err / (geometric_factor * energy)
     # After combine_calibration_products, it's combined in quadrature across cal prods
     # We just verify it's finite and positive where expected
-    sys_err = result_ds["ena_intensity_sys_err"]
+    sys_err = result_ds["bg_intensity_sys_err"]
     assert np.all(sys_err.values[np.isfinite(sys_err.values)] >= 0)
 
 
@@ -1721,11 +1833,34 @@ def test_calculate_all_rates_adds_calibration_systematic(
 
     # Verify ena_intensity_sys_err includes calibration systematic
     assert "ena_intensity_sys_err" in result_ds
+    assert "bg_intensity_sys_err" in result_ds
+    assert "ena_intensity_calibration_sys_err" in result_ds
 
     # The sys_err should be larger than CALIBRATION_UNCERTAINTY_FRACTION * intensity
     # because it includes both bg systematic and calibration systematic in quadrature
     intensity = result_ds["ena_intensity"]
     sys_err = result_ds["ena_intensity_sys_err"]
+    bg_sys_err = result_ds["bg_intensity_sys_err"]
+    calib_sys_err = result_ds["ena_intensity_calibration_sys_err"]
+
+    # calib_sys_err should be exactly CALIBRATION_UNCERTAINTY_FRACTION of intensity
+    np.testing.assert_allclose(
+        calib_sys_err.values,
+        CALIBRATION_UNCERTAINTY_FRACTION * intensity.values,
+    )
+
+    # The combined sys_err should be the quadrature sum of the two components
+    valid_mask = (
+        np.isfinite(sys_err.values)
+        & np.isfinite(bg_sys_err.values)
+        & np.isfinite(calib_sys_err.values)
+    )
+    np.testing.assert_allclose(
+        sys_err.values[valid_mask],
+        np.sqrt(
+            bg_sys_err.values[valid_mask] ** 2 + calib_sys_err.values[valid_mask] ** 2
+        ),
+    )
 
     # Minimum expected sys_err is 22% of intensity (if bg systematic were zero)
     min_expected = CALIBRATION_UNCERTAINTY_FRACTION * np.abs(intensity)