Index: trunk/ippconfig/recipes/ppStack.config
===================================================================
--- trunk/ippconfig/recipes/ppStack.config	(revision 20496)
+++ trunk/ippconfig/recipes/ppStack.config	(revision 20497)
@@ -3,4 +3,5 @@
 ITER		S32	1		# Number of rejection iterations
 COMBINE.REJ	F32	3.0		# Rejection threshold in combination (sigma)
+COMBINE.SYS	F32	0.08		# Relative systematic error in combination
 MASK.VAL	STR	MASK.VALUE,BAD.WARP	# Mask value of input bad pixels
 MASK.BAD	STR	BLANK		# Mask value to give bad pixels
Index: trunk/ppStack/src/ppStackArguments.c
===================================================================
--- trunk/ppStack/src/ppStackArguments.c	(revision 20496)
+++ trunk/ppStack/src/ppStackArguments.c	(revision 20497)
@@ -145,4 +145,6 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0,
                      "Combination rejection thresold (sigma)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-sys", 0,
+                     "Relative systematic error in combination", NAN);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-val", 0, "Mask value of input bad pixels", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Mask value to give bad pixels", NULL);
@@ -226,4 +228,5 @@
     VALUE_ARG_RECIPE_INT("-iter",             "ITER",           S32, 0);
     VALUE_ARG_RECIPE_FLOAT("-combine-rej",    "COMBINE.REJ",    F32);
+    VALUE_ARG_RECIPE_FLOAT("-combine-sys",    "COMBINE.SYS",    F32);
     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
     VALUE_ARG_RECIPE_FLOAT("-image-rej",      "IMAGE.REJ",      F32);
Index: trunk/ppStack/src/ppStackReadout.c
===================================================================
--- trunk/ppStack/src/ppStackReadout.c	(revision 20496)
+++ trunk/ppStack/src/ppStackReadout.c	(revision 20497)
@@ -96,4 +96,5 @@
     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
+    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
     bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
     bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
@@ -133,6 +134,6 @@
     }
 
-    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, true,
-                        useVariance, safe)) {
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys,
+                        true, useVariance, safe)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
         psFree(stack);
Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 20496)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 20497)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-31 21:50:59 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-11-01 02:59:33 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -219,4 +219,5 @@
                           int numIter, // Number of rejection iterations
                           float rej, // Number of standard deviations at which to reject
+                          float sys,    // Relative systematic error
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
@@ -311,11 +312,17 @@
           if (useVariance && numIter > 0) {
               // Use variance to check that the two are consistent
-              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
+              float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
+              float var1 = pixelVariances->data.F32[0]; // Variance of first
+              float var2 = pixelVariances->data.F32[1]; // Variance of second
 #ifdef VARIANCE_FACTORS
-              float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] +
-                  pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]];
-#else
-              float sigma2 = pixelVariances->data.F32[0] + pixelVariances->data.F32[1];
+              // Variance factor contributes to the rejection level
+              var1 *= varFactors->data.F32[pixelSources->data.U16[0]];
+              var2 *= varFactors->data.F32[pixelSources->data.U16[1]];
 #endif
+              // Systematic error contributes to the rejection level
+              var1 += PS_SQR(sys * pixelData->data.F32[0]);
+              var2 += PS_SQR(sys * pixelData->data.F32[1]);
+
+              float sigma2 = var1 + var2; //
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
@@ -346,9 +353,11 @@
                   float rej2 = PS_SQR(rej); // Rejection level squared
                   for (int i = 0; i < num; i++) {
+                      pixelVariances->data.F32[i] *= rej2;
 #ifdef VARIANCE_FACTORS
-                      pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]];
-#else
-                      pixelVariances->data.F32[i] *= rej2;
+                      // Variance factor contributes to the rejection level
+                      pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]];
 #endif
+                      // Systematic error contributes to the rejection level
+                      pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]);
                   }
               }
@@ -554,5 +563,6 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad,
-                    int kernelSize, int numIter, float rej, bool entire, bool useVariance, bool safe)
+                    int kernelSize, int numIter, float rej, float sys,
+                    bool entire, bool useVariance, bool safe)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
@@ -660,5 +670,5 @@
                                                      x, y); // Reject these images
                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                                  reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                                  reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
                 }
             }
@@ -674,5 +684,5 @@
                                                  x, y); // Reject these images
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }
@@ -701,5 +711,5 @@
             for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              NULL, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }
Index: trunk/psModules/src/imcombine/pmStack.h
===================================================================
--- trunk/psModules/src/imcombine/pmStack.h	(revision 20496)
+++ trunk/psModules/src/imcombine/pmStack.h	(revision 20497)
@@ -8,6 +8,6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-09-12 04:12:42 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-11-01 02:59:33 $
  *
  * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -48,4 +48,5 @@
                     int numIter,        ///< Number of iterations
                     float rej,          ///< Rejection limit (standard deviations)
+                    float sys,          ///< Relative systematic error
                     bool entire,        ///< Combine entire image even if rejection lists provided?
                     bool useVariance,   ///< Use variance values for rejection?