Index: trunk/psphot/src/psphotFullForce.SourceStats.c
===================================================================
--- trunk/psphot/src/psphotFullForce.SourceStats.c	(revision 36630)
+++ trunk/psphot/src/psphotFullForce.SourceStats.c	(revision 36633)
@@ -85,5 +85,6 @@
     psAssert (markVal, "missing MARK.PSPHOT");
 
-    psArray *sources = detections->allSources;
+    // psArray *sources = detections->allSources;
+    psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
 
     // generate the array of sources, define the associated pixel
Index: trunk/psphot/src/psphotFullForceReadout.c
===================================================================
--- trunk/psphot/src/psphotFullForceReadout.c	(revision 36630)
+++ trunk/psphot/src/psphotFullForceReadout.c	(revision 36633)
@@ -52,7 +52,10 @@
     psphotLoadExtSources (config, view, filerule);
 
+#ifdef notmoved
+    // XXX: moved down below
     // merge the newly selected sources into the existing list (detections->allSources)
     // NOTE: merge OLD and NEW
     psphotMergeSources (config, view, filerule);
+#endif
 
     // construct sources and measure moments and other basic stats (saved on detections->allSources)
@@ -63,4 +66,16 @@
         return psphotReadoutCleanup (config, view, filerule);
     }
+    
+    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
+    // clump defined for it is used not measured (detections->newSources)
+    if (!psphotRoughClass (config, view, filerule)) { // pass 1
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
+    // merge the newly selected sources into the existing list (detections->allSources)
+    // NOTE: merge OLD and NEW
+    psphotMergeSources (config, view, filerule);
+
 
     // generate a psf model for any readouts which need one
Index: trunk/psphot/src/psphotFullForceSummary.c
===================================================================
--- trunk/psphot/src/psphotFullForceSummary.c	(revision 36630)
+++ trunk/psphot/src/psphotFullForceSummary.c	(revision 36633)
@@ -78,4 +78,18 @@
         usage();
     }
+    if ((N = psArgumentGet(argc, argv, "-cff"))) {
+        if (argc <= N+1) {
+          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+          usage();
+        }
+        psArgumentRemove(N, &argc, argv);
+
+        // We read CFF file as table not as sources. All we want are the nominal parameters input to psphotFullForce
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "CFF_FILE", 0, "CFF file with nominal parameters", argv[N]);
+        psArgumentRemove(N, &argc, argv);
+    } else {
+        psErrorStackPrint(stderr, "-cff must be supplied.");
+        usage();
+    }
  
     if (argc < 2) {
@@ -148,4 +162,45 @@
     }
     output->save = true;
+
+    // We read cff file as a table not as sources
+    psString cffName = psMetadataLookupStr(NULL, config->arguments, "CFF_FILE");
+    if (!cffName) {
+        psError(PSPHOT_ERR_CONFIG, true, "CFF_FILE is missing from arguments");
+        return false;
+    }
+    psString resolvedName = pmConfigConvertFilename(cffName, config, false, false);
+    if (!resolvedName) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to resolve CFF_FILE %s", cffName);
+        return false;
+    }
+    psFits *fits = psFitsOpen(resolvedName, "r");
+    if (!fits) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to open fits CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psArray *inTable = psFitsReadTable(fits);
+    if (!inTable) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to read cff fits table from CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psFitsClose(fits);
+
+    // Convert to an array indexed by ID
+    psArray *sortedTable = psArrayAlloc(4*inTable->n);
+    for (int i=0; i<inTable->n; i++) {
+        psMetadata *row = inTable->data[i];
+        psS32 ID = psMetadataLookupS32(&status, row, "ID");
+        if (ID >= sortedTable->n) {
+            sortedTable = psArrayRealloc(sortedTable, 2*ID);
+        }
+        if (sortedTable->data[ID]) {
+            psError(PSPHOT_ERR_CONFIG, true, "Duplicate row with ID %d", ID);
+            return false;
+        }
+        sortedTable->data[ID] = psMemIncrRefCounter(row);
+    }
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
+    psFree(inTable);
+    psFree(sortedTable);
 
     return true;
Index: trunk/psphot/src/psphotFullForceSummaryReadout.c
===================================================================
--- trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36630)
+++ trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36633)
@@ -14,4 +14,5 @@
     psArray *zeroPt;
     psArray *exptime;
+    psArray *cffTable;
 } galaxyShapeOptions;
 
@@ -69,4 +70,9 @@
             if (!setOptions(&options, readout, recipe, true)) {
                 psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options.");
+                return false;
+            }
+            options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
+            if (!options.cffTable) {
+                psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
                 return false;
             }
@@ -309,6 +315,30 @@
                 // copy the best fit params to the model
                 psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
-                axes.major *= options->rMajor->data.F32[min_j];
-                axes.minor *= options->rMinor->data.F32[min_j];
+                if (outSrc->seq >= options->cffTable->n) {
+                    psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
+                        outSrc->seq, options->cffTable->n);
+                    return NULL;
+                }
+                psMetadata *row = options->cffTable->data[outSrc->seq];
+                if (!row) {
+                    psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                bool mdok;
+                psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR");
+                if (!mdok) {
+                    psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR");
+                if (!mdok) {
+                    psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
+                    psFree(outSrc);
+                    return NULL;
+                }
+                axes.major = rMajor * options->rMajor->data.F32[min_j];
+                axes.minor = rMinor * options->rMinor->data.F32[min_j];
                 pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
                 outSrc->modelEXT->chisq = minChisq; 
Index: trunk/psphot/src/psphotGalaxyShape.c
===================================================================
--- trunk/psphot/src/psphotGalaxyShape.c	(revision 36630)
+++ trunk/psphot/src/psphotGalaxyShape.c	(revision 36633)
@@ -205,4 +205,5 @@
 	if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
 	if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        if (!source->modelPSF) continue;
 
 	// psphotSetRadiusMomentsExact sets the radius based on Mrf
@@ -274,13 +275,19 @@
 
             int i = source->galaxyFits->chisq->n - 1;
-            float thisChisq = source->galaxyFits->chisq->data.F32[i];
-            if (isfinite(thisChisq) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
-                chisqBest = thisChisq;
-                fRmajorBest = fRmajor;
-                fRminorBest = fRminor;
+            float flux = source->galaxyFits->Flux->data.F32[i];
+            if (isfinite(flux)) {
+                float thisChisq = source->galaxyFits->chisq->data.F32[i];
+                if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                    chisqBest = thisChisq;
+                    fRmajorBest = fRmajor;
+                    fRminorBest = fRminor;
+                }
             }
+            // reset I0 to avoid potential problems on the next iteration
+            PAR[PM_PAR_I0] = 1.0;
 	}
     }
 
+#define SAVE_BEST_MODEL
 #ifdef SAVE_BEST_MODEL
     // Save model with smallest chisq
@@ -301,15 +308,16 @@
         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
 
-        // Replace modelEXT with this model
-        // XXX: only do this if the model is good
-        psFree (source->modelEXT);
-
-        source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
-        source->mode |= PM_SOURCE_MODE_EXTMODEL;
-        source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-        // cache the model flux
-        pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        // Replace modelEXT with this model, if the model is good
+        if (isfinite(PAR[PM_PAR_I0])) {
+            psFree (source->modelEXT);
+
+            source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            source->mode |= PM_SOURCE_MODE_EXTMODEL;
+            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+
+            // cache the model flux
+            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        }
     }