Index: /trunk/ippconfig/gpc1/camera.config
===================================================================
--- /trunk/ippconfig/gpc1/camera.config	(revision 25382)
+++ /trunk/ippconfig/gpc1/camera.config	(revision 25383)
@@ -132,4 +132,5 @@
   CMF.XSRC STR {CHIP.NAME}.xsrc # use .PSF and .EXT?
   CMF.XFIT STR {CHIP.NAME}.xfit # use .PSF and .EXT?
+  CMF.DETEFF STR {CHIP.NAME}.deteff
 
   PSF.HEAD  STR	{CHIP.NAME}.hdr
Index: /trunk/ippconfig/recipes/psphot.config
===================================================================
--- /trunk/ippconfig/recipes/psphot.config	(revision 25382)
+++ /trunk/ippconfig/recipes/psphot.config	(revision 25383)
@@ -265,4 +265,7 @@
 PSPHOT.CR.NSIGMA.SOFTEN             F32   0.0025          # Softening parameter for weights
 
+# Detection efficiency
+EFF.NUM                             S32   500			# Number of fake sources per bin
+@EFF.MAG			    F32	  -2.0 -1.0 -0.5 -0.25 -0.1 -0.05 0.0 0.05 0.1 0.25 0.5 1.0 2.0	# Magnitude of fake sources relative to limit
 
 # Recipe overrides for CHIP
Index: /trunk/psLib/src/fits/psFits.c
===================================================================
--- /trunk/psLib/src/fits/psFits.c	(revision 25382)
+++ /trunk/psLib/src/fits/psFits.c	(revision 25383)
@@ -344,6 +344,8 @@
 // Therefore, we implement our own version of moving to an extension specified by name.  The pure cfitsio
 // version is used if "conventions.compression" handling is turned off in the psFits structure.
-bool psFitsMoveExtName(const psFits* fits,
-                       const char* extname)
+static bool fitsMoveExtName(const psFits* fits, // FITS file
+                            const char* extname, // Extension name
+                            bool errors // Generate errors?
+    )
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -356,5 +358,7 @@
         // User wants to use cfitsio.  Good luck to them!
         if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
-            psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            }
             return false;
         }
@@ -378,5 +382,7 @@
         if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
             // We've run off the end
-            psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            }
             return false;
         }
@@ -406,4 +412,15 @@
     }
     psAbort("Should never reach here.");
+}
+
+
+bool psFitsMoveExtName(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, true);
+}
+
+bool psFitsMoveExtNameClean(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, false);
 }
 
Index: /trunk/psLib/src/fits/psFits.h
===================================================================
--- /trunk/psLib/src/fits/psFits.h	(revision 25382)
+++ /trunk/psLib/src/fits/psFits.h	(revision 25383)
@@ -245,4 +245,14 @@
 );
 
+/** Moves the FITS HDU to the specified extension name without generating errors.
+ *
+ *  @return bool        TRUE if the extension name was found and move was
+ *                      successful, otherwise FALSE
+ */
+bool psFitsMoveExtNameClean(
+    const psFits* fits,                ///< the psFits object to move
+    const char* extname                ///< the extension name
+);
+
 /** Moves the FITS HDU to the specified extension number
  *
Index: /trunk/psLib/src/imageops/psImageConvolve.c
===================================================================
--- /trunk/psLib/src/imageops/psImageConvolve.c	(revision 25382)
+++ /trunk/psLib/src/imageops/psImageConvolve.c	(revision 25383)
@@ -678,4 +678,137 @@
 }
 
+static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
+                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
+                                  const psVector *gaussNorm, float minGauss, int size, int start, int stop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+    for (int i = start; i < stop; i++) {
+        int xPix = x->data.S32[i], yPix = y->data.S32[i]; // Pixel coordinates for smoothing
+
+        int yMin = PS_MAX(yPix - size, 0);
+        int yMax = PS_MIN(yPix + size, yLast);
+        int xMin = PS_MAX(xPix - size, 0);
+        int xMax = PS_MIN(xPix + size, xLast);
+
+        const float *yGauss = &gauss[yMin - yPix];
+
+        double ySumIG = 0.0, ySumG = 0.0;
+        for (int v = yMin; v <= yMax; v++, yGauss++) {
+            const float *xGauss = &gauss[xMin - xPix];
+            double xSumIG = 0.0, xSumG = 0.0;
+            const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[v][xMin];
+            const psF32 *imageData = &image->data.F32[v][xMin];
+            for (int u = xMin; u <= xMax; u++, xGauss++, imageData++, maskData++) {
+                if (*maskData & maskVal) {
+                    continue;
+                }
+                xSumIG += *imageData * *xGauss;
+                xSumG += *xGauss;
+            }
+            if (xSumG > minGauss) {
+                ySumIG += xSumIG * *yGauss;
+                ySumG += xSumG * *yGauss;
+            }
+        }
+
+        out->data.F32[i] = ySumG > minGauss ? ySumIG / ySumG : NAN;
+    }
+
+    return true;
+}
+
+static bool psImageSmoothMaskPixelsThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "programming error: no job arguments");
+    psAssert(job->args->n == 11, "programming error: wrong number of job arguments");
+
+    psVector *out = job->args->data[0]; // Output vector
+    const psImage *image  = job->args->data[1]; // Input image
+    const psImage *mask   = job->args->data[2]; // Input mask
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3], PS_TYPE_IMAGE_MASK_DATA);
+    const psVector *x = job->args->data[4];
+    const psVector *y = job->args->data[5];
+    const psVector *gaussNorm = job->args->data[6];
+    float minGauss = PS_SCALAR_VALUE(job->args->data[7], F32);
+    int size = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int start = PS_SCALAR_VALUE(job->args->data[9], S32);
+    int stop = PS_SCALAR_VALUE(job->args->data[10], S32);
+    return imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm,
+                                 minGauss, size, start, stop);
+}
+
+
+psVector *psImageSmoothMaskPixels(const psImage *image, const psImage *mask, psImageMaskType maskVal,
+                                  const psVector *x, const psVector *y,
+                                  float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    int num = x->n;                     // Number of pixels to smooth
+    psVector *out = psVectorAlloc(num, PS_TYPE_F32); // Output results
+
+    // Generate normalized gaussian
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    // Columns
+    if (threaded) {
+        int numThreads = psThreadPoolSize(); // Number of threads
+        int delta = (numThreads) ? num / numThreads + 1 : num; // Block of cols to do at once
+        for (int start = 0; start < num; start += delta) {
+            int stop = PS_MIN(start + delta, num);  // End of block
+
+            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
+            psArrayAdd(job->args, 1, out);
+            psArrayAdd(job->args, 1, (psImage*)image);
+            psArrayAdd(job->args, 1, (psImage*)mask);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            psArrayAdd(job->args, 1, (psVector*)x);
+            psArrayAdd(job->args, 1, (psVector*)y);
+            psArrayAdd(job->args, 1, gaussNorm);
+            PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(gaussNorm);
+                psFree(out);
+                return NULL;
+            }
+            psFree(job);
+        }
+    } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y,
+                                      gaussNorm, minGauss, size, 0, num)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        psFree(gaussNorm);
+        psFree(out);
+        return NULL;
+    }
+
+    psFree(gaussNorm);
+    return out;
+}
+
+
 // Smooth an image with masked pixels
 // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
@@ -1325,7 +1458,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
+        int deltaRows = (numThreads) ? numRows / numThreads + 1 : numRows; // Block of rows to do at once
         for (int start = 0; start < numRows; start += deltaRows) {
-            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
+            int stop = PS_MIN(start + deltaRows, numRows);  // end of row block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1363,7 +1496,7 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
+        int deltaCols = (numThreads) ? numCols / numThreads + 1 : numCols; // Block of cols to do at once
         for (int start = 0; start < numCols; start += deltaCols) {
-            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
+            int stop = PS_MIN(start + deltaCols, numCols);  // end of col block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1486,4 +1619,10 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
+            task->function = &psImageSmoothMaskPixelsThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");
Index: /trunk/psLib/src/imageops/psImageConvolve.h
===================================================================
--- /trunk/psLib/src/imageops/psImageConvolve.h	(revision 25382)
+++ /trunk/psLib/src/imageops/psImageConvolve.h	(revision 25383)
@@ -222,4 +222,18 @@
     );
 
+/// Smooth particular pixels on an image, allowing for masked pixels
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
+psVector *psImageSmoothMaskPixels(
+    const psImage *image,               ///< Input image (F32)
+    const psImage *mask,                ///< Mask image
+    psImageMaskType maskVal,            ///< Value to mask
+    const psVector *x,                  ///< x coordinates
+    const psVector *y,                  ///< y coordinates
+    float sigma,                        ///< Width of the smoothing kernel (pixels)
+    float numSigma,                     ///< Size of the smoothing box (sigma)
+    float minGauss                      ///< Minimum fraction of Gaussian to accept
+    );
 
 psImage *psImageSmoothMask_Threaded(psImage *output,
Index: /trunk/psLib/src/types/psMetadata.c
===================================================================
--- /trunk/psLib/src/types/psMetadata.c	(revision 25382)
+++ /trunk/psLib/src/types/psMetadata.c	(revision 25383)
@@ -996,5 +996,4 @@
     if (metadataItem->type == PS_DATA_METADATA_MULTI) {
         // if multiple keys found, use the first.
-        //        metadataItem = (psMetadataItem*)((metadataItem->data.list)->head);
         metadataItem = (psMetadataItem*)(metadataItem->data.list->head->data);
         if (status) {
@@ -1082,5 +1081,4 @@
     } \
     \
-    /* psFree(metadataItem); currently, the lookup doesn't increment the ref count */ \
     return value; \
 }
@@ -1101,4 +1099,44 @@
 psMetadataLookupNumTYPE(VectorMaskType,VectorMask)
 psMetadataLookupNumTYPE(ImageMaskType,ImageMask)
+
+#define psMetadataLookupPtrTYPE(TYPENAME,NAME,TYPE,VAL) \
+TYPENAME psMetadataLookup##NAME(bool *status, const psMetadata *md, const char *key) \
+{ \
+    PS_ASSERT_METADATA_NON_NULL(md, NULL); \
+    PS_ASSERT_STRING_NON_EMPTY(key, NULL); \
+    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); /* The item of interest */ \
+    if (!item) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (item->type == PS_DATA_METADATA_MULTI) { \
+        /* if multiple keys found, use the first. */ \
+        item = item->data.list->head->data; \
+    } \
+    if (item->type != TYPE) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (status) { \
+        *status = true; \
+    } \
+    return item->data.VAL; \
+}
+
+psMetadataLookupPtrTYPE(psMetadata*, Metadata, PS_DATA_METADATA, md)
+psMetadataLookupPtrTYPE(psString, Str, PS_DATA_STRING, str)
+psMetadataLookupPtrTYPE(psTime*, Time, PS_DATA_TIME, V)
+psMetadataLookupPtrTYPE(psVector*, Vector, PS_DATA_VECTOR, V)
+
 
 psMetadataItem* psMetadataGet(const psMetadata *md,
@@ -1257,105 +1295,7 @@
 }
 
-psMetadata *psMetadataLookupMetadata(bool *status,
-                                     const psMetadata *md,
-                                     const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    psMetadata *value = NULL;  // The value to return
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_METADATA) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.md; // The requested metadata
-    }
-    return value;
-}
-
-psTime *psMetadataLookupTime(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key);
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-        return NULL;
-    }
-
-    if (item->type != PS_DATA_TIME) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_TIME, as expected.\n", key);
-        }
-        return NULL;
-    }
-
-    // We have the requested metadata
-    if (status) {
-        *status = true;
-    }
-
-    return item->data.V;
-}
-
-
-psString psMetadataLookupStr(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    //    char *value = NULL;   // The value to return
-    psString value = NULL;
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_STRING) {
-        // The value at the key isn't of the desired type
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_STRING, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.V;
-    }
-    return value;
-}
+
+
+
 
 psList *psMetadataKeys(psMetadata *md)
Index: /trunk/psLib/src/types/psMetadata.h
===================================================================
--- /trunk/psLib/src/types/psMetadata.h	(revision 25382)
+++ /trunk/psLib/src/types/psMetadata.h	(revision 25383)
@@ -717,4 +717,5 @@
 PS_METADATA_LOOKUP_TYPE_DECL(Ptr, psPtr);
 PS_METADATA_LOOKUP_TYPE_DECL(Str, psString);
+PS_METADATA_LOOKUP_TYPE_DECL(Vector, psVector*);
 PS_METADATA_LOOKUP_TYPE_DECL(Metadata, psMetadata*);
 PS_METADATA_LOOKUP_TYPE_DECL(Time, psTime*);
Index: /trunk/psModules/src/camera/pmFPAMosaic.c
===================================================================
--- /trunk/psModules/src/camera/pmFPAMosaic.c	(revision 25382)
+++ /trunk/psModules/src/camera/pmFPAMosaic.c	(revision 25383)
@@ -740,5 +740,5 @@
 static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                        psImage **mosaicMask, // The mosaic mask, to be returned
-                       psImage **mosaicVariances, // The mosaic variances, to be returned
+                       psImage **mosaicVariance, // The mosaic variance, to be returned
                        int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
                        const pmChip *chip, // Chip to mosaic
@@ -749,5 +749,5 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinChip);
     assert(yBinChip);
@@ -826,5 +826,5 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank);
     }
@@ -847,5 +847,5 @@
 static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                       psImage **mosaicMask, // The mosaic mask, to be returned
-                      psImage **mosaicVariances, // The mosaic variances, to be returned
+                      psImage **mosaicVariance, // The mosaic variance, to be returned
                       int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
                       const pmFPA *fpa,  // FPA to mosaic
@@ -857,5 +857,5 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinFPA);
     assert(yBinFPA);
@@ -960,5 +960,5 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank);
     }
@@ -1025,5 +1025,5 @@
     psImage *mosaicImage   = NULL;      // The mosaic image
     psImage *mosaicMask    = NULL;      // The mosaic mask
-    psImage *mosaicVariances = NULL;      // The mosaic variances
+    psImage *mosaicVariance = NULL;      // The mosaic variances
 
     // Find the HDU
@@ -1052,6 +1052,6 @@
         }
         if (hdu->variances) {
-            mosaicVariances = psImageSubset(hdu->variances->data[0], bounds);
-            if (!mosaicVariances) {
+            mosaicVariance = psImageSubset(hdu->variances->data[0], bounds);
+            if (!mosaicVariance) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n");
                 return false;
@@ -1061,5 +1061,5 @@
         // Case 2 --- we need to mosaic by cut and paste
         psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) {
+        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariance, &xBin, &yBin, source, targetCell, blank)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
             return false;
@@ -1069,4 +1069,5 @@
     }
     psTrace("psModules.camera", 1, "xBin,yBin: %d,%d\n", xBin, yBin);
+
 
     // Set the concepts for the target cell
@@ -1090,9 +1091,31 @@
     target->parent->concepts = psMetadataCopy(target->parent->concepts, source->parent->concepts); // FPA lvl
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < source->cells->n; i++) {
+        pmCell *cell = source->cells->data[i]; // Cell of interest
+        if (!cell || !cell->data_exists) {
+            continue;
+        }
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+        if (!ro || !ro->covariance) {
+            continue;
+        }
+        psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+    }
+    psKernel *mosaicCovariance = NULL;  // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariance = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the target cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
-    newReadout->variance = mosaicVariances;
+    newReadout->variance = mosaicVariance;
+    newReadout->covariance = mosaicCovariance;
     psFree(newReadout);                 // Drop reference
 
@@ -1334,4 +1357,32 @@
     target->concepts = psMetadataCopy(target->concepts, source->concepts);
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < covariances->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chip || !chip->data_exists) {
+            continue;
+        }
+        psArray *cells = chip->cells;   // Cells in chip
+        for (long j = 0; j < cells->n; j++) {
+            pmCell *cell = cells->data[i]; // Cell of interest
+            if (!cell || !cell->data_exists) {
+                continue;
+            }
+            pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+            if (!ro || !ro->covariance) {
+                continue;
+            }
+            psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+        }
+    }
+    psKernel *mosaicCovariances = NULL; // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariances = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the new cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
@@ -1339,4 +1390,5 @@
     newReadout->mask   = mosaicMask;
     newReadout->variance = mosaicVariances;
+    newReadout->covariance = mosaicCovariances;
     psFree(newReadout);                 // Drop reference
 
Index: /trunk/psModules/src/camera/pmReadoutFake.c
===================================================================
--- /trunk/psModules/src/camera/pmReadoutFake.c	(revision 25382)
+++ /trunk/psModules/src/camera/pmReadoutFake.c	(revision 25383)
@@ -82,8 +82,4 @@
     }
     PS_ASSERT_PTR_NON_NULL(psf, false);
-    if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius.");
-        return false;
-    }
 
     readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
@@ -128,6 +124,11 @@
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
         fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
-        float fakeRadius = radius > 0 ? radius :
-            PS_MAX(1.0, fakeModel->modelRadius(fakeModel->params, minFlux)); // Radius of fake source
+        float fakeRadius = 1.0;         // Radius of fake source
+        if (isfinite(minFlux)) {
+            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
+        }
+        if (radius > 0) {
+            fakeRadius = PS_MAX(fakeRadius, radius);
+        }
 
         if (xOffset) {
Index: /trunk/psModules/src/camera/pmReadoutFake.h
===================================================================
--- /trunk/psModules/src/camera/pmReadoutFake.h	(revision 25382)
+++ /trunk/psModules/src/camera/pmReadoutFake.h	(revision 25383)
@@ -13,6 +13,20 @@
 #include <pmSourceMasks.h>
 
+/// Generate a fake readout from vectors
+bool pmReadoutFakeFromVectors(pmReadout *readout, ///< Output readout
+                              int numCols, int numRows, ///< Dimension of image
+                              const psVector *x, const psVector *y, ///< Source coordinates
+                              const psVector *mag, ///< Source magnitudes
+                              const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
+                              const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL
+                              const pmPSF *psf, ///< PSF for sources
+                              float minFlux, ///< Minimum flux to bother about; for setting source radius
+                              int radius, ///< Fixed radius for sources
+                              bool circularise, ///< Circularise PSF model?
+                              bool normalisePeak ///< Normalise the peak value?
+    );
+
 /// Generate a fake readout from an array of sources
-bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL
+bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout
                               int numCols, int numRows, ///< Dimension of image
                               const psArray *sources, ///< Array of pmSource
Index: /trunk/psModules/src/objects/Makefile.am
===================================================================
--- /trunk/psModules/src/objects/Makefile.am	(revision 25382)
+++ /trunk/psModules/src/objects/Makefile.am	(revision 25383)
@@ -53,5 +53,6 @@
 	pmGrowthCurveGenerate.c \
 	pmGrowthCurve.c \
-	pmSourceMatch.c
+	pmSourceMatch.c \
+	pmDetEff.c
 
 EXTRA_DIST = \
@@ -90,5 +91,6 @@
 	pmTrend2D.h \
 	pmGrowthCurve.h \
-	pmSourceMatch.h
+	pmSourceMatch.h \
+	pmDetEff.h
 
 CLEANFILES = *~
Index: /trunk/psModules/src/objects/pmDetEff.c
===================================================================
--- /trunk/psModules/src/objects/pmDetEff.c	(revision 25383)
+++ /trunk/psModules/src/objects/pmDetEff.c	(revision 25383)
@@ -0,0 +1,168 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <pslib.h>
+
+#include "pmDetEff.h"
+
+
+static void detEffFree(pmDetEff *de)
+{
+    psFree(de->magOffsets);
+    psFree(de->counts);
+    psFree(de->magDiffMean);
+    psFree(de->magDiffStdev);
+    psFree(de->magErrMean);
+}
+
+
+pmDetEff *pmDetEffAlloc(float magRef, int numSources, int numBins)
+{
+    pmDetEff *de = psAlloc(sizeof(pmDetEff)); // Detection efficiency, to return
+    psMemSetDeallocator(de, (psFreeFunc)detEffFree);
+
+    de->magRef = magRef;
+    de->numSources = numSources;
+    de->numBins = numBins;
+
+    de->magOffsets = NULL;
+    de->counts = NULL;
+    de->magDiffMean = NULL;
+    de->magDiffStdev = NULL;
+    de->magErrMean = NULL;
+
+    return de;
+}
+
+
+bool pmDetEffWrite(psFits *fits, pmDetEff *de, const psMetadata *header, const char *extname)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    PM_ASSERT_DETEFF_RESULTS(de, false);
+
+    psArray *table = psArrayAlloc(de->numBins); // Table to write
+    for (int i = 0; i < de->numBins; i++) {
+        psMetadata *row = psMetadataAlloc(); // Table row
+        psMetadataAddF32(row, PS_LIST_TAIL, "OFFSET", 0, "Magnitude offset from reference",
+                         de->magOffsets->data.F32[i]);
+        psMetadataAddS32(row, PS_LIST_TAIL, "COUNTS", 0, "Number of sources recovered",
+                         de->counts->data.S32[i]);
+        psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.MEAN", 0, "Mean magnitude difference",
+                         de->magDiffMean->data.F32[i]);
+        psMetadataAddF32(row, PS_LIST_TAIL, "DIFF.STDEV", 0, "Stdev magnitude difference",
+                         de->magDiffStdev->data.F32[i]);
+        psMetadataAddF32(row, PS_LIST_TAIL, "ERR.MEAN", 0, "Mean magnitude error",
+                         de->magErrMean->data.F32[i]);
+        table->data[i] = row;
+    }
+
+    psMetadata *deHeader = psMetadataCopy(NULL, header); // Header for detection efficiency
+    psMetadataAddF32(deHeader, PS_LIST_TAIL, "DETEFF.MAGREF", PS_META_REPLACE, "Magnitude reference",
+                     de->magRef);
+    psMetadataAddS32(deHeader, PS_LIST_TAIL, "DETEFF.NUM", PS_META_REPLACE, "Number of fake sources injected",
+                     de->numSources);
+
+    if (!psFitsWriteTable(fits, deHeader, table, extname)) {
+        psError(PS_ERR_IO, false, "Unable to write detection efficiency table.");
+        psFree(table);
+        psFree(deHeader);
+        return false;
+    }
+
+    psFree(table);
+    psFree(deHeader);
+
+    return true;
+}
+
+bool pmReadoutWriteDetEff(psFits *fits, const pmReadout *readout,
+                          const psMetadata *header, const char *extname)
+{
+    PM_ASSERT_READOUT_NON_NULL(readout, false);
+
+    bool mdok;                          // Status of MD lookup
+    pmDetEff *de = psMetadataLookupPtr(&mdok, readout->analysis, PM_DETEFF_ANALYSIS); // Detection efficiency
+    if (!mdok || !de) {
+        // Wrote everything there was to write
+        return true;
+    }
+    return pmDetEffWrite(fits, de, header, extname);
+}
+
+
+pmDetEff *pmDetEffRead(psFits *fits, const char *extname)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_STRING_NON_EMPTY(extname, false);
+
+    if (!psFitsMoveExtNameClean(fits, extname)) {
+        // Nothing to read
+        return NULL;
+    }
+
+    psMetadata *header = psFitsReadHeader(NULL, fits); // Header for table
+    if (!header) {
+        psError(PS_ERR_IO, false, "Unable to read FITS header");
+        return NULL;
+    }
+
+    int numBins = psFitsTableSize(fits);// Size of table
+    bool mdok;                          // Status of MD lookup
+    int numSources = psMetadataLookupS32(&mdok, header, "DETEFF.NUM"); // Number of fake sources injected
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find number of sources");
+        psFree(header);
+        return NULL;
+    }
+    float magRef = psMetadataLookupF32(&mdok, header, "DETEFF.MAGREF"); // Magnitude reference
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find magnitude reference");
+        psFree(header);
+        return NULL;
+    }
+    psFree(header);
+
+    pmDetEff *de = pmDetEffAlloc(magRef, numSources, numBins); // Detection efficiency
+    de->magOffsets = psVectorAlloc(numBins, PS_TYPE_F32);
+    de->counts = psVectorAlloc(numBins, PS_TYPE_S32);
+    de->magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32);
+    de->magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32);
+    de->magErrMean = psVectorAlloc(numBins, PS_TYPE_F32);
+
+    psArray *table = psFitsReadTable(fits); // FITS table
+    if (!table) {
+        psError(PS_ERR_IO, false, "Unable to read detection efficiency table.");
+        psFree(de);
+        return false;
+    }
+
+    for (int i = 0; i < numBins; i++) {
+        psMetadata *row = table->data[i]; // Table row
+        de->magOffsets->data.F32[i] = psMetadataLookupF32(NULL, row, "OFFSET");
+        de->counts->data.S32[i] = psMetadataLookupS32(NULL, row, "COUNTS");
+        de->magDiffMean->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.MEAN");
+        de->magDiffStdev->data.F32[i] = psMetadataLookupF32(NULL, row, "DIFF.STDEV");
+        de->magErrMean->data.F32[i] = psMetadataLookupF32(NULL, row, "ERR.MEAN");
+    }
+
+    psFree(table);
+    return de;
+}
+
+bool pmReadoutReadDetEff(psFits *fits, const pmReadout *readout, const char *extname)
+{
+    PM_ASSERT_READOUT_NON_NULL(readout, false);
+
+    pmDetEff *de = pmDetEffRead(fits, extname);
+    if (!de) {
+        if (psErrorCodeLast() != PS_ERR_NONE) {
+            return false;
+        }
+        return true;
+    }
+
+    return psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PM_DETEFF_ANALYSIS,
+                            PS_META_REPLACE | PS_DATA_UNKNOWN, "Detection efficiency", de);
+}
Index: /trunk/psModules/src/objects/pmDetEff.h
===================================================================
--- /trunk/psModules/src/objects/pmDetEff.h	(revision 25383)
+++ /trunk/psModules/src/objects/pmDetEff.h	(revision 25383)
@@ -0,0 +1,81 @@
+#ifndef PM_DET_EFF_H
+#define PM_DET_EFF_H
+
+#include <pslib.h>
+#include <string.h>
+
+#include "pmFPA.h"
+
+#define PM_DETEFF_ANALYSIS "DETEFF"     // Location of detection efficiency on pmReadout.analysis
+
+// Detection efficiency characterisation
+typedef struct {
+    float magRef;                       // Reference magnitude
+    int numSources;                     // Number of sources
+    int numBins;                        // Number of bins
+    psVector *magOffsets;               // Magnitude offsets for each bin
+    psVector *counts;                   // Counts of sources recovered for each bin
+    psVector *magDiffMean;              // Mean magnitude difference for each bin
+    psVector *magDiffStdev;             // Stdev of magnitude difference for each bin
+    psVector *magErrMean;               // Mean magnitude error for each bin
+} pmDetEff;
+
+
+/// Allocator
+pmDetEff *pmDetEffAlloc(float magRef,   // Reference magnitude
+                        int numSources, // Number of sources
+                        int numBins     // Number of bins
+                        );
+
+/// Write detection efficiency to FITS file
+bool pmDetEffWrite(psFits *fits,        // FITS file to which to write
+                   pmDetEff *deteff,    // Detection efficiency to write
+                   const psMetadata *header, // Header to write
+                   const char *extname  // Extension name
+                   );
+
+/// Write detection efficiency from a readout to a FITS file
+bool pmReadoutWriteDetEff(psFits *fits,// FITS file to which to write
+                          const pmReadout *readout, // Readout with detection efficiency
+                          const psMetadata *header, // Header to write
+                          const char *extname // Extension name
+    );
+
+/// Read detection efficiency
+pmDetEff *pmDetEffRead(psFits *fits,    // FITS file from which to read
+                       const char *extname // Extension name
+                       );
+
+/// Read detection efficiency into a readout
+bool pmReadoutReadDetEff(psFits *fits,// FITS file to which to write
+                         const pmReadout *readout, // Readout with detection efficiency
+                         const char *extname // Extension name
+    );
+
+#define PM_ASSERT_DETEFF_NON_NULL(DE, RETURN) { \
+    if (!(DE)) { \
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Detection efficiency %s is NULL", #DE); \
+        return RETURN; \
+    } \
+}
+
+#define PM_ASSERT_DETEFF_RESULTS(DE, RETURN) { \
+    PM_ASSERT_DETEFF_NON_NULL(DE, RETURN); \
+    PS_ASSERT_VECTOR_NON_NULL((DE)->magOffsets, RETURN); \
+    PS_ASSERT_VECTOR_SIZE((DE)->magOffsets, (long)(DE)->numBins, RETURN); \
+    PS_ASSERT_VECTOR_TYPE((DE)->magOffsets, PS_TYPE_F32, RETURN); \
+    PS_ASSERT_VECTOR_NON_NULL((DE)->counts, RETURN); \
+    PS_ASSERT_VECTOR_SIZE((DE)->counts, (long)(DE)->numBins, RETURN); \
+    PS_ASSERT_VECTOR_TYPE((DE)->counts, PS_TYPE_S32, RETURN); \
+    PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffMean, RETURN); \
+    PS_ASSERT_VECTOR_SIZE((DE)->magDiffMean, (long)(DE)->numBins, RETURN); \
+    PS_ASSERT_VECTOR_TYPE((DE)->magDiffMean, PS_TYPE_F32, RETURN); \
+    PS_ASSERT_VECTOR_NON_NULL((DE)->magDiffStdev, RETURN); \
+    PS_ASSERT_VECTOR_SIZE((DE)->magDiffStdev, (long)(DE)->numBins, RETURN); \
+    PS_ASSERT_VECTOR_TYPE((DE)->magDiffStdev, PS_TYPE_F32, RETURN); \
+    PS_ASSERT_VECTOR_NON_NULL((DE)->magErrMean, RETURN); \
+    PS_ASSERT_VECTOR_SIZE((DE)->magErrMean, (long)(DE)->numBins, RETURN); \
+    PS_ASSERT_VECTOR_TYPE((DE)->magErrMean, PS_TYPE_F32, RETURN); \
+}
+
+#endif
Index: /trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO.c	(revision 25382)
+++ /trunk/psModules/src/objects/pmSourceIO.c	(revision 25383)
@@ -42,8 +42,82 @@
 #include "pmSource.h"
 #include "pmModelClass.h"
+#include "pmDetEff.h"
 #include "pmSourceIO.h"
 
 #define BLANK_HEADERS "BLANK.HEADERS"   // Name of metadata in camera configuration containing header names
                                         // for putting values into a blank PHU
+
+// lookup the EXTNAME values used for table data and image header segments
+static bool sourceExtensions(psString *headname, // Extension name for header
+                             psString *dataname, // Extension name for data
+                             psString *deteffname, // Extension name for detection efficiency
+                             psString *xsrcname, // Extension name for extended sources
+                             psString *xfitname, // Extension name for extended fits
+                             const pmFPAfile *file, // File of interest
+                             const pmFPAview *view // View to level of interest
+                             )
+{
+    bool status;                        // Status of MD lookup
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+        psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+        return false;
+    }
+
+    // EXTNAME for image header
+    if (headname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *headname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for table data
+    if (dataname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *dataname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for detection efficiency
+    if (deteffname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.DETEFF");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DETEFF in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *deteffname = pmFPAfileNameFromRule(rule, file, view);
+    }
+
+    // EXTNAME for extended source data table
+    if (xsrcname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xsrcname = pmFPAfileNameFromRule (rule, file, view);
+    }
+
+    if (xfitname) {
+        // EXTNAME for extended source data table
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xfitname = pmFPAfileNameFromRule (rule, file, view);
+    }
+
+    return true;
+}
+
 
 // translations between psphot object types and dophot object types
@@ -271,8 +345,4 @@
 
     char *exttype  = NULL;
-    char *dataname = NULL;
-    char *xsrcname = NULL;
-    char *xfitname = NULL;
-    char *headname = NULL;
 
     // if sources is NULL, write out an empty table
@@ -354,49 +424,12 @@
 
         // define the EXTNAME values for the different data segments:
-        {
-            // lookup the EXTNAME values used for table data and image header segments
-            char *rule = NULL;
-
-            // Menu of EXTNAME rules
-            psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-            if (!menu) {
-                psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-                return false;
-            }
-
-            // EXTNAME for image header
-            rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
-            if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
-                return false;
-            }
-            headname = pmFPAfileNameFromRule (rule, file, view);
-
-            // EXTNAME for table data
-            rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
-            if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
-                return false;
-            }
-            dataname = pmFPAfileNameFromRule (rule, file, view);
-
-            if (XSRC_OUTPUT) {
-              // EXTNAME for extended source data table
-              rule = psMetadataLookupStr(&status, menu, "CMF.XSRC");
-              if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XSRC in EXTNAME.RULES in camera.config");
-                return false;
-              }
-              xsrcname = pmFPAfileNameFromRule (rule, file, view);
-            }
-            if (XFIT_OUTPUT) {
-              // EXTNAME for extended source data table
-              rule = psMetadataLookupStr(&status, menu, "CMF.XFIT");
-              if (!rule) {
-                psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XFIT in EXTNAME.RULES in camera.config");
-                return false;
-              }
-              xfitname = pmFPAfileNameFromRule (rule, file, view);
-            }
+        psString headname = NULL;
+        psString dataname = NULL;
+        psString deteffname = NULL;
+        psString xsrcname = NULL;
+        psString xfitname = NULL;
+        if (!sourceExtensions(&headname, &dataname, &deteffname, XSRC_OUTPUT ? &xsrcname : NULL,
+                              XFIT_OUTPUT ? &xfitname : NULL, file, view)) {
+            return false;
         }
 
@@ -480,49 +513,54 @@
 
             // XXX these are case-sensitive since the EXTYPE is case-sensitive
-            status = false;
+            status = true;
             if (!strcmp (exttype, "SMPDATA")) {
-                status = pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_SMPDATA (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_DEV_0")) {
-                status = pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_DEV_0 (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_DEV_1")) {
-                status = pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_DEV_1 (file->fits, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_CAL_0")) {
-                status = pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_PS1_CAL_0 (file->fits, readout, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_V1")) {
-                status = pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
+                status &= pmSourcesWrite_CMF_PS1_V1 (file->fits, readout, sources, file->header, outhead, dataname);
             }
             if (!strcmp (exttype, "PS1_V2")) {
-                status = pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
-            }
+                status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
+            }
+
+            if (deteffname) {
+                status &= pmReadoutWriteDetEff(file->fits, readout, outhead, deteffname);
+            }
+
             if (xsrcname) {
               if (!strcmp (exttype, "PS1_DEV_1")) {
-                  status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_CAL_0")) {
-                  status = pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
+                  status &= pmSourcesWrite_PS1_CAL_0_XSRC (file->fits, readout, sources, file->header, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_V1")) {
-                  status = pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_CMF_PS1_V1_XSRC (file->fits, sources, xsrcname, recipe);
               }
               if (!strcmp (exttype, "PS1_V2")) {
-                  status = pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
+                  status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
               }
             }
             if (xfitname) {
               if (!strcmp (exttype, "PS1_DEV_1")) {
-                  status = pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_PS1_DEV_1_XFIT (file->fits, sources, xfitname);
               }
               if (!strcmp (exttype, "PS1_CAL_0")) {
-                  status = pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
+                  status &= pmSourcesWrite_PS1_CAL_0_XFIT (file->fits, readout, sources, file->header, xfitname);
               }
               if (!strcmp (exttype, "PS1_V1")) {
-                  status = pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_CMF_PS1_V1_XFIT (file->fits, sources, xfitname);
               }
               if (!strcmp (exttype, "PS1_V2")) {
-                  status = pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
+                  status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
               }
             }
@@ -572,6 +610,6 @@
     // not needed if only one chip
     if (file->fpa->chips->n == 1) {
-	pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
-	return true;
+        pmSourceIO_WriteMatchedRefs (file->fits, file->fpa, config);
+        return true;
     }
 
@@ -885,26 +923,11 @@
         hdu = pmFPAviewThisHDU (view, file->fpa);
 
-        // lookup the EXTNAME values used for table data and image header segments
-        char *rule = NULL;
-        // Menu of EXTNAME rules
-        psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-        if (!menu) {
-            psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-            return false;
-        }
-        // EXTNAME for image header
-        rule = psMetadataLookupStr(&status, menu, "CMF.HEAD");
-        if (!rule) {
-            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.HEAD in EXTNAME.RULES in camera.config");
-            return false;
-        }
-        char *headname = pmFPAfileNameFromRule (rule, file, view);
-        // EXTNAME for table data
-        rule = psMetadataLookupStr(&status, menu, "CMF.DATA");
-        if (!rule) {
-            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.DATA in EXTNAME.RULES in camera.config");
-            return false;
-        }
-        char *dataname = pmFPAfileNameFromRule (rule, file, view);
+        // define the EXTNAME values for the different data segments:
+        psString headname = NULL;
+        psString dataname = NULL;
+        psString deteffname = NULL;
+        if (!sourceExtensions(&headname, &dataname, &deteffname, NULL, NULL, file, view)) {
+            return false;
+        }
 
         // advance to the IMAGE HEADER extension
@@ -958,4 +981,9 @@
                 sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
             }
+
+            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
+                psError(PS_ERR_IO, false, "Unable to read detection efficiency");
+                return false;
+            }
         }
 
@@ -1070,3 +1098,3 @@
 }
 
-    
+
Index: /trunk/psModules/src/psmodules.h
===================================================================
--- /trunk/psModules/src/psmodules.h	(revision 25382)
+++ /trunk/psModules/src/psmodules.h	(revision 25383)
@@ -133,4 +133,5 @@
 #include <pmSourceVisual.h>
 #include <pmSourceMatch.h>
+#include <pmDetEff.h>
 
 // The following headers are from random locations, here because they cross bounds
Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 25382)
+++ /trunk/psphot/src/Makefile.am	(revision 25383)
@@ -128,5 +128,6 @@
 	psphotCheckStarDistribution.c  \
 	psphotThreadTools.c  	       \
-	psphotAddNoise.c
+	psphotAddNoise.c               \
+	psphotEfficiency.c
 
 # dropped? psphotGrowthCurve.c
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 25382)
+++ /trunk/psphot/src/psphot.h	(revision 25383)
@@ -49,9 +49,9 @@
 bool            psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
-bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
+bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 
 bool            psphotReplaceAllSources (psArray *sources, psMetadata *recipe);
-bool            psphotRemoveAllSources (psArray *sources, psMetadata *recipe);
+bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
 
 bool            psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
@@ -64,5 +64,5 @@
 bool            psphotGuessModel_Threaded (psThreadJob *job);
 
-bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf);
+bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf);
 bool            psphotMagnitudes_Threaded (psThreadJob *job);
 
@@ -76,4 +76,5 @@
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
+bool            psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources);
 
 // thread-related:
@@ -147,7 +148,7 @@
 pmPSF          *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
 bool            psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
-bool            psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe);
-bool            psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add);
+bool            psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
+bool            psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
+bool            psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add);
 bool            psphotRadialPlot (int *kapa, const char *filename, pmSource *source);
 bool            psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);
Index: /trunk/psphot/src/psphotAddNoise.c
===================================================================
--- /trunk/psphot/src/psphotAddNoise.c	(revision 25382)
+++ /trunk/psphot/src/psphotAddNoise.c	(revision 25383)
@@ -1,13 +1,13 @@
 # include "psphotInternal.h"
 
-bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
   return psphotAddOrSubNoise (readout, sources, recipe, true);
 }
 
-bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
   return psphotAddOrSubNoise (readout, sources, recipe, false);
 }
 
-bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) {
+bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) {
 
     bool status = false;
@@ -42,5 +42,5 @@
     PS_ASSERT (status, false);
     if (isfinite(GAIN)) {
-	FACTOR /= GAIN;
+        FACTOR /= GAIN;
     }
 
@@ -70,7 +70,7 @@
         oldshape.sxy = PAR[PM_PAR_SXY];
 
-	// XXX can this be done more intelligently?
-	if (oldI0 == 0.0) continue;
-	if (!isfinite(oldI0)) continue;
+        // XXX can this be done more intelligently?
+        if (oldI0 == 0.0) continue;
+        if (!isfinite(oldI0)) continue;
 
         // increase size and height of source
Index: /trunk/psphot/src/psphotEfficiency.c
===================================================================
--- /trunk/psphot/src/psphotEfficiency.c	(revision 25383)
+++ /trunk/psphot/src/psphotEfficiency.c	(revision 25383)
@@ -0,0 +1,476 @@
+# include "psphotInternal.h"
+
+#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
+                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
+
+//#define TESTING
+
+
+// Calculate the limiting magnitude for an image
+//
+// We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo
+// non-Gaussian PSF) to the limiting flux in the un-smoothed original image.
+static bool effLimit(float *magLim,           // Limiting magntiude, to return
+                     int *radius,             // Radius for fake sources, to return
+                     float *minFlux,          // Minimum flux for fake sources, to return
+                     float *norm,             // Normalisation of PSF (conversion: peak --> integrated flux)
+                     float *covarFactor,// Covariance factor
+                     const pmReadout *ro,     // Readout of interest
+                     const pmPSF *psf,        // Point-spread function
+                     float thresh,            // Threshold for source identification
+                     float smoothSigma,       // Gaussian smoothing sigma
+                     float smoothNsigma,      // Smoothing limit
+                     psImageMaskType maskVal  // Value to mask
+                     )
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_VARIANCE(ro, false);
+    assert(magLim);
+    assert(radius);
+    assert(minFlux);
+    assert(norm);
+    assert(covarFactor);
+
+    psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
+    psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix
+    psFree(kernel);
+    *covarFactor = psImageCovarianceFactor(newCovar); // Variance factor
+    psFree(newCovar);
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);        // Random number generator
+    if (!psImageBackground(stats, NULL, ro->variance, ro->mask, maskVal, rng) ||
+        !isfinite(stats->robustMedian)) {
+        psError(PSPHOT_ERR_DATA, false, "Unable to determine mean variance");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    psFree(rng);
+    float meanVar = stats->robustMedian; // Mean variance
+    psFree(stats);
+
+    // Need to normalise out difference between Gaussian and real PSF
+    pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols,
+                                             0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+    if (!normModel || (normModel->flags & MODEL_MASK)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model.");
+        psFree(normModel);
+        return false;
+    }
+    *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
+    psFree(normModel);
+
+    // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
+    // since the variance factor tells us the variance in the smoothed image.  Now, the trick is working
+    // out what the intensity in the smoothed image is, and how it is related to the flux.  We are
+    // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star,
+    // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the
+    // unsmoothed image, and w is the width of the Gaussian.  Now, a normalised 2D Gaussian convolved
+    // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2.  Therefore:
+    // I_smooth = Flux / 2*norm.
+    float peakLim = thresh * sqrtf(meanVar * *covarFactor); // Limiting peak value in smoothed image
+    float fluxLim = 2.0 * *norm * peakLim; // Limiting flux in original
+    *magLim = -2.5 * log10f(fluxLim);
+    psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim);
+    psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim);
+    psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim);
+
+    *radius = smoothSigma * smoothNsigma;
+
+    *minFlux = 0.1 * sqrtf(meanVar);
+
+    return true;
+}
+
+/// Generate a fake image and add it in to the existing readout
+static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
+                        const pmReadout *ro,            // Readout of interest
+                        const pmPSF *psf,               // Point-spread function
+                        const psVector *magOffsets,     // Magnitude offsets for fake sources
+                        int numSources,                 // Number of fake sources for each bin
+                        float refMag,                   // Reference magnitude
+                        int radius,                     // Radius for fake sources
+                        float minFlux                   // Minimum flux for fake sources
+                        )
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+    PS_ASSERT_VECTOR_NON_NULL(magOffsets, false);
+    PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false);
+    assert(xSrc);
+    assert(ySrc);
+
+    int numBins = magOffsets->n;                                    // Number of bins
+    int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
+
+    *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32);
+    *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32);
+
+    psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+    psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+    psVector *magAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    for (int i = 0, index = 0; i < numBins; i++) {
+        float mag = refMag + magOffsets->data.F32[i]; // Instrumental magnitude of sources
+
+        for (int j = 0; j < numSources; j++, index++) {
+            xAll->data.F32[index] = psRandomUniform(rng) * numCols;
+            yAll->data.F32[index] = psRandomUniform(rng) * numRows;
+            magAll->data.F32[index] = mag;
+        }
+        memcpy((*xSrc)->data.F32[i], &xAll->data.F32[i * numSources],
+               numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy((*ySrc)->data.F32[i], &yAll->data.F32[i * numSources],
+               numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+    }
+    psFree(rng);
+
+    pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
+    if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll,
+                                  NULL, NULL, psf, minFlux, radius, false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image");
+        psFree(fakeRO);
+        psFree(magAll);
+        return false;
+    }
+    psFree(magAll);
+
+    psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
+    psFree(fakeRO);
+
+    return true;
+}
+
+
+// Determine detection efficiency
+bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf,
+                psMetadata *recipe, const psArray *realSources)
+{
+    PM_ASSERT_READOUT_NON_NULL(readout, false);
+    PM_ASSERT_READOUT_IMAGE(readout, false);
+    PS_ASSERT_PTR_NON_NULL(psf, false);
+    PS_ASSERT_METADATA_NON_NULL(recipe, false);
+    PS_ASSERT_ARRAY_NON_NULL(realSources, false);
+
+    psTimerStart("psphot.fake");
+
+    // Collect recipe information
+    float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x
+    float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y
+    if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe");
+        return false;
+    }
+    float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
+    if (!isfinite(smoothNsigma)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe");
+        return false;
+    }
+    float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2");
+    if (!isfinite(thresh)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe");
+        return false;
+    }
+    psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets
+    if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe");
+        return NULL;
+    }
+    int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin
+    if (numSources == 0) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe");
+        return NULL;
+    }
+    float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
+    if (!isfinite(minGauss)) {
+        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
+        minGauss = 0.5;
+    }
+
+    float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma
+    int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
+    int numBins = magOffsets->n;                          // Number of bins
+
+    psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
+
+    // remove all sources, adding noise for subtracted sources
+    psphotRemoveAllSources(realSources, recipe);
+    //    psphotAddNoise(readout, realSources, recipe);
+
+
+#if defined(TESTING) && 0
+    {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                readout->image->data.F32[y][x] = psRandomGaussian(rng);
+                readout->variance->data.F32[y][x] = 1.0;
+                readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+            }
+        }
+        psFree(readout->covariance);
+        readout->covariance = NULL;
+        psFree(rng);
+    }
+#endif
+
+
+
+    float magLim;                       // Guess at limiting magnitude
+    int radius;                         // Radius for fake sources
+    float minFlux;                      // Minimum flux for fake sources
+    float norm;                         // Normalisation of PSF
+    float covarFactor;                  // Covariance factor
+    if (!effLimit(&magLim, &radius, &minFlux, &norm, &covarFactor, readout,
+                  psf, thresh, smoothSigma, smoothNsigma, maskVal)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image");
+        return false;
+    }
+
+#ifdef TESTING
+    psphotSaveImage(NULL, readout->image, "orig_image.fits");
+    psphotSaveImage(NULL, readout->variance, "orig_variance.fits");
+    psphotSaveImage(NULL, readout->mask, "orig_mask.fits");
+#endif
+
+    psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
+    if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets,
+                     numSources, magLim, radius, minFlux)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
+        psFree(xFake);
+        psFree(yFake);
+        return false;
+    }
+
+#ifdef TESTING
+    psphotSaveImage(NULL, readout->image, "fake_image.fits");
+    psphotSaveImage(NULL, readout->variance, "fake_variance.fits");
+    psphotSaveImage(NULL, readout->mask, "fake_mask.fits");
+#endif
+
+    // XXX Could speed this up significantly by only convolving the central pixels of each fake source
+    psVector *significance = NULL;       // Significance image
+    {
+        int num = numSources * numBins; // Total number of sources
+        // Pixel coordinates of sources
+        psVector *xPix = psVectorAlloc(num, PS_TYPE_S32);
+        psVector *yPix = psVectorAlloc(num, PS_TYPE_S32);
+        for (int i = 0, index = 0; i < numBins; i++) {
+            for (int j = 0; j < numSources; j++, index++) {
+                float x = xFake->data.F32[i][j];
+                float y = yFake->data.F32[i][j];
+                xPix->data.S32[index] = PS_MIN(x + 0.5, numCols - 1);
+                yPix->data.S32[index] = PS_MIN(y + 0.5, numRows - 1);
+            }
+        }
+
+        psVector *convImage = psImageSmoothMaskPixels(readout->image, readout->mask, maskVal, xPix, yPix,
+                                                      smoothSigma, smoothNsigma, minGauss); // Convolved image
+        if (!convImage) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth image pixels");
+            psFree(xPix);
+            psFree(yPix);
+            psFree(xFake);
+            psFree(yFake);
+            return false;
+        }
+        psVector *convVar = psImageSmoothMaskPixels(readout->variance, readout->mask, maskVal, xPix, yPix,
+                                                    smoothSigma, smoothNsigma, minGauss); // Convolved varianc
+        if (!convVar) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth variance pixels");
+            psFree(xPix);
+            psFree(yPix);
+            psFree(xFake);
+            psFree(yFake);
+            return false;
+        }
+
+        float factor = 1.0 / covarFactor; // Correction for covariance
+
+        const psImage *mask = readout->mask;  // Mask for readout
+        for (int i = 0; i < num; i++) {
+            float imageVal = convImage->data.F32[i]; // Image value
+            float varVal = convVar->data.F32[i]; // Variance value
+            int x = xPix->data.S32[i], y = yPix->data.S32[i]; // Coordinates
+            if (imageVal < 0 || varVal <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
+                convImage->data.F32[i] = 0.0;
+            } else {
+                convImage->data.F32[i] = factor * PS_SQR(imageVal) / varVal;
+            }
+        }
+
+        psFree(xPix);
+        psFree(yPix);
+
+        significance = convImage;
+        psFree(convVar);
+    }
+
+    thresh *= thresh;                   // "Significance" is actually significance-squared
+
+    psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin
+    psArray *fakeSources = psArrayAlloc(numSources);            // Fake sources in each bin
+    psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources
+    for (int i = 0, index = 0; i < numBins; i++) {
+        int numFound = 0;               // Number found
+
+        // Determine extraction size
+        float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin
+        float peak = powf(10.0, -0.4 * mag) / norm;   // Peak flux
+        pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source
+        if (!model || (model->flags & MODEL_MASK)) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to generate model for bin %d", i);
+            psFree(model);
+            return false;
+        }
+        float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source
+        psFree(model);
+
+        psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin
+        for (int j = 0; j < numSources; j++, index++) {
+            // Coordinates of interest
+            float sig = significance->data.F32[index]; // Significance of pixel
+            if (sig > thresh) {
+                pmSource *source = pmSourceAlloc(); // Fake source
+                float x = xFake->data.F32[i][j];
+                float y = yFake->data.F32[i][j];
+                source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE);
+                if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) {
+                    psErrorClear();
+                    continue;
+                }
+                source->peak->xf = x;
+                source->peak->yf = y;
+                source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig));
+                source->type = PM_SOURCE_TYPE_STAR;
+
+                numFound++;
+                psArrayAdd(sources, sources->n, source);
+                psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source);
+                psFree(source);
+            }
+        }
+        fakeSources->data[i] = sources;
+        count->data.S32[i] = numFound;
+    }
+    psFree(xFake);
+    psFree(yFake);
+    psFree(significance);
+
+    if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
+        psFree(fakeSources);
+        psFree(count);
+        return false;
+    }
+
+    // Disable aperture corrections; casting away const!
+    pmTrend2D *apTrend = psf->ApTrend;  // Aperture trend
+    ((pmPSF*)psf)->ApTrend = NULL;
+
+    if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
+        psFree(fakeSources);
+        psFree(count);
+        ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const!
+        return false;
+    }
+    psFree(fakeSourcesAll);
+
+    // Re-enable aperture corrections; casting away const!
+    ((pmPSF*)psf)->ApTrend = apTrend;
+
+    psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin
+    psVector *magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of diff in magnitude for each bin
+    psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin
+    psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32);   // Magnitude differences
+    psVector *magErr = psVectorAlloc(numSources, PS_TYPE_F32);   // Magnitude errors
+    psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors
+    psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV |
+                                     PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    for (int i = 0; i < numBins; i++) {
+        psStatsInit(magStats);
+        psVectorInit(magMask, 0);
+
+#ifdef TESTING
+        psString name = NULL;
+        psStringAppend(&name, "fake_%d.dat", i);
+        FILE *file = fopen(name, "w");
+        psFree(name);
+#endif
+
+        float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude
+        psArray *sources = fakeSources->data[i];         // Sources in bin
+        for (int j = 0; j < sources->n; j++) {
+            pmSource *source = sources->data[j]; // Source of interest
+            if (!source || !isfinite(source->psfMag)) {
+                magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;
+                continue;
+            }
+
+#ifdef TESTING
+            fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf,
+                    source->modelPSF->params->data.F32[PM_PAR_XPOS],
+                    source->modelPSF->params->data.F32[PM_PAR_YPOS],
+                    magRef, source->psfMag, source->errMag);
+#endif
+            magDiff->data.F32[j] = source->psfMag - magRef;
+            magErr->data.F32[j] = source->errMag;
+        }
+        magDiff->n = sources->n;
+        magMask->n = sources->n;
+        magErr->n = sources->n;
+        if (!psVectorStats(magStats, magDiff, NULL, magMask, 0xFF)) {
+            // Probably because we don't have enough sources
+            psErrorClear();
+        }
+
+        if (isfinite(magStats->robustMedian) && isfinite(magStats->robustStdev)) {
+            magDiffMean->data.F32[i] = magStats->robustMedian;
+            magDiffStdev->data.F32[i] = magStats->robustStdev;
+        } else {
+            magDiffMean->data.F32[i] = magStats->sampleMean;
+            magDiffStdev->data.F32[i] = magStats->sampleStdev;
+        }
+
+        if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) {
+            // Probably because we don't have enough sources
+            psErrorClear();
+        }
+
+        if (isfinite(magStats->robustMedian)) {
+            magErrMean->data.F32[i] = magStats->robustMedian;
+        } else {
+            magErrMean->data.F32[i] = magStats->sampleMean;
+        }
+
+#ifdef TESTING
+        fclose(file);
+#endif
+    }
+
+    psFree(magStats);
+    psFree(magDiff);
+    psFree(magMask);
+    psFree(magErr);
+
+    psFree(fakeSources);
+
+    pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency
+    de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32);
+    de->counts = count;
+    de->magDiffMean = magDiffMean;
+    de->magDiffStdev = magDiffStdev;
+    de->magErrMean = magErrMean;
+
+    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PM_DETEFF_ANALYSIS, PS_META_REPLACE | PS_DATA_UNKNOWN,
+                     "Detection efficiency", de);
+    psFree(de);
+
+    psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
+
+
+    return true;
+}
Index: /trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 25382)
+++ /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 25383)
@@ -12,5 +12,5 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
+bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
 
     bool status;
Index: /trunk/psphot/src/psphotMagnitudes.c
===================================================================
--- /trunk/psphot/src/psphotMagnitudes.c	(revision 25382)
+++ /trunk/psphot/src/psphotMagnitudes.c	(revision 25383)
@@ -1,5 +1,5 @@
 # include "psphotInternal.h"
 
-bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
+bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) {
 
     bool status = false;
@@ -64,5 +64,5 @@
 
             psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, (pmPSF*)psf);    // Casting away const
             psArrayAdd(job->args, 1, binning);
             psArrayAdd(job->args, 1, backModel);
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 25382)
+++ /trunk/psphot/src/psphotReadout.c	(revision 25383)
@@ -61,5 +61,5 @@
     // display the backsub and backgnd images
     psphotVisualShowBackground (config, view, readout);
-    
+
     // run a single-model test if desired (exits from here if test is run)
     psphotModelTest (config, view, recipe);
@@ -224,4 +224,9 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
+    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
+
     // replace failed sources?
     // psphotReplaceUnfitSources (sources);
Index: /trunk/psphot/src/psphotReadoutMinimal.c
===================================================================
--- /trunk/psphot/src/psphotReadoutMinimal.c	(revision 25382)
+++ /trunk/psphot/src/psphotReadoutMinimal.c	(revision 25383)
@@ -97,4 +97,9 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
+    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
+
     // drop the references to the image pixels held by each source
     psphotSourceFreePixels (sources);
Index: /trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /trunk/psphot/src/psphotReplaceUnfit.c	(revision 25382)
+++ /trunk/psphot/src/psphotReplaceUnfit.c	(revision 25383)
@@ -47,5 +47,5 @@
 }
 
-bool psphotRemoveAllSources (psArray *sources, psMetadata *recipe) {
+bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
 
     bool status;
Index: /trunk/psphot/src/psphotSignificanceImage.c
===================================================================
--- /trunk/psphot/src/psphotSignificanceImage.c	(revision 25382)
+++ /trunk/psphot/src/psphotSignificanceImage.c	(revision 25383)
@@ -22,20 +22,23 @@
     }
 
-    bool status_x, status_y;
-    float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
-    float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
-    if (status_x && status_y) {
-      // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
-      SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
-      NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
-      guess = false;
+    bool statusMajor, statusMinor;
+    float fwhmMajor = psMetadataLookupF32(&statusMajor, recipe, "FWHM_MAJ");
+    float fwhmMinor = psMetadataLookupF32(&statusMinor, recipe, "FWHM_MIN");
+    if (statusMajor && statusMinor) {
+        // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
+        if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
+            psWarning("fwhmMajor (%f) or fwhmMinor (%f) is bad!", fwhmMajor, fwhmMinor);
+        }
+        SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        guess = false;
     } else {
-      // if we do not know the FWHM, use the guess smoothing kernel supplied.
-      // it is a configuration error if these are not supplied
-      SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
-      PS_ASSERT (status, NULL);
-      NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
-      PS_ASSERT (status, NULL);
-      guess = true;
+        // if we do not know the FWHM, use the guess smoothing kernel supplied.
+        // it is a configuration error if these are not supplied
+        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
+        PS_ASSERT (status, NULL);
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        PS_ASSERT (status, NULL);
+        guess = true;
     }
     // record the actual smoothing sigma