Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 19486)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 19532)
@@ -15,5 +15,5 @@
                      PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 
-//#define TESTING                         // Enable debugging output
+#define TESTING                         // Enable debugging output
 
 
@@ -48,5 +48,5 @@
 #endif
 
-bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2,
                   const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
 {
@@ -111,5 +111,5 @@
                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
             float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(output->parent->concepts, "CELL.VARFACTOR");
+            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
             if (!isfinite(vf)) {
                 vf = 1.0;
@@ -341,4 +341,25 @@
     assert((*regions)->n == (*kernels)->n);
 
+    // Correct the variance for the chi^2
+    {
+        *chi2 = 0.0;
+        int num = 0;                    // Number of measurements of chi^2
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
+            *chi2 += kernels->mean;
+            num++;
+        }
+        psFree(iter);
+
+        float vf = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.VARFACTOR"); // Variance factor
+        *chi2 /= vf * num;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {