Index: trunk/psphot/src/psphotStackReadout.c
===================================================================
--- trunk/psphot/src/psphotStackReadout.c	(revision 32633)
+++ trunk/psphot/src/psphotStackReadout.c	(revision 32695)
@@ -82,8 +82,4 @@
     }
 
-    // XXX I think this is not defined correctly for an array of images.
-    // XXX I probably need to subtract the model (same model?) for both RAW and OUT.
-    // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level.
-
     // generate a background model (median, smoothed image)
     if (!psphotModelBackground (config, view, STACK_DET)) {
@@ -91,4 +87,7 @@
     }
     if (!psphotSubtractBackground (config, view, STACK_DET)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_SRC)) {
 	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
@@ -114,5 +113,7 @@
     }
 
-    // if DET and SRC are different images, copy the detections from DET to 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)) {
@@ -130,4 +131,5 @@
 	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+    // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
 
@@ -166,6 +168,4 @@
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    // XXX why do this as a stack operation?
-    // psphotFitSourcesLinearStack (config, objects, false);
     psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
@@ -186,6 +186,6 @@
     psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
-    // replace all sources
-    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    // replace all sources (do NOT ignore subtraction state)
+    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
@@ -194,4 +194,5 @@
     // 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); // pass 2 (detections->allSources)
 
@@ -200,4 +201,16 @@
     // 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); // ignore subtraction state for sources
+	}
+
 	// add noise for subtracted objects
 	psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
@@ -212,7 +225,11 @@
 
 	// 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)) {
-	    // XXX how does this handle 1st vs 2nd pass sources?
-	    if (!psphotCopySources (config, view, 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);
@@ -237,5 +254,5 @@
 	// 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); // pass 2
+	psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
 
 	// merge the newly selected sources into the existing list
@@ -312,5 +329,5 @@
 
 	// replace the flux in the image so it is returned to its original state
-	psphotReplaceAllSources (config, view, STACK_OUT);
+	psphotReplaceAllSources (config, view, STACK_OUT, false);
 
 	// smooth to the next FWHM, or set 'smoothAgain' to false if no more