DM-46582: Add Ex metric for TEx computation in visit / ccd table

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,

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,

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."""