Add MedianHistAction; calc for deltaSkyCorr

pipelines/visitQualityExtended.yaml

@@ -5,3 +5,9 @@ tasks:
     class: lsst.analysis.tools.tasks.astrometricCatalogMatch.AstrometricCatalogMatchVisitTask
   refCatSourceVisit:
     class: lsst.analysis.tools.tasks.refCatSourceAnalysis.RefCatSourceAnalysisTask
+  deltaSkyCorrVisit:
+    class: lsst.analysis.tools.tasks.deltaSkyCorrAnalysis.DeltaSkyCorrAnalysisTask
+    config:
+      atools.deltaSkyCorr: DeltaSkyCorrXYPlot
+      python: |
+        from lsst.analysis.tools.atools import *

python/lsst/analysis/tools/actions/scalar/scalarActions.py

@@ -246,3 +246,45 @@ 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 to median")
+    midKey = Field[str]("Key of bin midpoints Vector to median")
+
+    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) -> Scalar:
+        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

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 *

python/lsst/analysis/tools/atools/deltaSkyCorr.py

@@ -0,0 +1,58 @@
+# 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.vector import LoadVector, ConstantValue
+from ..interfaces import AnalysisTool
+from ..actions.scalar.scalarActions import MedianHistAction
+
+
+class DeltaSkyCorrXYPlot(AnalysisTool):
+    parametrizeBands = 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.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.metric.units = {
+            "median": "nJy",
+        }
+
+        self.produce.metric.newNames = {
+            "median": "delta_skyCorr_median",
+        }

python/lsst/analysis/tools/tasks/deltaSkyCorrAnalysis.py

@@ -23,8 +23,11 @@
 import numpy as np
 from lsst.pex.config import Field, ListField
 from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections
+from lsst.pipe.base import connectionTypes as ct
 from lsst.pipe.base.connectionTypes import Input, Output
 
+from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
+
 log = logging.getLogger(__name__)
 
 
@@ -55,6 +58,14 @@ class DeltaSkyCorrHistConnections(PipelineTaskConnections, dimensions=("instrume
         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",
@@ -69,11 +80,11 @@ class DeltaSkyCorrHistConfig(PipelineTaskConfig, pipelineConnections=DeltaSkyCor
     """Config class for DeltaSkyCorrHistTask."""
 
     bin_range = ListField[float](
-        doc="The lower and upper range for the histogram bins.",
+        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.",
+        doc="The width of each histogram bin, in nJy.",
         default=0.0001,
     )
 
@@ -96,7 +107,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
         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):
+    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).
@@ -118,6 +129,8 @@ def run(self, skyCorrs, injected_skyCorrs, num_initial_bgs):
             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
         -------
@@ -142,6 +155,7 @@ def run(self, skyCorrs, injected_skyCorrs, num_initial_bgs):
         # 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(
@@ -157,6 +171,8 @@ def run(self, skyCorrs, injected_skyCorrs, num_initial_bgs):
             # 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()
@@ -166,6 +182,7 @@ def run(self, skyCorrs, injected_skyCorrs, 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)
@@ -179,3 +196,26 @@ def run(self, skyCorrs, injected_skyCorrs, num_initial_bgs):
             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 = ct.Input(
+        doc="tmp",
+        name="delta_skyCorr_hist",
+        storageClass="ArrowNumpyDict",
+        deferLoad=True,
+        dimensions=("instrument", "visit"),
+    )
+
+
+class DeltaSkyCorrAnalysisConfig(AnalysisBaseConfig, pipelineConnections=DeltaSkyCorrAnalysisConnections):
+    pass
+
+
+class DeltaSkyCorrAnalysisTask(AnalysisPipelineTask):
+    ConfigClass = DeltaSkyCorrAnalysisConfig
+    _DefaultName = "deltaSkyCorrAnalysis"