From afc07ebc42ebdfb1a210a34f95b92733aa228fbd Mon Sep 17 00:00:00 2001 From: Pierre-Francois Leget Date: Tue, 1 Oct 2024 13:13:00 -0700 Subject: [PATCH] Add Ex metric for TEx computation in visit and ccd table. --- .../pipe/tasks/computeExposureSummaryStats.py | 108 +++++++++++++++++- python/lsst/pipe/tasks/postprocess.py | 12 +- tests/test_computeExposureSummaryStats.py | 16 +++ 3 files changed, 131 insertions(+), 5 deletions(-) diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py index ed9351434..e064ec974 100644 --- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py +++ b/python/lsst/pipe/tasks/computeExposureSummaryStats.py @@ -35,6 +35,7 @@ import lsst.afw.image as afwImage import lsst.geom as geom from lsst.meas.algorithms import ScienceSourceSelectorTask +from lsst.meas.algorithms.computeExPsf import ComputeExPsfTask, ComputeExPsfConfig from lsst.utils.timer import timeMethod import lsst.ip.isr as ipIsr @@ -137,6 +138,22 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config): doc="Signal-to-noise ratio for computing the magnitude limit depth.", default=5.0 ) + psfTE1TreecorrConfig = pexConfig.ConfigField( + dtype=ComputeExPsfConfig, + doc="Treecorr config for computing scalar value of TE1.", + ) + psfTE2TreecorrConfig = pexConfig.ConfigField( + dtype=ComputeExPsfConfig, + doc="Treecorr config for computing scalar value of TE1.", + ) + psfTE3TreecorrConfig = pexConfig.ConfigField( + dtype=ComputeExPsfConfig, + doc="Treecorr config for computing scalar value of TE1.", + ) + psfTE4TreecorrConfig = pexConfig.ConfigField( + dtype=ComputeExPsfConfig, + doc="Treecorr config for computing scalar value of TE1.", + ) def setDefaults(self): super().setDefaults() @@ -159,6 +176,22 @@ def setDefaults(self): self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux" self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr" + min_theta = [1e-6, 5.0, 1e-6, 5.0] + max_theta = [1.0, 100.0, 5.0, 20.0] + TExTreecorrConfig = [ + self.psfTE1TreecorrConfig, + self.psfTE2TreecorrConfig, + self.psfTE3TreecorrConfig, + self.psfTE4TreecorrConfig, + ] + + for texc, mint, maxt in zip(TExTreecorrConfig, min_theta, max_theta): + texc.treecorr.min_sep = mint / 60.0 + texc.treecorr.max_sep = maxt / 60.0 + texc.treecorr.nbins = 1 + texc.treecorr.bin_type = "Linear" + texc.treecorr.sep_units = "degree" + class ComputeExposureSummaryStatsTask(pipeBase.Task): """Task to compute exposure summary statistics. @@ -199,6 +232,18 @@ class ComputeExposureSummaryStatsTask(pipeBase.Task): - maxDistToNearestPsf - psfTraceRadiusDelta - psfApFluxDelta + - psfTE1E1 + - psfTE1E2 + - psfTE1Ex + - psfTE2E1 + - psfTE2E2 + - psfTE2Ex + - psfTE3E1 + - psfTE3E2 + - psfTE3Ex + - psfTE4E1 + - psfTE4E2 + - psfTE4Ex This quantity is computed based on the aperture correction map, the psfSigma, and the image mask to assess the robustness of the aperture @@ -325,6 +370,18 @@ def update_psf_stats( summary.psfTraceRadiusDelta = nan summary.psfApFluxDelta = nan summary.psfApCorrSigmaScaledDelta = nan + summary.psfTE1E1 = nan + summary.psfTE1E2 = nan + summary.psfTE1Ex = nan + summary.psfTE2E1 = nan + summary.psfTE2E2 = nan + summary.psfTE2Ex = nan + summary.psfTE3E1 = nan + summary.psfTE3E2 = nan + summary.psfTE3Ex = nan + summary.psfTE4E1 = nan + summary.psfTE4E2 = nan + summary.psfTE4Ex = nan if psf is None: return @@ -410,10 +467,13 @@ def update_psf_stats( psfE1 = (psfXX - psfYY)/(psfXX + psfYY) psfE2 = 2*psfXY/(psfXX + psfYY) - psfStarDeltaE1Median = np.median(starE1 - psfE1) - psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale='normal') - psfStarDeltaE2Median = np.median(starE2 - psfE2) - psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale='normal') + e1Residuals = starE1 - psfE1 + e2Residuals = starE2 - psfE2 + + psfStarDeltaE1Median = np.median(e1Residuals) + psfStarDeltaE1Scatter = sigmaMad(e1Residuals, scale='normal') + psfStarDeltaE2Median = np.median(e2Residuals) + psfStarDeltaE2Scatter = sigmaMad(e2Residuals, scale='normal') psfStarDeltaSizeMedian = np.median(starSize - psfSize) psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal') @@ -427,6 +487,46 @@ def update_psf_stats( summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter) summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter) + ra = psf_cat["coord_ra"] + dec = psf_cat["coord_dec"] + + # Comp TEx + + TExTreecorrConfig = [ + self.config.psfTE1TreecorrConfig, + self.config.psfTE2TreecorrConfig, + self.config.psfTE3TreecorrConfig, + self.config.psfTE4TreecorrConfig, + ] + + TExOutput = [ + [summary.psfTE1E1, summary.psfTE1E2, summary.psfTE1Ex], + [summary.psfTE2E1, summary.psfTE2E2, summary.psfTE2Ex], + [summary.psfTE3E1, summary.psfTE3E2, summary.psfTE3Ex], + [summary.psfTE4E1, summary.psfTE4E2, summary.psfTE4Ex], + ] + + isNotNan = np.array([True] * len(ra)) + isNotNan &= np.isfinite(ra) + isNotNan &= np.isfinite(dec) + isNotNan &= np.isfinite(e1Residuals) + isNotNan &= np.isfinite(e2Residuals) + + if np.sum(isNotNan) > 1: + + for config, texoutput in zip(TExTreecorrConfig, TExOutput): + + task = ComputeExPsfTask(config) + output = task.run( + e1Residuals[isNotNan], e2Residuals[isNotNan], + ra[isNotNan], dec[isNotNan], + units="degree", + ) + + texoutput[0] = output.metric_E1 + texoutput[1] = output.metric_E2 + texoutput[2] = output.metric_Ex + if image_mask is not None: maxDistToNearestPsf = maximum_nearest_psf_distance( image_mask, diff --git a/python/lsst/pipe/tasks/postprocess.py b/python/lsst/pipe/tasks/postprocess.py index 2418056fa..5bce6bd88 100644 --- a/python/lsst/pipe/tasks/postprocess.py +++ b/python/lsst/pipe/tasks/postprocess.py @@ -1230,7 +1230,10 @@ def run(self, visitSummaryRefs): "maxDistToNearestPsf", "effTime", "effTimePsfSigmaScale", "effTimeSkyBgScale", "effTimeZeroPointScale", - "magLim"] + "magLim", "psfTE1E1", "psfTE1E2", "psfTE1Ex", + "psfTE2E1", "psfTE2E2", "psfTE2Ex", + "psfTE3E1", "psfTE3E2", "psfTE3Ex", + "psfTE4E1", "psfTE4E2", "psfTE4Ex"] ccdEntry = summaryTable[selectColumns].to_pandas().set_index("id") # 'visit' is the human readable visit number. # 'visitId' is the key to the visitId table. They are the same. @@ -1364,6 +1367,13 @@ def run(self, visitSummaries): visitEntry["obsStart"] = visitEntry["expMidpt"] - 0.5 * pd.Timedelta(seconds=expTime) expTime_days = expTime / (60*60*24) visitEntry["obsStartMJD"] = visitEntry["expMidptMJD"] - 0.5 * expTime_days + TEx = [ + "psfTE1E1", "psfTE1E2", "psfTE1Ex", + "psfTE2E1", "psfTE2E2", "psfTE2Ex", + "psfTE3E1", "psfTE3E2", "psfTE3Ex", + ] + for tex in TEx: + visitEntry[tex] = visitRow[tex] visitEntries.append(visitEntry) # TODO: DM-30623, Add programId, exposureType, cameraTemp, diff --git a/tests/test_computeExposureSummaryStats.py b/tests/test_computeExposureSummaryStats.py index b2debe300..b0109fa27 100644 --- a/tests/test_computeExposureSummaryStats.py +++ b/tests/test_computeExposureSummaryStats.py @@ -162,6 +162,22 @@ def testComputeExposureSummary(self): self.assertFloatsAlmostEqual(summary.effTime, expTime, rtol=1e-3) self.assertFloatsAlmostEqual(summary.magLim, 26.584, rtol=1e-3) + # Check TExEx does exist, would need a more complex unit + # test with ~O(100) stars to check. + nan = float('nan') + self.assertFloatsAlmostEqual(summary.psfTE1E1, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE1E2, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE1Ex, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE2E1, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE2E2, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE2Ex, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE3E1, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE3E2, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE3Ex, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE4E1, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE4E2, nan, ignoreNaNs=True) + self.assertFloatsAlmostEqual(summary.psfTE4Ex, nan, ignoreNaNs=True) + def testComputeMagnitudeLimit(self): """Test the magnitude limit calculation using fiducials from SMTN-002 and syseng_throughputs."""