From 99a193b6de64816e9832fbc4038baf315bb9cf2b Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Fri, 18 Jul 2025 09:20:59 -0700 Subject: [PATCH 1/2] compute fiducial magnitude limit --- .../pipe/tasks/computeExposureSummaryStats.py | 120 +++++++++++++++++- tests/test_computeExposureSummaryStats.py | 51 +++++++- 2 files changed, 168 insertions(+), 3 deletions(-) diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py index 687b74e87..0f6f58496 100644 --- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py +++ b/python/lsst/pipe/tasks/computeExposureSummaryStats.py @@ -127,6 +127,13 @@ 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}, + ) maxEffectiveTransparency = pexConfig.Field( dtype=float, doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).", @@ -159,6 +166,52 @@ def setDefaults(self): self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux" self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr" + def fiducialMagnitudeLimit(self, band, exposureTime, 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 + exposureTime : `float` + The exposure time [seconds] + 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) + magLimSnr = self.magLimSnr + + # Derived fiducial quantities + fiducialPsfArea = psf_sigma_to_psf_area(fiducialPsfSigma, pixelScale) + fiducialZeroPoint += 2.5*np.log10(exposureTime) + fiducialSkyBg = fiducialSkyBackground * exposureTime + 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. @@ -674,6 +727,17 @@ def update_effective_time_stats(self, summary, exposure): # Effective exposure time (seconds) effectiveTime = t_eff * exposureTime + # Effective exposure time (seconds) + # md = exposure.getMetadata() + # if md.get("LSST ISR UNITS", "adu") == "electron": + # gain = 1.0 + # else: + # gainList = list(ipIsr.getExposureGains(exposure).values()) + # gain = np.nanmean(gainList) + # + # fiducialMagLim = self.config.fiducialMagLim(band, exposureTime, summary.pixelScale, gain) + # effectiveTime = compute_effective_time(summary.magLim, fiducialMagLim, exposureTime) + # Output quantities summary.effTime = float(effectiveTime) summary.effTimePsfSigmaScale = float(f_eff) @@ -953,7 +1017,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 +1036,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, + exposureTime +): + """Compute the effective exposure time from m5 following the prescription described in SMTN-296. + + teff = 10**(0.8 * (maglim - maglim_fid) ) * expTime + + where magLim is the magnitude limit, magLim_fid is the fiducial magnitude limit calculated from + the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and expTime + is the exposure time (s). + + Parameters + ---------- + magLim : `float` + The measured magnitude limit [mag]. + fiducialMagLim : `float` + The fiducial magnitude limit [mag]. + exposureTime : `float` + Exposure time [s]. + + Returns + ------- + effectiveTime : `float` + The effective exposure time. + """ + effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * exposureTime + return effectiveTime diff --git a/tests/test_computeExposureSummaryStats.py b/tests/test_computeExposureSummaryStats.py index e9015d114..dc2abf337 100644 --- a/tests/test_computeExposureSummaryStats.py +++ b/tests/test_computeExposureSummaryStats.py @@ -34,6 +34,7 @@ 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 # Values from syseng_throughputs notebook assuming 30s exposure @@ -192,12 +193,16 @@ def testComputeMagnitudeLimit(self): gain = 1.0 # For testing non-unity gain gain_test = 2.0 + exposure_time = 30 + pixel_scale = 0.2 + # Output magnitude limit from syseng_throughputs notebook + m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10} 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 @@ -221,6 +226,48 @@ 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, + } + # Other properties + exposureTime = 30 # second + 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, exposureTime, 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", exposureTime, pixelScale, gain) + self.assertFloatsAlmostEqual(m5fid, nan, ignoreNaNs=True) + class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase): pass From 708448d41e66a30cb350f1b931670851824f5608 Mon Sep 17 00:00:00 2001 From: Alex Drlica-Wagner Date: Tue, 10 Mar 2026 15:32:22 -0700 Subject: [PATCH 2/2] update fiducial and effective time calculation --- .../pipe/tasks/computeExposureSummaryStats.py | 67 +++++++++++-------- tests/test_computeExposureSummaryStats.py | 32 ++++++--- 2 files changed, 62 insertions(+), 37 deletions(-) diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py index 0f6f58496..6e6d06bea 100644 --- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py +++ b/python/lsst/pipe/tasks/computeExposureSummaryStats.py @@ -134,6 +134,20 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config): "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).", @@ -166,7 +180,7 @@ def setDefaults(self): self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux" self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr" - def fiducialMagnitudeLimit(self, band, exposureTime, pixelScale, gain): + 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. @@ -174,8 +188,6 @@ def fiducialMagnitudeLimit(self, band, exposureTime, pixelScale, gain): ---------- band : `str` The band of interest - exposureTime : `float` - The exposure time [seconds] pixelScale : `float` The pixel scale [arcsec/pix] gain : `float` @@ -194,12 +206,13 @@ def fiducialMagnitudeLimit(self, band, exposureTime, pixelScale, gain): 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(exposureTime) - fiducialSkyBg = fiducialSkyBackground * exposureTime + fiducialZeroPoint += 2.5*np.log10(fiducialExpTime) + fiducialSkyBg = fiducialSkyBackground * fiducialExpTime fiducialReadNoise /= gain # Calculate the fiducial magnitude limit @@ -721,22 +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 - - # Effective exposure time (seconds) - # md = exposure.getMetadata() - # if md.get("LSST ISR UNITS", "adu") == "electron": - # gain = 1.0 - # else: - # gainList = list(ipIsr.getExposureGains(exposure).values()) - # gain = np.nanmean(gainList) - # - # fiducialMagLim = self.config.fiducialMagLim(band, exposureTime, summary.pixelScale, gain) - # effectiveTime = compute_effective_time(summary.magLim, fiducialMagLim, 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) @@ -1065,15 +1076,15 @@ def psf_sigma_to_psf_area(psfSigma, pixelScale): def compute_effective_time( magLim, fiducialMagLim, - exposureTime + fiducialExpTime ): """Compute the effective exposure time from m5 following the prescription described in SMTN-296. - teff = 10**(0.8 * (maglim - maglim_fid) ) * expTime + teff = 10**(0.8 * (magLim - fiducialMagLim) ) * fiducialExpTime - where magLim is the magnitude limit, magLim_fid is the fiducial magnitude limit calculated from - the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and expTime - is the exposure time (s). + 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 ---------- @@ -1081,13 +1092,13 @@ def compute_effective_time( The measured magnitude limit [mag]. fiducialMagLim : `float` The fiducial magnitude limit [mag]. - exposureTime : `float` - Exposure time [s]. + fiducialExpTime : `float` + The fiducial exposure time [s]. Returns ------- effectiveTime : `float` The effective exposure time. """ - effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * exposureTime + effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * fiducialExpTime return effectiveTime diff --git a/tests/test_computeExposureSummaryStats.py b/tests/test_computeExposureSummaryStats.py index dc2abf337..19cab69ab 100644 --- a/tests/test_computeExposureSummaryStats.py +++ b/tests/test_computeExposureSummaryStats.py @@ -36,6 +36,7 @@ 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- @@ -126,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) @@ -181,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 @@ -195,8 +202,6 @@ def testComputeMagnitudeLimit(self): gain_test = 2.0 exposure_time = 30 pixel_scale = 0.2 - # Output magnitude limit from syseng_throughputs notebook - m5_fid = {'g': 24.90, 'r': 24.48, 'i': 24.10} for band in ['g', 'r', 'i']: # Translate into DM quantities @@ -215,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) @@ -246,8 +257,11 @@ def testComputeFiducialMagnitudeLimit(self): 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 - exposureTime = 30 # second pixelScale = 0.2 # arcsec/pix gain = 1.0 # [e-/ADU] @@ -259,13 +273,13 @@ def testComputeFiducialMagnitudeLimit(self): # Compare fiducial m5 calculated from config to SMTN-002 values for band in fiducialMagLim.keys(): print(f"\n{band}") - m5fid = config.fiducialMagnitudeLimit(band, exposureTime, pixelScale, gain) + 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", exposureTime, pixelScale, gain) + m5fid = config.fiducialMagnitudeLimit("block", pixelScale, gain) self.assertFloatsAlmostEqual(m5fid, nan, ignoreNaNs=True)