diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 12524628..62b14189 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -100,6 +100,12 @@ class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections,
         storageClass="MatchingKernel",
         name="{fakesType}{coaddName}Diff_psfMatchKernel",
     )
+    kernelSources = connectionTypes.Output(
+        doc="Final selection of sources used for psf matching.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="SourceCatalog",
+        name="{fakesType}{coaddName}Diff_psfMatchSources"
+    )
 
 
 class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
@@ -117,6 +123,12 @@ class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
         storageClass="MatchingKernel",
         name="{fakesType}{coaddName}Diff_psfScoreMatchKernel",
     )
+    kernelSources = connectionTypes.Output(
+        doc="Final selection of sources used for psf matching.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="SourceCatalog",
+        name="{fakesType}{coaddName}Diff_psfScoreMatchSources"
+    )
 
 
 class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
@@ -247,7 +259,7 @@ def setDefaults(self):
         self.makeKernel.kernel.active.spatialKernelOrder = 1
         self.makeKernel.kernel.active.spatialBgOrder = 2
         self.sourceSelector.doFlags = True
-        self.sourceSelector.flags.good = ["calib_psf_used",]
+        self.sourceSelector.flags.good = ["calib_psf_used", ]
         self.sourceSelector.doSignalToNoise = True  # apply signal to noise filter
         self.sourceSelector.signalToNoise.minimum = 10
         self.sourceSelector.signalToNoise.maximum = 500
@@ -529,7 +541,8 @@ def runConvolveTemplate(self, template, science, selectSources):
                                      matchedTemplate=matchedTemplate,
                                      matchedScience=science,
                                      backgroundModel=kernelResult.backgroundModel,
-                                     psfMatchingKernel=kernelResult.psfMatchingKernel)
+                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
+                                     kernelSources=kernelSources)
 
     def runConvolveScience(self, template, science, selectSources):
         """Convolve the science image with a PSF-matching kernel and subtract
@@ -595,7 +608,8 @@ def runConvolveScience(self, template, science, selectSources):
                                      matchedTemplate=matchedTemplate,
                                      matchedScience=matchedScience,
                                      backgroundModel=kernelResult.backgroundModel,
-                                     psfMatchingKernel=kernelResult.psfMatchingKernel,)
+                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
+                                     kernelSources=kernelSources)
 
     def finalize(self, template, science, difference, kernel,
                  templateMatched=True,
@@ -1127,7 +1141,8 @@ def runPreconvolve(self, template, science, matchedScience, selectSources, preCo
                                      matchedTemplate=matchedTemplate,
                                      matchedScience=matchedScience,
                                      backgroundModel=kernelResult.backgroundModel,
-                                     psfMatchingKernel=kernelResult.psfMatchingKernel)
+                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
+                                     kernelSources=kernelSources)
 
 
 def checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0.,