Index: trunk/ppSub/src/ppSubReadout.c
===================================================================
--- trunk/ppSub/src/ppSubReadout.c	(revision 17798)
+++ trunk/ppSub/src/ppSubReadout.c	(revision 17813)
@@ -184,4 +184,47 @@
 #endif
 
+    // "Subtract" the mask and weight map
+    outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
+    if (inConv->weight && refConv->weight) {
+        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
+    }
+    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+
+    // Photometry is to be performed in two stages:
+    // 1. Measure the PSF using the PSF-matched images
+    // 2. Find and measure sources on the subtracted image
+
+    // Photometry stage 1: measure the PSF
+    if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
+        // We use a summed image as the basis for the PSF: this will have the maximum S/N.
+        outRO->image = (psImage*)psBinaryOp(outRO->image, inConv->image, "+", refConv->image);
+
+        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
+        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
+
+        // We have a list of sources, so we could use those to redetermine the PSF model.
+        // Until Gene makes the necessary adaptations to psphot, we will simply redetermine the PSF model from
+        // scratch.
+
+        // Need to ensure aperture residual is not calculated
+        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
+        if (!recipe) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find %s recipe.", PSPHOT_RECIPE);
+            return false;
+        }
+        const char *breakpoint = psMetadataLookupStr(NULL, recipe, "BREAK_POINT"); // Current break point
+        psMetadataAddStr(recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
+                         "ALTERED break point for psphot operations", "PSFMODEL");
+
+        if (!psphotReadout(config, view)) {
+            psWarning("Unable to perform photometry on subtracted image.");
+            psErrorStackPrint(stderr, "Error stack from photometry:");
+            psErrorClear();
+        }
+
+        psMetadataAddStr(recipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
+                         "RESTORED break point for psphot operations", breakpoint);
+    }
+
     // Do the subtraction
     {
@@ -197,9 +240,4 @@
 
         outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-        outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
-        if (minuend->weight && subtrahend->weight) {
-            outRO->weight = (psImage*)psBinaryOp(outRO->weight, minuend->weight, "+", subtrahend->weight);
-        }
-        outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
 
 #ifdef TESTING
@@ -249,7 +287,8 @@
     }
 
+    // Photometry stage 2: find and measure sources on the subtracted image
     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
+        // The PSF should already be stored for the readout
         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
 
         // Need to ensure aperture residual is not calculated