Index: trunk/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- trunk/psphot/src/psphotExtendedSourceFits.c	(revision 30624)
+++ trunk/psphot/src/psphotExtendedSourceFits.c	(revision 31154)
@@ -43,4 +43,5 @@
     int NplainPass = 0;
     int Nfaint = 0;
+    int Nfail = 0;
 
     psTimerStart ("psphot.extended");
@@ -145,5 +146,5 @@
 
     // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
@@ -176,32 +177,38 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
-
-            if (false && !psThreadJobAddPending(job)) {
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+// set this to 0 to run without threading
+# if (1)	    
+            if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
 		psFree(AnalysisRegion);
                 return false;
-            } else {
-		if (!psphotExtendedSourceFits_Threaded(job)) {
-		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-		    psFree(AnalysisRegion);
-		    return false;
-		}
-                psScalar *scalar = NULL;
-                scalar = job->args->data[7];
-                Next += scalar->data.S32;
-                scalar = job->args->data[8];
-                Nconvolve += scalar->data.S32;
-                scalar = job->args->data[9];
-                NconvolvePass += scalar->data.S32;
-                scalar = job->args->data[10];
-                Nplain += scalar->data.S32;
-                scalar = job->args->data[11];
-                NplainPass += scalar->data.S32;
-                scalar = job->args->data[12];
-                Nfaint += scalar->data.S32;
-		psFree(job);
+            } 
+# else
+	    if (!psphotExtendedSourceFits_Threaded(job)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		psFree(AnalysisRegion);
+		return false;
 	    }
-        }
+	    psScalar *scalar = NULL;
+	    scalar = job->args->data[7];
+	    Next += scalar->data.S32;
+	    scalar = job->args->data[8];
+	    Nconvolve += scalar->data.S32;
+	    scalar = job->args->data[9];
+	    NconvolvePass += scalar->data.S32;
+	    scalar = job->args->data[10];
+	    Nplain += scalar->data.S32;
+	    scalar = job->args->data[11];
+	    NplainPass += scalar->data.S32;
+	    scalar = job->args->data[12];
+	    Nfaint += scalar->data.S32;
+	    scalar = job->args->data[13];
+	    Nfail += scalar->data.S32;
+	    psFree(job);
+# endif
+	}
 
         // wait for the threads to finish and manage results
@@ -231,4 +238,6 @@
                 scalar = job->args->data[12];
                 Nfaint += scalar->data.S32;
+                scalar = job->args->data[13];
+                Nfail += scalar->data.S32;
             }
             psFree(job);
@@ -238,8 +247,8 @@
     psFree(AnalysisRegion);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
     return true;
 }
@@ -250,8 +259,9 @@
     bool status;
     int Next = 0;
+    int Nfaint = 0;
+    int Nfail = 0;
     int Nconvolve = 0;
     int NconvolvePass = 0;
     int Nplain = 0;
-    int Nfaint = 0;
     int NplainPass = 0;
     bool savePics = false;
@@ -268,7 +278,12 @@
     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
 
+    // pthread_t tid = pthread_self();     // Thread identifier
+
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode          = PM_SOURCE_FIT_EXT;
+    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+
     // XXX for now, use the defaults for the rest:
     // fitOptions->nIter         = fitIter;
@@ -296,4 +311,6 @@
         // XXX use the parameters guessed from moments
         // if (source->modelEXT == NULL) continue;
+
+	// fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
 
         // replace object in image
@@ -366,5 +383,5 @@
           // limit selection to some SN limit
           assert (source->peak); // how can a source not have a peak?
-          if (source->peak->SN < SNlim) {
+          if (sqrt(source->peak->detValue) < SNlim) {
 	      Nfaint ++;
 	      continue;
@@ -388,4 +405,5 @@
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
+		  Nfail ++;
                   continue;
               }
@@ -402,4 +420,5 @@
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
+		  Nfail ++;
                   continue;
               }
@@ -519,4 +538,7 @@
     scalar->data.S32 = Nfaint;
 
+    scalar = job->args->data[13];
+    scalar->data.S32 = Nfail;
+
     return true;
 }