diff --git a/pipelines/visitQualityExtended.yaml b/pipelines/visitQualityExtended.yaml
index 9a396b09d..0fdaca864 100644
--- a/pipelines/visitQualityExtended.yaml
+++ b/pipelines/visitQualityExtended.yaml
@@ -5,3 +5,11 @@ tasks:
     class: lsst.analysis.tools.tasks.astrometricCatalogMatch.AstrometricCatalogMatchVisitTask
   refCatSourceVisit:
     class: lsst.analysis.tools.tasks.refCatSourceAnalysis.RefCatSourceAnalysisTask
+  deltaSkyCorrHist:
+    class: lsst.analysis.tools.tasks.deltaSkyCorrAnalysis.DeltaSkyCorrHistTask
+  deltaSkyCorrAnalysis:
+    class: lsst.analysis.tools.tasks.deltaSkyCorrAnalysis.DeltaSkyCorrAnalysisTask
+    config:
+      atools.deltaSkyCorr: DeltaSkyCorrXYPlot
+      python: |
+        from lsst.analysis.tools.atools import *
diff --git a/python/lsst/analysis/tools/actions/scalar/scalarActions.py b/python/lsst/analysis/tools/actions/scalar/scalarActions.py
index 4756e45d4..56ad716fb 100644
--- a/python/lsst/analysis/tools/actions/scalar/scalarActions.py
+++ b/python/lsst/analysis/tools/actions/scalar/scalarActions.py
@@ -36,6 +36,8 @@
     "FracInRange",
     "FracNan",
     "SumAction",
+    "MedianHistAction",
+    "IqrHistAction",
 )
 
 import operator
@@ -246,3 +248,90 @@ def __call__(self, data: KeyedData, **kwargs) -> Scalar:
         mask = self.getMask(**kwargs)
         arr = cast(Vector, data[self.vectorKey.format(**kwargs)])[mask]
         return cast(Scalar, np.nansum(arr))
+
+
+class MedianHistAction(ScalarAction):
+    """Calculates the median of the given histogram data."""
+
+    histKey = Field[str]("Key of frequency Vector")
+    midKey = Field[str]("Key of bin midpoints Vector")
+
+    def getInputSchema(self) -> KeyedDataSchema:
+        return (
+            (self.histKey, Vector),
+            (self.midKey, Vector),
+        )
+
+    def histMedian(self, hist, bin_mid):
+        """Calculates the median of a histogram with binned values
+
+        Parameters
+        ----------
+        hist : `numpy.ndarray`
+            Frequency array
+        bin_mid : `numpy.ndarray`
+            Bin midpoints array
+
+        Returns
+        -------
+        median : `float`
+            Median of histogram with binned values
+        """
+        cumulative_sum = np.cumsum(hist)
+        median_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 2)
+        median = bin_mid[median_index]
+        return median
+
+    def __call__(self, data: KeyedData, **kwargs):
+        if len(data[self.histKey.format(**kwargs)]) != 0:
+            hist = cast(Vector, data[self.histKey.format(**kwargs)])
+            bin_mid = cast(Vector, data[self.midKey.format(**kwargs)])
+            med = cast(Scalar, float(self.histMedian(hist, bin_mid)))
+        else:
+            med = np.NaN
+        return med
+
+
+class IqrHistAction(ScalarAction):
+    """Calculates the interquartile range of the given histogram data."""
+
+    histKey = Field[str]("Key of frequency Vector")
+    midKey = Field[str]("Key of bin midpoints Vector")
+
+    def getInputSchema(self) -> KeyedDataSchema:
+        return (
+            (self.histKey, Vector),
+            (self.midKey, Vector),
+        )
+
+    def histIqr(self, hist, bin_mid):
+        """Calculates the interquartile range of a histogram with binned values
+
+        Parameters
+        ----------
+        hist : `numpy.ndarray`
+            Frequency array
+        bin_mid : `numpy.ndarray`
+            Bin midpoints array
+
+        Returns
+        -------
+        iqr : `float`
+            Inter-quartile range of histogram with binned values
+        """
+        cumulative_sum = np.cumsum(hist)
+        liqr_index = np.searchsorted(cumulative_sum, cumulative_sum[-1] / 4)
+        uiqr_index = np.searchsorted(cumulative_sum, (3 / 4) * cumulative_sum[-1])
+        liqr = bin_mid[liqr_index]
+        uiqr = bin_mid[uiqr_index]
+        iqr = uiqr - liqr
+        return iqr
+
+    def __call__(self, data: KeyedData, **kwargs):
+        if len(data[self.histKey.format(**kwargs)]) != 0:
+            hist = cast(Vector, data[self.histKey.format(**kwargs)])
+            bin_mid = cast(Vector, data[self.midKey.format(**kwargs)])
+            iqr = cast(Scalar, float(self.histIqr(hist, bin_mid)))
+        else:
+            iqr = np.NaN
+        return iqr
diff --git a/python/lsst/analysis/tools/atools/__init__.py b/python/lsst/analysis/tools/atools/__init__.py
index 10c7668b1..048b3a7d6 100644
--- a/python/lsst/analysis/tools/atools/__init__.py
+++ b/python/lsst/analysis/tools/atools/__init__.py
@@ -1,6 +1,7 @@
 from .astrometricRepeatability import *
 from .coveragePlots import *
 from .deblenderMetric import *
+from .deltaSkyCorr import *
 from .diaSolarSystemObjectMetrics import *
 from .diaSourceMetrics import *
 from .diffimMetadataMetrics import *
diff --git a/python/lsst/analysis/tools/atools/deltaSkyCorr.py b/python/lsst/analysis/tools/atools/deltaSkyCorr.py
new file mode 100644
index 000000000..27f21e957
--- /dev/null
+++ b/python/lsst/analysis/tools/atools/deltaSkyCorr.py
@@ -0,0 +1,67 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+__all__ = ("DeltaSkyCorrXYPlot",)
+
+from ..actions.plot.xyPlot import XYPlot
+from ..actions.scalar.scalarActions import IqrHistAction, MedianHistAction
+from ..actions.vector import ConstantValue, LoadVector
+from ..interfaces import AnalysisTool
+
+
+class DeltaSkyCorrXYPlot(AnalysisTool):
+    parameterizedBand: bool = False
+
+    def setDefaults(self):
+        super().setDefaults()
+        self.process.buildActions.x = LoadVector()
+        self.process.buildActions.y = LoadVector()
+        self.process.buildActions.x.vectorKey = "bin_mid"
+        self.process.buildActions.y.vectorKey = "hist"
+        self.process.buildActions.xerr = ConstantValue(value=0)
+        self.process.buildActions.yerr = ConstantValue(value=0)
+
+        self.process.calculateActions.median = MedianHistAction()
+        self.process.calculateActions.median.histKey = "hist"
+        self.process.calculateActions.median.midKey = "bin_mid"
+        self.process.calculateActions.iqr = IqrHistAction()
+        self.process.calculateActions.iqr.histKey = "hist"
+        self.process.calculateActions.iqr.midKey = "bin_mid"
+
+        self.produce.plot = XYPlot()
+        self.produce.plot.xAxisLabel = r"$\Delta$skyCorr Flux (nJy)"
+        self.produce.plot.yAxisLabel = "Frequency"
+        self.produce.plot.yScale = "log"
+        self.produce.plot.xLine = 0
+        self.produce.plot.strKwargs = {
+            "fmt": "-",
+            "color": "b",
+        }
+
+        self.produce.metric.units = {
+            "median": "nJy",
+            "iqr": "nJy",
+        }
+
+        self.produce.metric.newNames = {
+            "median": "delta_skyCorr_median",
+            "iqr": "delta_skyCorr_iqr",
+        }
diff --git a/python/lsst/analysis/tools/tasks/deltaSkyCorrAnalysis.py b/python/lsst/analysis/tools/tasks/deltaSkyCorrAnalysis.py
new file mode 100644
index 000000000..31f5541ae
--- /dev/null
+++ b/python/lsst/analysis/tools/tasks/deltaSkyCorrAnalysis.py
@@ -0,0 +1,222 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+import logging
+
+import numpy as np
+from lsst.pex.config import Field, ListField
+from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections
+from lsst.pipe.base.connectionTypes import Input, Output
+
+from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
+
+log = logging.getLogger(__name__)
+
+
+class DeltaSkyCorrHistConnections(PipelineTaskConnections, dimensions=("instrument", "visit")):
+    """Connections class for DeltaSkyCorrHistTask."""
+
+    skyCorrs = Input(
+        name="skyCorr",
+        storageClass="Background",
+        doc="Sky correction background models from a run without any synthetic source injection.",
+        multiple=True,
+        dimensions=("instrument", "visit", "detector"),
+        deferLoad=True,
+    )
+    injected_skyCorrs = Input(
+        name="injected_skyCorr",
+        storageClass="Background",
+        doc="Sky correction background models from a run with synthetic sources injected into the data.",
+        multiple=True,
+        dimensions=("instrument", "visit", "detector"),
+        deferLoad=True,
+    )
+    calexpBackgrounds = Input(
+        name="calexpBackground",
+        storageClass="Background",
+        doc="Initial per-detector background models associated with the calibrated exposure.",
+        multiple=True,
+        dimensions=("instrument", "visit", "detector"),
+        deferLoad=True,
+    )
+    photoCalib = Input(
+        name="calexp.photoCalib",
+        storageClass="PhotoCalib",
+        doc="Photometric calibration associated with the calibrated exposure.",
+        multiple=True,
+        dimensions=("instrument", "visit", "detector"),
+        deferLoad=True,
+    )
+    delta_skyCorr_hist = Output(
+        name="delta_skyCorr_hist",
+        storageClass="ArrowNumpyDict",
+        doc="A dictionary containing the histogram values, bin mid points, and bin lower/upper edges for the "
+        "aggregated skyCorr difference dataset, i.e., the difference between the injected and non-injected "
+        "sky correction background models.",
+        dimensions=("instrument", "visit"),
+    )
+
+
+class DeltaSkyCorrHistConfig(PipelineTaskConfig, pipelineConnections=DeltaSkyCorrHistConnections):
+    """Config class for DeltaSkyCorrHistTask."""
+
+    bin_range = ListField[float](
+        doc="The lower and upper range for the histogram bins, in nJy.",
+        default=[-1, 1],
+    )
+    bin_width = Field[float](
+        doc="The width of each histogram bin, in nJy.",
+        default=0.0001,
+    )
+
+
+class DeltaSkyCorrHistTask(PipelineTask):
+    """A task for generating a histogram of counts in the difference image
+    between an injected sky correction frame and a non-injected sky correction
+    frame (i.e., injected_skyCorr - skyCorr).
+    """
+
+    ConfigClass = DeltaSkyCorrHistConfig
+    _DefaultName = "deltaSkyCorrHist"
+
+    def __init__(self, initInputs=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def runQuantum(self, butlerQC, inputRefs, outputRefs):
+        inputs = butlerQC.get(inputRefs)
+        inputs["num_initial_bgs"] = len(inputs["calexpBackgrounds"][0].get())
+        delta_skyCorr_hist = self.run(**{k: v for k, v in inputs.items() if k != "calexpBackgrounds"})
+        butlerQC.put(delta_skyCorr_hist, outputRefs.delta_skyCorr_hist)
+
+    def run(self, skyCorrs, injected_skyCorrs, num_initial_bgs, photoCalib):
+        """Generate a histogram of counts in the difference image between an
+        injected sky correction frame and a non-injected sky correction frame
+        (i.e., injected_skyCorr - skyCorr).
+
+        Parameters
+        ----------
+        skyCorrs : `list`[`~lsst.daf.butler.DeferredDatasetHandle`]
+            Sky correction background models from a run without any synthetic
+            source injection.
+            These deferred dataset handles should normally resolve to
+            `~lsst.afw.math.BackgroundList` objects.
+        injected_skyCorrs : `list`[`~lsst.daf.butler.DeferredDatasetHandle`]
+            Sky correction background models from a run with synthetic sources
+            injected into the data.
+            These deferred dataset handles should normally resolve to
+            `~lsst.afw.math.BackgroundList` objects.
+        num_initial_bgs : `int`
+            The length of the initial per-detector background model list.
+            This number of background models will be skipped from the start of
+            each skyCorr/injected_skyCorr background model list.
+            See the Notes section for more details.
+        photoCalib : `list`[`~lsst.daf.butler.DeferredDatasetHandle`]
+            Photometric calibration, for conversion from counts to nJy.
+
+        Returns
+        -------
+        delta_skyCorr_hist : `dict`[`str`, `~numpy.ndarray`]
+            A dictionary containing the histogram values and bin lower/upper
+            edges for the skyCorr difference dataset.
+
+        Notes
+        -----
+        The first N background elements in the skyCorr/injected_skyCorr
+        background list are the inverse of the initial per-detector background
+        solution.
+        The effect of this is that adding a sky correction frame to a
+        background-subtracted calibrated exposure will undo the per-detector
+        background solution and apply the full focal plane sky correction in
+        its place.
+
+        For this task, we only want to compare the extra (subtractive) sky
+        correction components, so we skip the first N background models from
+        the sky frame.
+        """
+        # Generate lookup tables for the skyCorr/injected_skyCorr data.
+        lookup_skyCorrs = {x.dataId: x for x in skyCorrs}
+        lookup_injected_skyCorrs = {x.dataId: x for x in injected_skyCorrs}
+        lookup_photoCalib = {x.dataId: x for x in photoCalib}
+
+        # Set up the global histogram.
+        bin_edges = np.arange(
+            self.config.bin_range[0],
+            self.config.bin_range[1] + self.config.bin_width,
+            self.config.bin_width,
+        )
+        hist = np.zeros(len(bin_edges) - 1)
+        log.info("Generating a histogram containing %d bins.", len(hist))
+
+        # Loop over the skyCorr/injected_skyCorr data.
+        for dataId in lookup_injected_skyCorrs.keys():
+            # Get the skyCorr/injected_skyCorr data.
+            skyCorr = lookup_skyCorrs[dataId].get()
+            injected_skyCorr = lookup_injected_skyCorrs[dataId].get()
+            # And the photometric calibration
+            instFluxToNanojansky = lookup_photoCalib[dataId].get().instFluxToNanojansky(1)
+
+            # Isolate the extra (subtractive) sky correction components.
+            skyCorr_extras = skyCorr.clone()
+            skyCorr_extras._backgrounds = skyCorr_extras._backgrounds[num_initial_bgs:]
+            injected_skyCorr_extras = injected_skyCorr.clone()
+            injected_skyCorr_extras._backgrounds = injected_skyCorr_extras._backgrounds[num_initial_bgs:]
+
+            # Create the delta_skyCorr array.
+            delta_skyCorr_det = injected_skyCorr_extras.getImage().array - skyCorr_extras.getImage().array
+            delta_skyCorr_det *= instFluxToNanojansky  # Convert image to nJy
+
+            # Compute the per-detector histogram; update the global histogram.
+            hist_det, _ = np.histogram(delta_skyCorr_det, bins=bin_edges)
+            hist += hist_det
+
+        # Return results.
+        num_populated_bins = len([x for x in hist if x == 0])
+        log.info("Populated %d of %d histogram bins.", len(hist) - num_populated_bins, len(hist))
+        bin_mid = bin_edges[:-1] + (self.config.bin_width / 2)
+        delta_skyCorr_hist = dict(
+            hist=hist, bin_lower=bin_edges[:-1], bin_upper=bin_edges[1:], bin_mid=bin_mid
+        )
+        return delta_skyCorr_hist
+
+
+class DeltaSkyCorrAnalysisConnections(
+    AnalysisBaseConnections,
+    dimensions=("instrument", "visit"),
+    defaultTemplates={"outputName": "deltaSkyCorr"},
+):
+    data = Input(
+        name="delta_skyCorr_hist",
+        storageClass="ArrowNumpyDict",
+        doc="A dictionary containing the histogram values, bin mid points, and bin lower/upper edges for the "
+        "aggregated skyCorr difference dataset, i.e., the difference between the injected and non-injected "
+        "sky correction background models.",
+        deferLoad=True,
+        dimensions=("instrument", "visit"),
+    )
+
+
+class DeltaSkyCorrAnalysisConfig(AnalysisBaseConfig, pipelineConnections=DeltaSkyCorrAnalysisConnections):
+    pass
+
+
+class DeltaSkyCorrAnalysisTask(AnalysisPipelineTask):
+    ConfigClass = DeltaSkyCorrAnalysisConfig
+    _DefaultName = "deltaSkyCorrAnalysis"