DM-51841: Compute effective time from magnitude limit

python/lsst/pipe/tasks/computeExposureSummaryStats.py

@@ -127,6 +127,27 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
         "Keyed by band.",
         default={'u': 25.0, 'g': 25.0, 'r': 25.0, 'i': 25.0, 'z': 25.0, 'y': 25.0},
     )
+    fiducialReadNoise = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial readnoise (electrons) assumed when calculating effective exposure time. "
+        "Keyed by band.",
+        default={'u': 9.0, 'g': 9.0, 'r': 9.0, 'i': 9.0, 'z': 9.0, 'y': 9.0},
+    )
+    fiducialExpTime = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial exposure time (seconds). "
+        "Keyed by band.",
+        default={'u': 30.0, 'g': 30.0, 'r': 30.0, 'i': 30.0, 'z': 30.0, 'y': 30.0},
+    )
+    fiducialMagLim = pexConfig.DictField(
+        keytype=str,
+        itemtype=float,
+        doc="Fiducial magnitude limit depth at SNR=5. "
+        "Keyed by band.",
+        default={'u': 25.0, 'g': 25.0, 'r': 25.0, 'i': 25.0, 'z': 25.0, 'y': 25.0},
+    )
     maxEffectiveTransparency = pexConfig.Field(
         dtype=float,
         doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).",
@@ -159,6 +180,51 @@ def setDefaults(self):
         self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux"
         self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr"
 
+    def fiducialMagnitudeLimit(self, band, pixelScale, gain):
+        """Compute the fiducial point-source magnitude limit based on config values.
+        This follows the conventions laid out in SMTN-002, LSE-40, and DMTN-296.
+
+        Parameters
+        ----------
+        band : `str`
+            The band of interest
+        pixelScale : `float`
+            The pixel scale [arcsec/pix]
+        gain : `float`
+            The instrumental gain for the exposure [e-/ADU]. The gain should
+            be 1.0 if the image units are [e-].
+
+        Returns
+        -------
+        magnitude_limit : `float`
+            The fiducial magnitude limit calculated from fiducial values.
+        """
+        nan = float("nan")
+
+        # Fiducials from config
+        fiducialPsfSigma = self.fiducialPsfSigma.get(band, nan)
+        fiducialSkyBackground = self.fiducialSkyBackground.get(band, nan)
+        fiducialZeroPoint = self.fiducialZeroPoint.get(band, nan)
+        fiducialReadNoise = self.fiducialReadNoise.get(band, nan)
+        fiducialExpTime = self.fiducialExpTime.get(band, nan)
+        magLimSnr = self.magLimSnr
+
+        # Derived fiducial quantities
+        fiducialPsfArea = psf_sigma_to_psf_area(fiducialPsfSigma, pixelScale)
+        fiducialZeroPoint += 2.5*np.log10(fiducialExpTime)
+        fiducialSkyBg = fiducialSkyBackground * fiducialExpTime
+        fiducialReadNoise /= gain
+
+        # Calculate the fiducial magnitude limit
+        fiducialMagLim = compute_magnitude_limit(fiducialPsfArea,
+                                                 fiducialSkyBg,
+                                                 fiducialZeroPoint,
+                                                 fiducialReadNoise,
+                                                 gain,
+                                                 magLimSnr)
+
+        return fiducialMagLim
+
 
 class ComputeExposureSummaryStatsTask(pipeBase.Task):
     """Task to compute exposure summary statistics.
@@ -668,11 +734,20 @@ def update_effective_time_stats(self, summary, exposure):
             fiducialSkyBackground = self.config.fiducialSkyBackground[band]
             b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime)
 
-        # Effective exposure time scale factor
-        t_eff = f_eff * c_eff * b_eff
+        # Old effective exposure time scale factor
+        # t_eff = f_eff * c_eff * b_eff
 
-        # Effective exposure time (seconds)
-        effectiveTime = t_eff * exposureTime
+        # New effective exposure time (seconds) from magnitude limit
+        if band not in self.config.fiducialMagLim:
+            self.log.debug(f"Fiducial magnitude limit not found for {band}")
+            effectiveTime = nan
+        elif band not in self.config.fiducialExpTime:
+            self.log.debug(f"Fiducial exposure time not found for {band}")
+            effectiveTime = nan
+        else:
+            effectiveTime = compute_effective_time(summary.magLim,
+                                                   self.config.fiducialMagLim[band],
+                                                   self.config.fiducialExpTime[band])
 
         # Output quantities
         summary.effTime = float(effectiveTime)
@@ -953,7 +1028,6 @@ def compute_magnitude_limit(
     -------
     magnitude_limit : `float`
         The expected magnitude limit at the given signal to noise.
-
     """
     # Effective number of pixels within the PSF calculated directly
     # from the PSF model.
@@ -973,3 +1047,58 @@ def compute_magnitude_limit(
     magnitude_limit = -2.5*np.log10(source) + zeroPoint
 
     return magnitude_limit
+
+
+def psf_sigma_to_psf_area(psfSigma, pixelScale):
+    """Convert psfSigma [pixels] to an approximate psfArea [pixel^2] as defined in LSE-40.
+    This is the same implementation followed by SMTN-002 to estimate SNR=5 magnitude limits.
+
+    Parameters
+    ----------
+    psfSigma : `float`
+        The PSF sigma value [pix].
+    pixelScale : `float`
+        The pixel scale [arcsec/pix].
+
+    Returns
+    -------
+    psf_area : `float`
+        Approximation of the PSF area [pix^2]
+    """
+    # Follow SMTN-002 to convert to geometric and effective FWHM
+    fwhm_geom = psfSigma * 2.355 * pixelScale
+    fwhm_eff = (fwhm_geom - 0.052) / 0.822
+    # Area prefactor comes from LSE-40
+    psf_area = 2.266 * (fwhm_eff / pixelScale)**2
+    return psf_area
+
+
+def compute_effective_time(
+        magLim,
+        fiducialMagLim,
+        fiducialExpTime
+):
+    """Compute the effective exposure time from m5 following the prescription described in SMTN-296.
+
+      teff = 10**(0.8 * (magLim - fiducialMagLim) ) * fiducialExpTime
+
+    where `magLim` is the magnitude limit, `fiducialMagLim` is the fiducial magnitude limit calculated from
+    the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and
+    `fiducialExpTime` is the fiducial exposure time (s).
+
+    Parameters
+    ----------
+    magLim : `float`
+        The measured magnitude limit [mag].
+    fiducialMagLim : `float`
+        The fiducial magnitude limit [mag].
+    fiducialExpTime : `float`
+        The fiducial exposure time [s].
+
+    Returns
+    -------
+    effectiveTime : `float`
+        The effective exposure time.
+    """
+    effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * fiducialExpTime
+    return effectiveTime

tests/test_computeExposureSummaryStats.py

@@ -34,7 +34,9 @@
 from lsst.afw.geom import makeCdMatrix, makeSkyWcs
 from lsst.afw.cameraGeom.testUtils import DetectorWrapper
 from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsTask
+from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsConfig
 from lsst.pipe.tasks.computeExposureSummaryStats import compute_magnitude_limit
+from lsst.pipe.tasks.computeExposureSummaryStats import compute_effective_time
 
 # Values from syseng_throughputs notebook assuming 30s exposure
 # consisting of 2x15s snaps each with readnoise of 9e-
@@ -125,10 +127,16 @@ def testComputeExposureSummary(self):
 
         # Configure and run the task
         expSummaryTask = ComputeExposureSummaryStatsTask()
+        config = expSummaryTask.config
+
         # Configure nominal values for effective time calculation (normalized to 1s exposure)
-        expSummaryTask.config.fiducialZeroPoint = {band: float(zp - 2.5*np.log10(expTime))}
-        expSummaryTask.config.fiducialPsfSigma = {band: float(psfSize)}
-        expSummaryTask.config.fiducialSkyBackground = {band: float(skyMean/expTime)}
+        config.fiducialZeroPoint = {band: float(zp - 2.5*np.log10(expTime))}
+        config.fiducialPsfSigma = {band: float(psfSize)}
+        config.fiducialSkyBackground = {band: float(skyMean/expTime)}
+        config.fiducialReadNoise = {band: float(readNoise)}
+        config.fiducialExpTime = {band: expTime}
+        config.fiducialMagLim = {band: 26.584}
+
         # Run the task
         summary = expSummaryTask.run(exposure, None, background)
 
@@ -180,8 +188,8 @@ def testComputeExposureSummary(self):
         self.assertFloatsAlmostEqual(summary.zenithDistance, 30.57112, atol=1e-5)
 
         # Effective exposure time and depth
+        self.assertFloatsAlmostEqual(summary.magLim, config.fiducialMagLim[band], rtol=1e-3)
         self.assertFloatsAlmostEqual(summary.effTime, expTime, rtol=1e-3)
-        self.assertFloatsAlmostEqual(summary.magLim, 26.584, rtol=1e-3)
 
     def testComputeMagnitudeLimit(self):
         """Test the magnitude limit calculation using fiducials from SMTN-002
@@ -192,12 +200,14 @@ def testComputeMagnitudeLimit(self):
         gain = 1.0
         # For testing non-unity gain
         gain_test = 2.0
+        exposure_time = 30
+        pixel_scale = 0.2
 
         for band in ['g', 'r', 'i']:
             # Translate into DM quantities
-            psfArea = 2.266 * (fwhm_eff_fid[band] / 0.2)**2
+            psfArea = 2.266 * (fwhm_eff_fid[band] / pixel_scale)**2
             skyBg = skycounts_fid[band]
-            zeroPoint = zeropoint_fid[band] + 2.5*np.log10(30)
+            zeroPoint = zeropoint_fid[band] + 2.5*np.log10(exposure_time)
             readNoise = readnoise_fid[band]
 
             # Calculate the M5 values
@@ -210,6 +220,12 @@ def testComputeMagnitudeLimit(self):
                                               readNoise/gain_test, gain_test, snr)
             self.assertFloatsAlmostEqual(m5_gain, m5, atol=1e-6)
 
+            # Calculate effective time values
+            eff_time = compute_effective_time(m5, m5_fid[band], exposure_time)
+            self.assertFloatsAlmostEqual(eff_time, exposure_time, rtol=1e-2)
+            eff_time_half = compute_effective_time(m5 - 0.375, m5_fid[band], exposure_time)
+            self.assertFloatsAlmostEqual(eff_time_half, exposure_time/2.0, rtol=1e-2)
+
         # Check that input NaN lead to output NaN
         nan = float('nan')
         m5 = compute_magnitude_limit(nan, skyBg, zeroPoint, readNoise, gain, snr)
@@ -221,6 +237,51 @@ def testComputeMagnitudeLimit(self):
         m5 = compute_magnitude_limit(psfArea, skyBg, zeroPoint, nan, gain, snr)
         self.assertFloatsAlmostEqual(m5, nan, ignoreNaNs=True)
 
+    def testComputeFiducialMagnitudeLimit(self):
+        """Test the magnitude limit calculation using fiducials from SMTN-002
+        and syseng_throughputs."""
+
+        # Setup the config
+        config = ComputeExposureSummaryStatsConfig()
+        # Configure fiducial values.
+        # These come from obs_lsst/config/lsstCam
+        config.fiducialZeroPoint = {
+            "u": 26.52, "g": 28.51, "r": 28.36, "i": 28.17, "z": 27.78, "y": 26.82,
+        }
+        config.fiducialPsfSigma = {
+            "u": 1.72, "g": 1.63, "r": 1.56, "i": 1.51, "z": 1.47, "y": 1.44,
+        }
+        config.fiducialSkyBackground = {
+            "u": 1.51, "g": 15.45, "r": 32.95, "i": 52.94, "z": 79.40, "y": 90.57,
+        }
+        config.fiducialReadNoise = {
+            "u": 9.0, "g": 9.0, "r": 9.0, "i": 9.0, "z": 9.0, "y": 9.0,
+        }
+        config.fiducialExpTime = {
+            "u": 30.0, "g": 30.0, "r": 30.0, "i": 30.0, "z": 30.0, "y": 30.0,
+        }
+
+        # Other properties
+        pixelScale = 0.2  # arcsec/pix
+        gain = 1.0  # [e-/ADU]
+
+        # Fiducial m5 from SMTN-002
+        fiducialMagLim = {
+            "u": 23.70, "g": 24.97, "r": 24.52, "i": 24.13, "z": 23.56, "y": 22.55
+        }
+
+        # Compare fiducial m5 calculated from config to SMTN-002 values
+        for band in fiducialMagLim.keys():
+            print(f"\n{band}")
+            m5fid = config.fiducialMagnitudeLimit(band, pixelScale, gain)
+            print(f"{fiducialMagLim[band]:.2f} {m5fid:.3f}")
+            self.assertFloatsAlmostEqual(m5fid, fiducialMagLim[band], atol=1e-2)
+
+        # Check that passing an unknown band outputs NaN
+        nan = float('nan')
+        m5fid = config.fiducialMagnitudeLimit("block", pixelScale, gain)
+        self.assertFloatsAlmostEqual(m5fid, nan, ignoreNaNs=True)
+
 
 class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
     pass