Index: /trunk/psphot/src/psphotStackImageLoop.c
===================================================================
--- /trunk/psphot/src/psphotStackImageLoop.c	(revision 36118)
+++ /trunk/psphot/src/psphotStackImageLoop.c	(revision 36119)
@@ -22,7 +22,5 @@
 
     pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
-    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
-    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
 
     if (!input) {
@@ -32,23 +30,7 @@
 
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    if (useRaw && inputRaw == NULL) {
-        psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
-        useRaw = false;
-    }
-
-    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
-    bool match_psfs = radial_apertures;
-    bool needConvolved = radial_apertures || !useRaw;
-    if (!needConvolved) {
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.CNV");
-    }
 
     // just load the full set of images up front except for EXPNUM which we defer
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
-    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
 
@@ -71,17 +53,8 @@
 		psMemDump("load");
 
-                if (match_psfs) {
-                    // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
-                    if (!psphotStackMatchPSFs (config, view)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
-                } else {
-                    if (!psphotStackAllocateOutput (config, view, recipe)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
+		if (!psphotStackAllocateOutput (config, view, recipe)) {
+		    psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+		    psFree (view);
+		    return false;
                 }
 		psMemDump("stackmatch");
@@ -126,5 +99,5 @@
 
     // Load the appropriate EXPNUM image
-    pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
+    pmFPAfileActivate (config->files, true, "PSPHOT.STACK.EXPNUM.RAW");
 
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
@@ -146,5 +119,5 @@
                 if (! readout->data_exists) { continue; }
 
-                if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
+                if (!psphotSetNFrames (config, view, input->name)) ESCAPE ("failed to setNFrames.");
             }
         }
@@ -167,7 +140,4 @@
 bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
 
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
 
@@ -179,7 +149,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = useRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
@@ -205,7 +173,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
@@ -242,7 +208,5 @@
 	psAssert (output, "missing file?");
 
-	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
 	psAssert (input, "missing input file");
 
Index: /trunk/psphot/src/psphotStackReadout.c
===================================================================
--- /trunk/psphot/src/psphotStackReadout.c	(revision 36118)
+++ /trunk/psphot/src/psphotStackReadout.c	(revision 36119)
@@ -1,16 +1,43 @@
 # include "psphotInternal.h"
 
+static bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF);
+static bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF, int index);
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
 static void logMemStats(const char *heading);
 
-// we have 3 possible real filesets:
+// relevant filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
-
-// we have 3 files on which we operate:
-// DET (detection image)       : nominally RAW (optionally CNV?)
-// SRC (source analysis image) : nominally CNV (optionally RAW)
-// OUT (psf-matched images)    : always OUT
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+
+// XXX STACK_OUT currently is a copy of STACK_RAW, but should be a pointer to is as in psphot (single)
+
+// TEST CODE, can be removed
+bool psphotDumpImages (pmConfig *config, const pmFPAview *view, const char *filerule, char *base) {
+
+    // XXX do nothing
+    return true;
+
+    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
+
+    for (int i = 0; i < num; i++) {
+	// find the currently selected readout
+	pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+	psAssert (file, "missing file?");
+
+	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+	psAssert (readout, "missing readout?");
+
+	char line[256];
+	snprintf (line, 256, "%s.%d.im.fits", base, i);
+	psphotSaveImage (NULL, readout->image, line);
+
+	snprintf (line, 256, "%s.%d.wt.fits", base, i);
+	psphotSaveImage (NULL, readout->variance, line);
+
+	snprintf (line, 256, "%s.%d.mk.fits", base, i);
+	psphotSaveImage (NULL, readout->mask, line);
+    }
+    return true;
+}
 
 bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -67,18 +94,12 @@
     psAssert (breakPt, "configuration error: set BREAK_POINT");
 
-    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
-    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
-    char *STACK_DET = STACK_RAW;
-
     // load WCS 
-    if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
-        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
+    if (!psphotStackLoadWCS(config, view, STACK_RAW)) {
+        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_RAW);
         return false;
     }
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
+    if (!psphotAddPhotcode (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -86,77 +107,49 @@
 
     // Generate the mask and weight images (if not supplied) and set mask bits. 
-    // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
-    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_DET)) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (strcmp(STACK_SRC, STACK_DET)) {
-#define MODEL_BACKGROUND_SRC 1
-#ifdef MODEL_BACKGROUND_SRC
-        // work around the fact that the background levels on the convolved
-        // and unconvolved stacks can be different
-        if (!psphotModelBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
-#endif
-        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
 #ifdef MAKE_CHISQ_IMAGE
     // also make the chisq detection image
-    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
+    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!strcasecmp (breakPt, "CHISQ")) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 #endif
-    if (!strcasecmp (breakPt, "CHISQ")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
-        // this only happens if we had an error in psphotFindDetections
+    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
-    // is just a copy of the container pointer; the sources on both DET and SRC are the same
-    // memory objects
-    if (strcmp(STACK_SRC, STACK_DET)) {
-	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
-	    psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-	    return psphotReadoutCleanup (config, view, STACK_SRC);
-	}
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PEAKS")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    // psphotDumpTest (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    // psphotDumpTest (config, view, STACK_RAW);
     psMemDump("sourcestats");
     logMemStats("sourcestats");
@@ -164,60 +157,48 @@
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_SRC)) {
+    if (!psphotRoughClass (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, subtract radial profiles for the convolved
-    // image first.  The profiles found on the convolved image will be replaced by those
-    // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
-    // replace the the profiles (and re-subtract them) for the detection image, so we want
-    // to keep those versions of the profiles on the sources.
-    if (strcmp(STACK_SRC, STACK_DET)) {
-      // find and subtract radial profile models for saturated stars (XXX change name eventually)
-      if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
+	return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
 	psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
-      }
-    }
-    // find and subtract radial profile models for saturated stars (XXX change name eventually)
-    if (!psphotDeblendSatstars (config, view, STACK_DET)) {
-	psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // use bright stellar objects to measure PSF
-    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_SRC);
+    psphotMergeSources (config, view, STACK_RAW);
 
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
-    psphotStackVisualFilerule(config, view, STACK_SRC);
+    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
+    psphotStackVisualFilerule(config, view, STACK_RAW);
 
     // measure the radial profiles to the sky
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
@@ -225,16 +206,16 @@
     // but iterates to an appropriately larger size
     logMemStats("before.kron.1");
-    psphotKronIterate(config, view, STACK_SRC, 1);
+    psphotKronIterate(config, view, STACK_RAW, 1);
     logMemStats("after.kron.1");
 	
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_SRC, true);
+    psphotSourceSize (config, view, STACK_RAW, true);
 
     // non-linear PSF and EXT fit to brighter sources
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
-    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // replace all sources (do NOT ignore subtraction state)
-    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
 
     logMemStats("pass1");
@@ -245,6 +226,6 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
+    // NOTE 2 : this function subtracts the models from the given filerule
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -252,22 +233,10 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
-	// if DET and SRC are different images, generate children sources for all sources in
-	// the SRC image.  This operation replaces the existing DETECTION container on DET
-	// which is currently a view to the one on SRC).  children sources go to
-	// det->allSources
-	if (strcmp(STACK_SRC, STACK_DET)) {
-	    psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
-
-	    //  subtract all sources from DET (this will subtract using the psf model for SRC, which
-	    //  will somewhat oversubtract the sources -- this is OK
-	    psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
-	}
-
 	// add noise for subtracted objects
-	psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+	psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
 	// find fainter sources
 	// NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
-	psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+	psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
 
 	// remove noise for subtracted objects (ie, return to normal noise level)
@@ -276,43 +245,30 @@
         bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
         if (!footprintsUseUnsubtracted) {
-    	    psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+    	    psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
         }
-
-	// if DET and SRC are different images, copy the detections from DET to SRC 
-	// (this operation just ensures the metadata container has a view on SRC as well
-	if (strcmp(STACK_SRC, STACK_DET)) {
-	    // replace all sources in DET
-	    psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
-
-	    // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
-	    if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
-		psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-		return psphotReadoutCleanup (config, view, STACK_SRC);
-	    }
-	}
 
 	// define new sources based on only the new peaks & measure moments
 	// NOTE: new sources are saved on detections->newSources
-	psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+	psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
 
 	// set source type
 	// NOTE: apply only to detections->newSources
-	if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+	if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
 	    psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-	    return psphotReadoutCleanup (config, view, STACK_SRC);
+	    return psphotReadoutCleanup (config, view, STACK_RAW);
 	}
 
 	// replace all sources so fit below applies to all at once
 	// NOTE: apply only to OLD sources (which have been subtracted)
-	psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+	psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
 
 	// merge the newly selected sources into the existing list
 	// NOTE: merge OLD and NEW
 	// XXX check on free of sources...
-	psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+	psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
 
 	// Construct an initial model for each object, set the radius to fitRadius, set circular
 	// fit mask.  NOTE: only applied to sources without guess models
-	psphotGuessModels (config, view, STACK_SRC);
+	psphotGuessModels (config, view, STACK_RAW);
     }
 
@@ -325,7 +281,7 @@
     if (splitLinearFit) {
         psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
-        // Fit the detected sources separately from matched that wea are about to create.
+        // Fit the detected sources separately from matched ones that we are about to create.
         // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
     }
 
@@ -335,5 +291,5 @@
     // this just match the detections for the chisq image, and not bother measuring the source
     // stats in that case...?
-    objects = psphotMatchSources (config, view, STACK_SRC);
+    objects = psphotMatchSources (config, view, STACK_RAW);
     psMemDump("matchsources");
 
@@ -344,26 +300,26 @@
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     psphotStackObjectsUnifyPosition (objects);
 
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
 
     // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
     // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
     logMemStats("before.kron.2");
-    psphotKronIterate(config, view, STACK_SRC, 2);
+    psphotKronIterate(config, view, STACK_RAW, 2);
     logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources
-    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
+    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
@@ -371,22 +327,27 @@
     // drop matched sources without any useful measurements and set kron radii for the ones
     // we decide to keep
-    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
 
     // measure kron fluxes for the matched sources only
-    psphotKronIterate(config, view, STACK_SRC, 3);
+    psphotKronIterate(config, view, STACK_RAW, 3);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
     // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
-    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // create source children for the OUT filerule (for radial aperture photometry and output) 
-    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    // NOTE: The new source children have image arrays pointing to the readout associated with 
+    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT, 
+    // but keep the pointers the same so we do not break these source image references
+    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
+    // sourcesSubtracted argument
+    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
     if (!objectsOut) {
 	psFree(objects);
 	psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
@@ -396,98 +357,82 @@
         // this forces photometry on the undetected sources from other images
 
-        // NOTE: we always do the radial apertures analysis on the convolved image since
-        // those are the ones that are psf matched and are the source of STACK_OUT's pixels
-        // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
-        // but in that case we don't get past the psf matching step because there is no
-        // target psf for the RAW inputs
-
-        // If useRaw copy the sources to the convolved readout
-        if (strcmp(STACK_SRC, STACK_CNV)) {
-            if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
-        // mark any inputs that we want to skip the matched apertures for
-        psphotStackSetInputsToSkip(config, view, STACK_CNV, true);
-        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
-        psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
-        psMemDump("extmeas");
-
+	// set up the FWHM vector
+	psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+	psphotDumpImages (config, view, STACK_RAW, "raw.t0");
+	psphotDumpImages (config, view, STACK_OUT, "out.t0");
 
         int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
-        for (int entry = 1; entry < nRadialEntries; entry++) {
+
+        for (int entry = 0; entry < nRadialEntries; entry++) {
             // NOTE: entry 0 is the unmatched image set
 
-            // re-measure the PSF for the smoothed image (using entries in 'allSources')
-            psphotChoosePSF (config, view, STACK_OUT, false);
-
-            // this is necessary to update the models based on the new PSF
-            psphotResetModels (config, view, STACK_OUT);
-
-            // this is necessary to get the right normalization for the new models
-            psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+	    char line[256];
 
             // measure circular, radial apertures (objects sorted by S/N)
             psphotRadialApertures (config, view, STACK_OUT, entry); 
+	    snprintf (line, 256, "%s.%d", "out.t1", entry);
+	    psphotDumpImages (config, view, STACK_OUT, line);
 
             // replace the flux in the image so it is returned to its original state
             psphotReplaceAllSources (config, view, STACK_OUT, false);
-
-            // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
-            psMemDump("matched");
+	    snprintf (line, 256, "%s.%d", "out.t2", entry);
+	    psphotDumpImages (config, view, STACK_OUT, line);
+
+	    if (entry < nRadialEntries - 1) {
+		// smooth to the next FWHM
+		// this function does nothing if the targetFWHM is smaller than the currentFWHM
+		psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
+		snprintf (line, 256, "%s.%d", "out.t3", entry);
+		psphotDumpImages (config, view, STACK_OUT, line);
+		psMemDump("matched");
+
+		// re-measure the PSF for the smoothed image (using entries in 'allSources')
+		psphotChoosePSF (config, view, STACK_OUT, false);
+
+		// this is necessary to update the models based on the new PSF
+		psphotResetModels (config, view, STACK_OUT);
+
+		// this is necessary to get the right normalization for the new models
+		// and to subtract the sources
+		psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+		snprintf (line, 256, "%s.%d", "out.t4", entry);
+		psphotDumpImages (config, view, STACK_OUT, line);
+	    }
         }
     }
-    psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
-    psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_SRC)) {
+    if (!psphotApResid (config, view, STACK_RAW)) {
 	psFree (objects);
 	psFree (objectsOut);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-	return psphotReadoutCleanup (config, view, STACK_SRC);
+	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
-
-    if (!useRaw) {
-        // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
-        // the convolved images we need to generate PSFs for the DET images
-        if (!psphotChoosePSF (config, view, STACK_DET, false)) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
-            return psphotReadoutCleanup (config, view, STACK_DET);
-        }
-    }
-    if (!psphotEfficiency(config, view, STACK_DET)) {
+    psphotMagnitudes(config, view, STACK_RAW);
+
+    if (!psphotEfficiency(config, view, STACK_RAW)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
     }
-    psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
+    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
 
     logMemStats("final");
 #if (1)
-    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_RAW);
     psphotSourceMemory(config, view, STACK_OUT);
 #endif
 
-    // replace failed sources?
-    // psphotReplaceUnfitSources (sources);
-
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_DET);
+    psphotSkyReplace (config, view, STACK_RAW);
 
     // drop the references to the image pixels held by each source
+    psphotSourceFreePixels (config, view, STACK_RAW);
     psphotSourceFreePixels (config, view, STACK_OUT);
-    psphotSourceFreePixels (config, view, STACK_SRC);
 
 #ifdef MAKE_CHISQ_IMAGE
     // remove chisq image from config->file:PSPHOT.INPUT
-    psphotStackRemoveChisqFromInputs(config, STACK_DET);
-    if (strcmp(STACK_SRC, STACK_DET)) {
-	psphotStackRemoveChisqFromInputs(config, STACK_SRC);
-    }
+    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
 #endif
 
@@ -496,5 +441,5 @@
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_SRC);
+    return psphotReadoutCleanup (config, view, STACK_RAW);
 }
 
@@ -552,11 +497,8 @@
 }
 
-
-
 /* here is the process:
 
- * we have three(*) images:
- * RAW : unconvolved image stack 
- * CNV : input convolved image
+ * we have two image sets:
+ * RAW : unconvolved image stacks
 
  * OUT : psf-matched output image (there may be more than one of
@@ -649,2 +591,124 @@
 
    */
+
+
+// generate a vector fwhmValues where the first has the fwhm of the raw image and the
+// successive entries have the target values
+
+bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
+
+    bool status;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+	if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
+	    return false;
+	}
+    }
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
+
+// copy the pixels from RAW to OUT (
+
+bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
+
+    bool status;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // find the currently selected readout
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
+    psAssert (fileRaw, "missing file?");
+
+    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
+    psAssert (readoutRaw, "missing readout?");
+
+    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
+    if (readoutRaw->variance) {
+	readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
+    }
+    if (readoutRaw->mask) {
+	readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
+    }
+
+    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
+    // psAssert (chipRaw, "missing chip");
+
+    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
+    if (!psf) {
+	// we should have a PSF by this point in psphot
+	psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
+	return false;
+    }
+
+    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
+
+    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    psVectorAppend(fwhmValues, fwhmRaw);
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+    if (!status) {
+	psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+	psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+	for (int i = 0; i < targetSeeing->n; i++) {
+	    psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
+	}	    
+    } else {
+        psVectorAppend(fwhmValues, targetSeeing);
+    }
+
+    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree (fwhmValues);
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
+
+    psImage *image = readout->image;
+
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    float sumFWHM = 0.0;		  // FWHM for image
+    int numFWHM = 0;			  // Number of FWHM measurements
+    for (float x = 0; x < Nx; x += 0.25*Nx) {
+	for (float y = 0; y < Ny; y += 0.25*Ny) {
+	    float fwhm = pmPSFtoFWHM(psf, x, y);
+	    if (isfinite(fwhm)) {
+		sumFWHM += fwhm;
+		numFWHM++;
+	    }
+	}
+    }
+    if (numFWHM == 0) {
+	psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
+	return NAN;
+    } 
+
+    float fwhm = sumFWHM / (float) numFWHM;
+    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
+
+    return fwhm;
+}