Merge pull request #281 from lsst/tickets/DM-47425

DM-47425: Add task that can use arbitrary inputs (e.g. on-sky data) to fix up amp gain ratios.

pipelines/LATISS/cpPtcFixupGainRatios.yaml

@@ -0,0 +1,4 @@
+description: cp_pipe LATISS ptc fixup.
+instrument: lsst.obs.lsst.Latiss
+imports:
+  - location: $CP_PIPE_DIR/pipelines/_ingredients/cpPtcFixupGainRatios.yaml

pipelines/LATISS/cpPtcRename.yaml

@@ -0,0 +1,4 @@
+description: cp_pipe ptc renaming (LATISS).
+instrument: lsst.obs.lsst.Latiss
+imports:
+  - location: $CP_PIPE_DIR/pipelines/_ingredients/cpPtcRename.yaml

pipelines/LSSTComCam/cpPtcFixupGainRatios.yaml

@@ -0,0 +1,4 @@
+description: cp_pipe LSSTComCam ptc fixup.
+instrument: lsst.obs.lsst.LsstComCam
+imports:
+  - location: $CP_PIPE_DIR/pipelines/_ingredients/cpPtcFixupGainRatios.yaml

pipelines/LSSTComCam/cpPtcRename.yaml

@@ -0,0 +1,4 @@
+description: cp_pipe ptc renaming (ComCam).
+instrument: lsst.obs.lsst.LsstComCam
+imports:
+  - location: $CP_PIPE_DIR/pipelines/_ingredients/cpPtcRename.yaml

pipelines/_ingredients/cpPtcFixupGainRatios.yaml

@@ -0,0 +1,28 @@
+description: cp_pipe ptc fixup.
+tasks:
+  cpPtcFixupGainRatiosIsr:
+    class: lsst.ip.isr.IsrTaskLSST
+    config:
+      connections.ccdExposure: "raw"
+      connections.outputExposure: "cpPtcFixupGainRatiosIsrExp"
+      python: |
+        from lsst.cp.pipe import configureIsrTaskLSSTForCalibrations
+
+        configureIsrTaskLSSTForCalibrations(config)
+
+        config.doBootstrap = True
+        config.doCrosstalk = True
+        config.crosstalk.doQuadraticCrosstalkCorrection = False
+        config.doLinearize = True
+        config.doDefect = True
+        config.doAmpOffset = True
+        config.ampOffset.ampEdgeMaxOffset = 100000.0
+        config.ampOffset.ampEdgeInset = 10
+        config.ampOffset.doBackground = False
+        config.ampOffset.doDetection = False
+        config.ampOffset.doApplyAmpOffset = False
+  cpPtcFixupGainRatios:
+    class: lsst.cp.pipe.ptc.PhotonTransferCurveFixupGainRatiosTask
+
+contracts:
+  - cpPtcFixupGainRatiosIsr.doBootstrap == True

pipelines/_ingredients/cpPtcRename.yaml

@@ -0,0 +1,4 @@
+description: cp_pipe ptc renaming.
+tasks:
+  cpPtcRename:
+    class: lsst.cp.pipe.PhotonTransferCurveRenameTask

python/lsst/cp/pipe/ptc/__init__.py

@@ -24,3 +24,4 @@
 from .cpPlotPtc import *
 from .cpPtcExtract import *
 from .cpPtcSolve import *
+from .cpPtcFixupGainRatios import *

python/lsst/cp/pipe/ptc/cpPtcFixupGainRatios.py

@@ -0,0 +1,222 @@
+# This file is part of cp_pipe.
+#
+# 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 copy
+import numpy as np
+
+from lsst.ip.isr import PhotonTransferCurveDataset
+import lsst.pex.config as pexConfig
+import lsst.pipe.base as pipeBase
+import lsst.pipe.base.connectionTypes as cT
+
+from ..utils import ampOffsetGainRatioFixup
+
+
+__all__ = [
+    "PhotonTransferCurveFixupGainRatiosConfig",
+    "PhotonTransferCurveFixupGainRatiosTask",
+    "PhotonTransferCurveRenameConfig",
+    "PhotonTransferCurveRenameTask",
+]
+
+
+class PhotonTransferCurveFixupGainRatiosConnections(
+    pipeBase.PipelineTaskConnections,
+    dimensions=("instrument", "detector")
+):
+    exposureMetadata = cT.Input(
+        name="cpPtcFixupGainRatiosIsrExp.metadata",
+        doc="Input exposures for gain ratio fixup.",
+        storageClass="PropertyList",
+        dimensions=("instrument", "exposure", "detector"),
+        multiple=True,
+    )
+    inputPtc = cT.PrerequisiteInput(
+        name="ptc",
+        doc="Input PTC to modify.",
+        storageClass="PhotonTransferCurveDataset",
+        dimensions=("instrument", "detector"),
+        isCalibration=True,
+    )
+    outputPtc = cT.Output(
+        name="ptcFixed",
+        doc="Output modified PTC.",
+        storageClass="PhotonTransferCurveDataset",
+        dimensions=("instrument", "detector"),
+        multiple=False,
+        isCalibration=True,
+    )
+
+
+class PhotonTransferCurveFixupGainRatiosConfig(
+    pipeBase.PipelineTaskConfig,
+    pipelineConnections=PhotonTransferCurveFixupGainRatiosConnections,
+):
+    ampOffsetGainRatioMinAdu = pexConfig.Field(
+        dtype=float,
+        doc="Minimum number of adu to use for amp offset gain ratio fixup.",
+        default=1000.0,
+    )
+    ampOffsetGainRatioMaxAdu = pexConfig.Field(
+        dtype=float,
+        doc="Maximum number of adu to use for amp offset gain ratio fixup.",
+        default=40000.0,
+    )
+
+
+class PhotonTransferCurveFixupGainRatiosTask(pipeBase.PipelineTask):
+    """Task to use on-sky amp ratios to fix up gain ratios in a PTC.
+
+    This uses the ampOffsetGainRatioFixup with on-sky data (preferably
+    twilight flats or similar) to update gain ratios.
+    """
+    ConfigClass = PhotonTransferCurveFixupGainRatiosConfig
+    _DefaultName = "cpPhotonTransferCurveFixupGainRatios"
+
+    def runQuantum(self, butlerQC, inputRefs, outputRefs):
+        # docstring inherited.
+        inputs = butlerQC.get(inputRefs)
+        outputs = self.run(inputPtc=inputs["inputPtc"], exposureMetadata=inputs["exposureMetadata"])
+        butlerQC.put(outputs, outputRefs)
+
+    def run(self, *, inputPtc, exposureMetadata):
+        """Run the gain ratio fixup task.
+
+        Parameters
+        ----------
+        inputPtc : `lsst.ip.isr.PhotonTransferCurveDataset`
+            Input PTC to modify.
+        exposureMetadata: `list` [`lsst.daf.base.PropertyList`]
+            Input exposure metadata.
+
+        Returns
+        -------
+        results : `lsst.pipe.base.Struct`
+            The output struct contains:
+
+            ``outputPtc``
+                The output modified ptc.
+        """
+        ampNames = inputPtc.ampNames
+
+        # Create a set of fake partial PTC datasets.
+        fakePtc = PhotonTransferCurveDataset(
+            ampNames=ampNames,
+            ptcFitType="FAKEPTC",
+            covMatrixSide=1,
+            covMatrixSideFullCovFit=1,
+        )
+
+        for i, metadata in enumerate(exposureMetadata):
+            fakePartialPtc = PhotonTransferCurveDataset(ampNames=ampNames, ptcFitType="PARTIAL")
+
+            for ampName in ampNames:
+                fakePartialPtc.setAmpValuesPartialDataset(
+                    ampName,
+                    inputExpIdPair=(2*i, 2*i + 1),
+                    rawExpTime=float(i),
+                    rawMean=metadata[f"LSST ISR FINAL MEDIAN {ampName}"],
+                    rawVar=metadata[f"LSST ISR FINAL STDEV {ampName}"]**2.,
+                    ampOffset=metadata[f"LSST ISR AMPOFFSET PEDESTAL {ampName}"],
+                    expIdMask=True,
+                    gain=inputPtc.gainUnadjusted[ampName],
+                    noise=metadata[f"LSST ISR READNOISE {ampName}"]*inputPtc.gain[ampName],
+                    covariance=np.zeros((1, 1)),
+                    covSqrtWeights=np.zeros((1, 1)),
+                )
+
+            fakePtc.appendPartialPtc(fakePartialPtc)
+
+        detectorMeans = np.zeros(len(exposureMetadata))
+
+        for i in range(len(detectorMeans)):
+            arr = np.asarray([fakePtc.rawMeans[ampName][i] for ampName in ampNames])
+            detectorMeans[i] = np.nanmean(arr)
+
+        index = np.argsort(detectorMeans)
+        fakePtc.sort(index)
+
+        for ampName in ampNames:
+            fakePtc.finalMeans[ampName][:] = fakePtc.rawMeans[ampName].copy()
+            fakePtc.finalVars[ampName][:] = fakePtc.rawVars[ampName].copy()
+            fakePtc.gainUnadjusted[ampName] = inputPtc.gainUnadjusted[ampName]
+            fakePtc.gain[ampName] = inputPtc.gainUnadjusted[ampName]
+
+        ampOffsetGainRatioFixup(
+            fakePtc,
+            self.config.ampOffsetGainRatioMinAdu,
+            self.config.ampOffsetGainRatioMaxAdu,
+            log=self.log,
+        )
+
+        outputPtc = copy.copy(inputPtc)
+
+        # Replace the gain (leave gainUnadjusted alone).
+        for ampName in ampNames:
+            outputPtc.gain[ampName] = fakePtc.gain[ampName]
+
+        return pipeBase.Struct(
+            outputPtc=outputPtc,
+        )
+
+
+class PhotonTransferCurveRenameConnections(
+    pipeBase.PipelineTaskConnections,
+    dimensions=("instrument", "detector")
+):
+    inputPtc = cT.PrerequisiteInput(
+        name="ptcFixed",
+        doc="Input PTC to rename.",
+        storageClass="PhotonTransferCurveDataset",
+        dimensions=("instrument", "detector"),
+        isCalibration=True,
+    )
+    outputPtc = cT.Output(
+        name="ptc",
+        doc="Output PTC that has been renamed.",
+        storageClass="PhotonTransferCurveDataset",
+        dimensions=("instrument", "detector"),
+        multiple=False,
+        isCalibration=True,
+    )
+
+
+class PhotonTransferCurveRenameConfig(
+    pipeBase.PipelineTaskConfig,
+    pipelineConnections=PhotonTransferCurveRenameConnections,
+):
+    pass
+
+
+class PhotonTransferCurveRenameTask(pipeBase.PipelineTask):
+    """Task to rename a ptcFixed into a ptc."""
+    ConfigClass = PhotonTransferCurveRenameConfig
+    _DefaultName = "cpPhotonTransferCurveRename"
+
+    def runQuantum(self, butlerQC, inputRefs, outputRefs):
+        # docstring inherited.
+        inputs = butlerQC.get(inputRefs)
+
+        outputs = pipeBase.Struct(ptc=inputs["inputPtc"])
+        butlerQC.put(outputs, outputRefs)
+
+    def run(self):
+        pass

tests/test_pipelines.py

@@ -77,6 +77,8 @@ def _get_pipelines(self, exclude=[]):
             "cpDarkBootstrap.yaml",
             "cpFlatBootstrap.yaml",
             "cpSpectroFlat.yaml",
+            "cpPtcFixupGainRatios.yaml",
+            "cpPtcRename.yaml",
         }
 
         for ex in exclude:
@@ -157,6 +159,8 @@ def test_lsstcam_pipelines(self):
                 "cpBiasBootstrap.yaml",
                 "cpDarkBootstrap.yaml",
                 "cpFlatBootstrap.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "LSSTCam", pipeline))
 
@@ -170,6 +174,8 @@ def test_lsstcam_imsim_pipelines(self):
                 "cpBiasBootstrap.yaml",
                 "cpDarkBootstrap.yaml",
                 "cpFlatBootstrap.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "LSSTCam-imSim", pipeline))
 
@@ -204,6 +210,8 @@ def test_lsstcomcamsim_pipelines(self):
                 "cpLinearizer.yaml",
                 "cpCrosstalk.yaml",
                 "cpCti.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "LSSTComCamSim", pipeline))
 
@@ -216,6 +224,8 @@ def test_lsst_ts8_pipelines(self):
                 "cpBiasBootstrap.yaml",
                 "cpDarkBootstrap.yaml",
                 "cpFlatBootstrap.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "LSST-TS8", pipeline))
 
@@ -229,6 +239,8 @@ def test_decam_pipelines(self):
                 "cpBiasBootstrap.yaml",
                 "cpDarkBootstrap.yaml",
                 "cpFlatBootstrap.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "DECam", pipeline))
 
@@ -242,6 +254,8 @@ def test_hsc_pipelines(self):
                 "cpBiasBootstrap.yaml",
                 "cpDarkBootstrap.yaml",
                 "cpFlatBootstrap.yaml",
+                "cpPtcFixupGainRatios.yaml",
+                "cpPtcRename.yaml",
         ]):
             self._check_pipeline(os.path.join(self.pipeline_path, "HSC", pipeline))
 

tests/test_ptc.py

@@ -877,7 +877,8 @@ def test_ptcFitBootstrap(self):
             self.ptcFitAndCheckPtc(fitType=fitType, order=order, doFitBootstrap=True)
 
     def test_ampOffsetGainRatioFixup(self):
-        """Test the ampOffsetGainRatioFixup utility code."""
+        """Test the ampOffsetGainRatioFixup code via
+        PhotonTransferCurveFixupGainRatiosTask."""
         rng = np.random.RandomState(12345)
 
         gainsTruth = rng.normal(loc=1.7, scale=0.05, size=len(self.ampNames))
@@ -894,6 +895,12 @@ def test_ampOffsetGainRatioFixup(self):
         nFlat = 20
         meansElectron = gainsMedian * np.linspace(500.0, 30000.0, nFlat)
 
+        # Set up the fixup task.
+        fixupConfig = cpPipe.ptc.PhotonTransferCurveFixupGainRatiosConfig()
+        fixupConfig.ampOffsetGainRatioMinAdu = 1000.0
+        fixupConfig.ampOffsetGainRatioMaxAdu = 20000.0
+        fixupTask = cpPipe.ptc.PhotonTransferCurveFixupGainRatiosTask(config=fixupConfig)
+
         for testMode in ["full", "badamp"]:
             badAmp = None
 
@@ -928,17 +935,23 @@ def test_ampOffsetGainRatioFixup(self):
             ampOffset = AmpOffsetTask(config=config)
 
             detector = self.flatExp1.getDetector()
+            metadatas = []
             for i in range(nFlat):
                 exp = self.flatExp1.clone()
                 for amp in detector:
-                    exp[amp.getBBox()].image.array[:, :] = ptc.finalMeans[amp.getName()][i]
+                    amp_name = amp.getName()
+                    exp[amp.getBBox()].image.array[:, :] = ptc.finalMeans[amp_name][i]
+                    md = exp.metadata
 
-                res = ampOffset.run(exp)
+                    md[f"LSST ISR FINAL MEDIAN {amp_name}"] = ptc.finalMeans[amp_name][i]
+                    md[f"LSST ISR FINAL STDEV {amp_name}"] = np.sqrt(ptc.finalMeans[amp_name][i])
+                    md[f"LSST ISR READNOISE {amp_name}"] = 0.0
 
-                for j, amp in enumerate(detector):
-                    ptc.ampOffsets[amp.getName()][i] = res.pedestals[j]
+                ampOffset.run(exp)
+                metadatas.append(exp.metadata)
 
-            ampOffsetGainRatioFixup(ptc, 1000.0, 20000.0)
+            result = fixupTask.run(inputPtc=ptc, exposureMetadata=metadatas)
+            ptc = result.outputPtc
 
             # Check that the flats are flat after adjustment.
             for i in range(nFlat):