Index: /branches/pap/EFFICIENCY.txt
===================================================================
--- /branches/pap/EFFICIENCY.txt	(revision 25062)
+++ /branches/pap/EFFICIENCY.txt	(revision 25062)
@@ -0,0 +1,51 @@
+Strategy for efficiency analysis (a.k.a. "fake" stage)
+
+Want:
+* Detection efficiency as a function of instrumental magnitude
+* Masked fraction
+
+Given:
+* Image
+* Mask
+* Variance
+* PSF
+* Sources on image
+* Recipe:
+  + Magnitude bins (relative to guessed detection limit) for fake sources
+  + Number of fake sources for each bin
+
+
+
+Algorithm:
+* Determine mean instrumental magnitude detection limit
+  + Have:
+    - Mean variance
+    - Smoothing size
+    - Covariance
+    - Threshold
+  + Calculate mean peak flux of source at threshold
+  + Integrate PSF to determine magnitude
+* Remove all real sources from image
+* Add fake sources into image
+* Smooth image, variance (maybe only at positions of fake sources?)
+* Count number of sources masked
+* Count number of fake sources above significance level as a function of fake source magnitude
+
+
+If the required density of fake sources is so high that they begin to
+overlap, we will need to do add sources, smooth, and count
+independently for each magnitude bin.  This could be optimised
+using the distributive property of convolutions (f*[g+h]=f*g+f*h)
+if we smooth the image first, and then add in smoothed PSFs.
+
+
+For a 5% statistical measurement of the detection efficiency per bin, we want
+dN/N = 0.05 ==> N ~ 400 for Poisson statistics.
+Assume a 80% masked fraction, so want 500 sources per bin.
+For 5 magnitude bins, want 2500 sources
+Assume worst case PSF is 15 pix FWHM --> ~ 51 pixels for 4-sigma either side
+Total area is 6502500 pixels.
+Chip images are 4846x4868 = 23590328 pixels = 3.6 times the worst total area.
+Not too bad (especially since this is pretty much worst case, and we only
+care about the central pixel, which is much smaller)
+==> Let's just throw all sources on the same image.
Index: /branches/pap/psphot/src/psphotFake.c
===================================================================
--- /branches/pap/psphot/src/psphotFake.c	(revision 25061)
+++ /branches/pap/psphot/src/psphotFake.c	(revision 25062)
@@ -9,5 +9,5 @@
 
 
-// Calculate the limiting flux for an image
+// Calculate the limiting magnitude for an image
 //
 // We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo
@@ -62,6 +62,67 @@
     }
 
-    return thresh * sqrtf(meanVar * factor); // Limiting peak flux in smoothed image
+    return -2.5 * log10(thresh * sqrtf(meanVar * factor));
 }
+
+static bool fakeGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
+                         const pmReadout *ro,            // Readout of interest
+                         const psMetadata *recipe, // psphot recipe
+                         float refMag              // Reference magnitude
+                         )
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro, NULL);
+    PS_ASSERT_METADATA_NON_NULL(recipe, NULL);
+
+    psMetadata *fake = psMetadataLookupMetadata(NULL, recipe, "FAKE"); // Fake information
+    if (!fake) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE information in recipe");
+        return NULL;
+    }
+
+    psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "FAKE.MAG"); // Magnitude offsets
+    if (!magOffset || magOffset->type.type != PS_TYPE_F32) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE.MAG F32 vector in recipe");
+        return NULL;
+    }
+    int numBins = magOffset->n;                                     // Number of bins
+    int numSources = psMetadataLookupS32(NULL, recipe, "FAKE.NUM"); // Number of sources for each bin
+
+    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);
+
+    // Master list, for image creation
+    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; i <= numBins; i++) {
+        psArray *data = psArrayAlloc(3); // Sources for bin
+
+        psVector *x = psVectorAlloc(numSources, PS_TYPE_F32);
+        psVector *y = psVectorAlloc(numSources, PS_TYPE_F32);
+
+        float mag = refMag + magOffset->data.F32[i]; // Instrumental magnitude of sources
+
+        for (int j = 0; j <= numSources; j++) {
+            x->data.F32[j] = psRandomUniform(rng) * numCols;
+            y->data.F32[j] = psRandomUniform(rng) * numRows;
+            magAll->data.F32[i * numSources + j] = mag;
+        }
+        memcpy(xSrc->data.F32[i], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(ySrc->data.F32[i], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(&xAll->data.F32[i * numSources], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(&yAll->data.F32[i * numSources], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+    }
+    psFree(rng);
+
+    pmReadout *fakeRO = pmReadoutAlloc(); // Fake readout
+    pmReadoutFakeFromVectors(fakeRO, xAll, yAll, mag,
+
 
 
@@ -80,5 +141,5 @@
 
     // Generate array of fake sources and define the pixels with pmSourceDefinePixels
-    psArray *fakeOrig = psphotFakeGenerate(config, readout, recipe);
+    psArray *fakeOrig = fakeGenerate(config, readout, recipe);
 
     psArray *fakeMeasured = psArrayAlloc(fakeOrig->n); // Measured values for the sources in fakeOrig