Index: branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c
===================================================================
--- branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c	(revision 26691)
+++ branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c	(revision 26748)
@@ -42,4 +42,12 @@
     psAssert (!detections->newSources, "new sources already defined?");
 
+    // XXX TEST:
+    if (detections->allSources) {
+	psphotMaskCosmicRayFootprintCheck(detections->allSources);
+    }
+    if (detections->newSources) {
+	psphotMaskCosmicRayFootprintCheck(detections->newSources);
+    }
+
     // determine the number of allowed threads
     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
@@ -57,4 +65,15 @@
     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
     psAssert (status, "missing BREAK_POINT?");
+
+    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); psAssert (status, "missing SKY_INNER_RADIUS");
+    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    psAssert (maskVal, "missing MASK.PSPHOT");
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    psAssert (markVal, "missing MARK.PSPHOT");
 
     psArray *peaks = detections->peaks;
@@ -112,4 +131,14 @@
     }
 
+    // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis
+    float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
+    if (!status) {
+	SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    }
+    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
+    if (!status) {
+	RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    }
+
     // threaded measurement of the source magnitudes
     int Nfail = 0;
@@ -133,5 +162,13 @@
 
             psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, recipe);
+
+            PS_ARRAY_ADD_SCALAR(job->args, INNER,   PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, MIN_SN,  PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
+
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
@@ -152,4 +189,137 @@
 	    psFree(sources);
             return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[7];
+                Nmoments += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nfail += scalar->data.S32;
+                scalar = job->args->data[9];
+                Nfaint += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+
+    psFree (cellGroups);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
+
+    psphotVisualShowMoments (sources);
+
+    // save the new sources on the detection structure:
+    detections->newSources = sources;
+
+
+    if (detections->allSources) {
+	psphotMaskCosmicRayFootprintCheck(detections->allSources);
+    }
+    if (detections->newSources) {
+	psphotMaskCosmicRayFootprintCheck(detections->newSources);
+    }
+
+    return true;
+}
+
+// this function is currently only called by psphotCheckExtSources 
+bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
+
+    bool status = false;
+
+    psTimerStart ("psphot.stats");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // determine properties (sky, moments) of initial sources
+    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    if (!status) return false;
+
+    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
+
+    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
+    if (!status) return false;
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+        // allocate space for moments
+	if (!source->moments) {
+	    source->moments = pmMomentsAlloc();
+	}
+
+        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+        source->peak->assigned = true;
+    }
+
+    if (!strcasecmp (breakPt, "PEAKS")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
+        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
+        psphotVisualShowMoments (sources);
+        return sources;
+    }
+
+    // XXX how else could we get the window size in?
+    if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
+	psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
+	return NULL;
+    }
+
+    // threaded measurement of the source magnitudes
+    int Nfail = 0;
+    int Nmoments = 0;
+    int Nfaint = 0;
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
+
+	    // XXX: this must match the above
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, recipe);
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
+                psFree (job);
+                return NULL;
+            }
+            psFree(job);
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
+            return NULL;
         }
 
@@ -178,127 +348,4 @@
     psphotVisualShowMoments (sources);
 
-    // save the new sources on the detection structure:
-    detections->newSources = sources;
-    return true;
-}
-
-// XXX fix the return status of this function...
-bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
-
-    bool status = false;
-
-    psTimerStart ("psphot.stats");
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    assert (recipe);
-
-    // determine the number of allowed threads
-    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
-    if (!status) {
-        nThreads = 0;
-    }
-
-    // determine properties (sky, moments) of initial sources
-    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    if (!status) return false;
-
-    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
-
-    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
-    if (!status) return false;
-
-    for (int i = 0; i < sources->n; i++) {
-
-        pmSource *source = sources->data[i];
-        if (!source->peak) continue; // XXX how can we have a peak-less source?
-
-        // allocate space for moments
-	if (!source->moments) {
-	    source->moments = pmMomentsAlloc();
-	}
-
-        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
-        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
-        source->peak->assigned = true;
-    }
-
-    if (!strcasecmp (breakPt, "PEAKS")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
-        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
-        psphotVisualShowMoments (sources);
-        return sources;
-    }
-
-    // XXX how else could we get the window size in?
-    if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
-	psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
-	return NULL;
-    }
-
-    // threaded measurement of the source magnitudes
-    int Nfail = 0;
-    int Nmoments = 0;
-    int Nfaint = 0;
-
-    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
-    int Cx = 1, Cy = 1;
-    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
-
-    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
-
-    for (int i = 0; i < cellGroups->n; i++) {
-
-        psArray *cells = cellGroups->data[i];
-
-        for (int j = 0; j < cells->n; j++) {
-
-            // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
-
-            psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, recipe);
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
-
-            if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
-                psFree (job);
-                return NULL;
-            }
-            psFree(job);
-        }
-
-        // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
-            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
-            return NULL;
-        }
-
-        // we have only supplied one type of job, so we can assume the types here
-        psThreadJob *job = NULL;
-        while ((job = psThreadJobGetDone()) != NULL) {
-            if (job->args->n < 1) {
-                fprintf (stderr, "error with job\n");
-            } else {
-                psScalar *scalar = NULL;
-                scalar = job->args->data[2];
-                Nmoments += scalar->data.S32;
-                scalar = job->args->data[3];
-                Nfail += scalar->data.S32;
-                scalar = job->args->data[4];
-                Nfaint += scalar->data.S32;
-            }
-            psFree(job);
-        }
-    }
-
-    psFree (cellGroups);
-
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
-
-    psphotVisualShowMoments (sources);
-
     return (sources);
 }
@@ -311,22 +358,12 @@
 
     psArray *sources                = job->args->data[0];
-    psMetadata *recipe              = job->args->data[1];
-
-    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return false;
-    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-    if (!status) return false;
-    float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
-    if (!status) return false;
-    float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
-    if (!status) return false;
-
-    // bit-masks to test for good/bad pixels
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
-
-    // bit-mask to mark pixels not used in analysis
-    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
+
+    float INNER                     = PS_SCALAR_VALUE(job->args->data[1],F32);
+    float MIN_SN                    = PS_SCALAR_VALUE(job->args->data[2],F32);
+    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
+    float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
+
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -392,11 +429,11 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[2];
+    scalar = job->args->data[7];
     scalar->data.S32 = Nmoments;
 
-    scalar = job->args->data[3];
+    scalar = job->args->data[8];
     scalar->data.S32 = Nfail;
 
-    scalar = job->args->data[4];
+    scalar = job->args->data[9];
     scalar->data.S32 = Nfaint;
 
@@ -410,6 +447,7 @@
     bool status;
 
-    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-    if (!status) return false;
+    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
+    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
+    psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
 
     // XXX move this to a config file?
@@ -441,10 +479,10 @@
 
         // choose a grid scale that is a fixed fraction of the psf sigma^2
-        psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
-        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
-        psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+        float PSF_CLUMP_GRID_SCALE = 0.1*PS_SQR(sigma[i]);
+        float MOMENTS_SX_MAX = 2.0*PS_SQR(sigma[i]);
+        float MOMENTS_SY_MAX = 2.0*PS_SQR(sigma[i]);
 
         // determine the PSF parameters from the source moment values
-        pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
+        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
@@ -485,9 +523,9 @@
 
     // choose a grid scale that is a fixed fraction of the psf sigma^2
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
 
     psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);