diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py index 872ffcf0e..1880b0ed4 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. @@ -661,6 +714,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) @@ -940,7 +1004,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. @@ -959,3 +1022,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 8d10e6c49..894abd314 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 @@ -176,14 +177,16 @@ def testComputeMagnitudeLimit(self): # Assumed values from SMTN-002 snr = 5 gain = 1.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 @@ -201,6 +204,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