Index: trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- trunk/Ohana/src/relphot/include/relphot.h	(revision 39640)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 39641)
@@ -276,4 +276,5 @@
 int    PARALLEL_MANUAL;
 int    PARALLEL_SERIAL;
+char  *MANUAL_UNIQUER;
 
 int    PARALLEL_REGIONS_MANUAL;
@@ -286,5 +287,6 @@
 int    MOSAIC_ZEROPT;
 
-int    TEST_IMAGE;
+int    TEST_IMAGE1;
+int    TEST_IMAGE2;
 
 int    NLOOP;
@@ -309,4 +311,5 @@
 int    FREEZE_IMAGES;
 int    FREEZE_MOSAICS;
+int    CALIBRATE_STACKS_AND_WARPS;
 int    USE_GRID;
 int    KEEP_UBERCAL;
Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 39641)
@@ -1,3 +1,5 @@
 # include "relphot.h"
+# define BASIC_STATS 0
+
 void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 
@@ -614,7 +616,7 @@
 
   off_t i, j, m, c, n, N, Nmax;
-  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys, Nbright;
+  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys;
   float Msys, Mrel, Mmos, Mgrid, Mflat;
-  double *list, *dlist, *Mlist, *dMlist;
+  // double *list, *dlist, *Mlist, *dMlist;
 
   StatType stats;
@@ -637,8 +639,6 @@
     Nmax = MAX (Nmax, N_onImage[i]);
   }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-  ALLOCATE (Mlist, double, Nmax);
-  ALLOCATE (dMlist, double, Nmax);
+  StatDataSet *refStars = StatDataSetAlloc (1, Nmax);
+  StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
 
   Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
@@ -647,4 +647,11 @@
   for (i = 0; i < Nimage; i++) {
     
+    if (image[i].imageID == TEST_IMAGE1) {
+      fprintf (stderr, "test image 1\n");
+    }
+    if (image[i].imageID == TEST_IMAGE2) {
+      fprintf (stderr, "test image 1\n");
+    }
+
     /* on PoorImages run, skip good images */
     if (PoorImages) {
@@ -664,9 +671,8 @@
 
     int minUbercalDist = 1000;
-
-    // number of stars to measure the bright-end scatter
-    Nbright = 0;
-
-    N = 0;
+    
+    off_t Nref = 0;  // total number of reference stars on the image
+    int Nbright = 0; // number of stars to measure the bright-end scatter
+
     for (j = 0; j < N_onImage[i]; j++) {
       
@@ -716,20 +722,24 @@
 
     skip:
-      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
-      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      if (fabs(list[N]) > 0.03) {
+      refStars->flxlist[Nref] = Msys - Mrel - Mmos - Mgrid + Mflat;
+      refStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+      refStars->wgtlist[Nref] = 1;
+      refStars->msklist[Nref] = 0;
+      if (fabs(refStars->flxlist[Nref]) > 0.03) {
 	// fprintf (stderr, "deviant\n");
       }
 
-      if (image[i].imageID == TEST_IMAGE) {
-	fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, list[N]);
+      if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
+	fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, refStars->flxlist[Nref]);
       }
 
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-	Mlist[Nbright] = list[N];
-	dMlist[Nbright] = dlist[N];
+	brightStars->flxlist[Nbright] = refStars->flxlist[Nref];
+	brightStars->errlist[Nbright] = refStars->errlist[Nref];
+	brightStars->wgtlist[Nbright] = 1;
+	brightStars->msklist[Nbright] = 0;
 	Nbright ++;
       }
-      N++;
+      Nref++;
     }
     /* N_onImage[i] is all measurements, N is good measurements */
@@ -740,5 +750,5 @@
     mark = FALSE;
     if (!PoorImages) {
-      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onImage[i]);
+      mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]);
       if (mark) {
 	image[i].flags |= ID_IMAGE_PHOTOM_FEW;
@@ -750,12 +760,22 @@
 
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-    liststats (list, dlist, NULL, N, &stats);
-    image[i].Mcal  = stats.mean;
-    image[i].dMcal = stats.error;
-    image[i].nFitPhotom = N;
-    image[i].Xm    = 100.0*log10(stats.chisq);
+# if (BASIC_STATS) 
+    liststats (refStars->flxlist, refStars->errlist, NULL, Nref, &stats);
+# else
+    liststats_irls (refStars, Nref, &stats);
+# endif
+    image[i].Mcal       = stats.mean;
+    image[i].dMcal      = stats.error;
+    image[i].nFitPhotom = Nref;
+    image[i].Xm         = 100.0*log10(stats.chisq);
     Ncalibrated ++;
 
-    if (image[i].imageID == TEST_IMAGE) {
+    if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
+      for (j = 0; j < Nref; j++) {
+	fprintf (stderr, "%1d, %8d : %6.3f  %6.3f  %6.3f  %d\n", (int) i, (int) image[i].imageID, refStars->flxlist[j], refStars->errlist[j], refStars->wgtlist[j], refStars->msklist[j]);
+      }
+    }
+
+    if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
       fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
     }
@@ -767,9 +787,9 @@
     if (PLOTSTUFF) {
       fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
-      plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+      plot_setMcal (refStars->flxlist, N, &stats, CLOUD_TOLERANCE);
     }
 
     // bright end scatter
-    liststats (Mlist, dMlist, NULL, Nbright, &stats);
+    liststats (brightStars->flxlist, brightStars->errlist, NULL, Nbright, &stats);
     image[i].dMagSys = stats.sigma;
 
@@ -782,8 +802,7 @@
     image[i].ubercalDist = minUbercalDist + 1;
   }
-  free (list);
-  free (dlist);
-  free (Mlist);
-  free (dMlist);
+
+  StatDataSetFree (brightStars, 1);
+  StatDataSetFree (refStars, 1);
 
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);
Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 39641)
@@ -128,4 +128,12 @@
   }
 
+  // elements needed for parallel regions / parallel images
+  MANUAL_UNIQUER = NULL;
+  if ((N = get_argument (argc, argv, "-manual-uniquer"))) {
+    remove_argument (N, &argc, argv);
+    MANUAL_UNIQUER = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   PLOTSTUFF = FALSE;
   if ((N = get_argument (argc, argv, "-plot"))) {
@@ -176,8 +184,18 @@
   }
 
-  TEST_IMAGE = -1;
-  if ((N = get_argument (argc, argv, "-test-image"))) {
-    remove_argument (N, &argc, argv);
-    TEST_IMAGE = atoi(argv[N]);
+  TEST_IMAGE1 = -1;
+  if ((N = get_argument (argc, argv, "-test-image1"))) {
+    char *endptr;
+    remove_argument (N, &argc, argv);
+    TEST_IMAGE1 = strtol(argv[N], &endptr, 0);
+    if (*endptr) relphot_usage (); 
+    remove_argument (N, &argc, argv);
+  }
+  TEST_IMAGE2 = -1;
+  if ((N = get_argument (argc, argv, "-test-image2"))) {
+    char *endptr;
+    remove_argument (N, &argc, argv);
+    TEST_IMAGE2 = strtol(argv[N], &endptr, 0);
+    if (*endptr) relphot_usage (); 
     remove_argument (N, &argc, argv);
   }
@@ -283,9 +301,13 @@
     FREEZE_IMAGES = TRUE;
   }
-
   FREEZE_MOSAICS = FALSE;
   if ((N = get_argument (argc, argv, "-mosfreeze"))) {
     remove_argument (N, &argc, argv);
     FREEZE_MOSAICS = TRUE;
+  }
+  CALIBRATE_STACKS_AND_WARPS = FALSE;
+  if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_STACKS_AND_WARPS = TRUE;
   }
 
@@ -620,5 +642,6 @@
   }
 
-  TEST_IMAGE = -1;
+  TEST_IMAGE1 = -1;
+  TEST_IMAGE2 = -1;
 
   RESET = FALSE;
Index: trunk/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/bcatalog.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/bcatalog.c	(revision 39641)
@@ -106,5 +106,6 @@
 
       // skip garbage measurements
-      if (catalog[0].measure[offset].psfQF < 0.85) { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF < 0.85))     { Npsfqf ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < 0.85)) { Npsfqf ++; continue; }
       if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
 
Index: trunk/Ohana/src/relphot/src/liststats.c
===================================================================
--- trunk/Ohana/src/relphot/src/liststats.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/liststats.c	(revision 39641)
@@ -210,5 +210,6 @@
   
   // XXX add to dataset elements?
-  ALLOCATE_PTR (wt, double, Npoints);
+  ALLOCATE_PTR (wtvals, double, Npoints);
+  ALLOCATE_PTR (wtlist, double, Npoints);
 
   int converged = FALSE;
@@ -217,9 +218,9 @@
     for (int i = 0; i < Npoints; i++) {
       // we are only including the formal error, not the weight in the definition of wt[]
-      wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+      wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     }
     
     double oldValue = value;
-    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wt, Npoints)) {
+    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wtvals, Npoints)) {
       value = oldValue;
       break;
@@ -235,10 +236,15 @@
       
   // calculate the weight thresholds to mask the bad points:
-  double Sum_W = 0.0;
+  // double Sum_W = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
-    Sum_W += wt[i];
-  }
-  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
+    wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+    wtlist[i] = wtvals[i];
+    // Sum_W += wtvals[i];
+  }
+  dsort (wtlist, Npoints);
+  int midpt = 0.5 * Npoints;
+  double WtMedian = (Npoints % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
+  // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
+  double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
 
   // generate the unmasked subset
@@ -255,5 +261,5 @@
   double dSig = 0.0;
   for (int i = 0; i < Npoints; i++) {
-    if ((wt[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
+    if ((wtvals[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
       dataset->msklist[i] = TRUE; // mark the masked points
       continue;
@@ -313,5 +319,6 @@
   free (dykeep);
   free (wtkeep);
-  free (wt);
+  free (wtvals);
+  free (wtlist);
 
   return TRUE;
Index: trunk/Ohana/src/relphot/src/load_catalogs.c
===================================================================
--- trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 39641)
@@ -115,5 +115,10 @@
   int TIME = time(NULL);
   int PID = getpid();
-  snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
+
+  if (PARALLEL_MANUAL && MANUAL_UNIQUER) {
+    snprintf (uniquer, 12, "%11s", MANUAL_UNIQUER);
+  } else {
+    snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
+  }
 
   int Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
Index: trunk/Ohana/src/relphot/src/relphot_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/relphot_images.c	(revision 39640)
+++ trunk/Ohana/src/relphot/src/relphot_images.c	(revision 39641)
@@ -53,5 +53,5 @@
   FlatCorrectionTable *flatcorr = NULL;
 
-  if (NLOOP > 0) {
+  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     /* load catalog data from region files (hostID is 0 since we are not a client */
     catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
@@ -154,6 +154,7 @@
 
     /* set Mcal & Mmos for bad images */
-    setMcal  (catalog, TRUE, flatcorr);
-    setMmos  (catalog, TRUE, flatcorr);
+    int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
+    setMcal  (catalog, onlyPoorImages, flatcorr);
+    setMmos  (catalog, onlyPoorImages, flatcorr);
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
@@ -248,6 +249,8 @@
      context, this function writes the image parameters as a subset table for the remote
      clients */
-  reload_catalogs (skylist, flatcorr, 0, NULL);
-  MARKTIME("-- updated all catalogs: %f sec\n", dtime);
+  if (!CALIBRATE_STACKS_AND_WARPS) {
+    reload_catalogs (skylist, flatcorr, 0, NULL);
+    MARKTIME("-- updated all catalogs: %f sec\n", dtime);
+  }
 
   if (UPDATE_CATFORMAT) {
@@ -255,5 +258,5 @@
     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
   }
-  if (NLOOP > 0) {
+  if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
     // do not save changes if we did not make changes. 
     dvo_image_update (&db, VERBOSE);