Changeset 41556

Timestamp:

Apr 27, 2021, 8:39:17 AM (5 years ago)

Author:

eugene

Message:

updates to the good/bad night logic; save image stats

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (2 diffs)
• include/relphot.h (modified) (3 diffs)
• src/ImageOps.c (modified) (11 diffs)
• src/MagResidSave.c (added)
• src/MosaicOps.c (modified) (26 diffs)
• src/StarOps.c (modified) (1 diff)
• src/TGroupOps.c (modified) (31 diffs)
• src/relphot_images.c (modified) (6 diffs)
• src/relphot_parallel_images.c (modified) (2 diffs)
• src/relphot_parallel_regions.c (modified) (2 diffs)
• src/save_images.c (added)

trunk/Ohana/src/relphot/Makefile

@@ -54 +54 @@
 $(SRC)/load_catalogs.$(ARCH).o   \
 $(SRC)/load_images.$(ARCH).o     \
+$(SRC)/save_images.$(ARCH).o     \
+$(SRC)/MagResidSave.$(ARCH).o    \
 $(SRC)/synthetic_mags.$(ARCH).o  \
 $(SRC)/plot_scatter.$(ARCH).o    \
@@ -117 +119 @@
 $(SRC)/select_images.$(ARCH).o   \
 $(SRC)/load_images.$(ARCH).o     \
+$(SRC)/save_images.$(ARCH).o     \
+$(SRC)/MagResidSave.$(ARCH).o    \
 $(SRC)/load_catalogs.$(ARCH).o   \
 $(SRC)/reload_catalogs.$(ARCH).o \

trunk/Ohana/src/relphot/include/relphot.h

@@ -643 +643 @@
 int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
 
-void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
-int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
-int           setMgrp             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
+void          setMcal             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+int           setMmos             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+int           setMgrp             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
 
 int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
@@ -841 +841 @@
 void RepairWarpMeasuresCloseLogfile ();
 
-int save_images_backup (FITS_DB *db);
-
 void dump_tgroups (Catalog *catalog, int Npass);
 void dump_catalog (Catalog *catalog, off_t c, int Npass);
@@ -853 +851 @@
 int UseStandardOLS (ZptFitModeType mode);
 int GetActivePhotcodeIndex (int photcode);
+
+int save_images_backup (FITS_DB *db);
+int save_images_updates (FITS_DB *db);
+
+int MagResidSave(char *filename, Catalog *catalog);

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -75 +75 @@
 Image *getimage (off_t N) {
   return (&image[N]);
+}
+
+off_t *get_N_onImage() {
+  return N_onImage;
+}
+
+IDX_T **get_ImageToCatalog() {
+  return ImageToCatalog;
+}
+
+IDX_T **get_ImageToMeasure() {
+  return ImageToMeasure;
 }
 
@@ -570 +582 @@
 
 /* determine Mcal values for all images */
-void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+void setMcal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
   off_t i, j, m, c, n;
-  int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
+  int mark, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
 
   if (IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_NONE) return;
@@ -580 +592 @@
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  if (PoorImages) {
-    // XXX use bad stars and measurements for PoorImages? or not?
-    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
-    IMAGE_BAD = 0;
-  }
 
   off_t Nmax = 0;
@@ -618 +624 @@
     }
 
-    /* on PoorImages run, skip good images */
-    if (PoorImages) {
-      bad = image[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
-      if (!bad) continue;
-    }      
-
     if (image[i].photcode == 0) continue; // skip the PHU images
 
@@ -629 +629 @@
     int badMosaic = (image[i].flags & ID_IMAGE_MOSAIC_POOR);
     int isMosaic  = isMosaicChip(image[i].photcode);
+    int useMcal   = TRUE;
 
     // unset at start by default (set if actually used below)
@@ -638 +639 @@
     // in BAD_NIGHT mode, we fit ONLY images in bad nights
     if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT)) {
-      if (!badNight) continue;
+      if (!badNight) useMcal = FALSE; // do not fit images from good nights
     }
 
     // in BAD_NIGHT_BAD_MOSAIC mode, we fit ONLY bad mosaics in bad nights
     if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC)) {
-      if (!badNight) continue;
-      if (!badMosaic && isMosaic) continue; // skip good mosaics (but not good images)
+      if (!badNight) useMcal = FALSE; // do not fit images from good nights
+      if (!badMosaic && isMosaic) useMcal = FALSE; // skip good mosaics (but not good images)
     }
 
     // in BAD_MOSAIC mode, we fit bad mosaics ignoring night state
     if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_MOSAIC)) {
-      if (!badMosaic && isMosaic) continue; // skip good mosaics (but not bad images)
+      if (!badMosaic && isMosaic) useMcal = FALSE; // skip good mosaics (but not bad images)
     }
 
@@ -759 +760 @@
 
     /* too few good measurements or too many bad measurements */
-    mark = FALSE;
-    if (!PoorImages) {
-      mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]);
-      if (mark) {
-        image[i].flags |= ID_IMAGE_PHOTOM_FEW;
-        Nfew ++;
-      } else {
-        image[i].flags &= ~ID_IMAGE_PHOTOM_FEW;
-      }      
-    }
+    mark = (Nref < IMAGE_TOOFEW) || (Nref < IMAGE_GOOD_FRACTION*N_onImage[i]);
+    if (mark) {
+      image[i].flags |= ID_IMAGE_PHOTOM_FEW;
+      image[i].dMcal      = NAN;
+      image[i].McalChiSq  = NAN;
+      image[i].nFitPhotom = 0;
+      Nfew ++;
+    } else {
+      image[i].flags &= ~ID_IMAGE_PHOTOM_FEW;
+    }      
     if (mark) continue;
 
+    // use liststats to find the 20-pct, median, 80-pct points
+    StatType stats;
+    liststats_setmode (&stats, "MEDIAN");
+    liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &stats);
+    double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+
+    // soften the individual errors with 10% of the scatter above
+    for (j = 0; j < Nref; j++) {
+      double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma);
+      psfStars.alldata->dyVector[j] = newSigma;
+    }
+
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
+    // XXX: apply airmass from above and fit only zp?
+
+    // XXX soften the stats?
+
     fit1d_irls (&psfStars, Nref);
-    image[i].McalPSF    = psfStars.bSaveArray[0][0];
+    fit1d_irls (&kronStars, Nkron);
+    fit1d_irls (&brightStars, Nbright);
+
+    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
+    if (useMcal) {
+      image[i].McalPSF    = psfStars.bSaveArray[0][0];
+      image[i].McalAPER   = kronStars.bSaveArray[0][0];
+      image[i].flags |= ID_IMAGE_IMAGE_PHOTCAL;  // set this flag (unset by default)
+    } else {
+      image[i].McalPSF    = 0.0;
+      image[i].McalAPER   = 0.0;
+    }
+
+    // record statistics on the fit regardless if it is used or not
     image[i].dMcal      = psfStars.bSigma[0];
     image[i].nFitPhotom = psfStars.Nmeas;
@@ -779 +809 @@
     Ncalibrated ++;
 
-    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-    // XXX: apply airmass from above and fit only zp?
-    // XXX EAM : temporarily use this to track calibration quality
-    if (0) {
-      fit1d_irls (&kronStars, Nkron);
-      image[i].McalAPER   = kronStars.bSaveArray[0][0];
-    }
+    // bright end scatter
+    image[i].dMagSys = brightStars.sigma;
 
     if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
@@ -796 +821 @@
     }
 
-    // bright end scatter
-    fit1d_irls (&brightStars, Nbright);
-    image[i].dMagSys = brightStars.sigma;
-
-    if (image[i].McalPSF < -CLOUD_TOLERANCE) {
-      image[i].McalPSF = 0.0;
-    }
-
-    image[i].flags |= ID_IMAGE_IMAGE_PHOTCAL;  // set this flag
-
     // minUbercalDist calculated here is the min value for any star owned by this image
     // since this particular image is tied to that star, bump its distance by 1
@@ -813 +828 @@
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);
   fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys);
-
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-    STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-  }
 
   FitDataSetFree (&brightStars);
@@ -1291 +1300 @@
   image = NULL;
 }
-
-void checkImages (char *name) {
-
-  off_t i;
-  fprintf (stderr, "--- %s ---\n", name);
-  
-  OhanaMemblock *ref = (OhanaMemblock *) image - 1;
-  fprintf (stderr, "file: %s\n", ref->file);
-
-  fprintf (stderr, "0x%08llx \n", (unsigned long long) image);
-
-  for (i = 0; i < Nimage; i++) {
-    fprintf (stderr, "%s : 0x%08x : %6.3f  %6.3f\n", image[i].name, image[i].flags, image[i].McalPSF, image[i].McalAPER);
-  }
-}

trunk/Ohana/src/relphot/src/MosaicOps.c

@@ -1 +1 @@
 # include "relphot.h"
 off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
+int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr);
 
 // see discussion in ImagesOps.c re: IDX_T
@@ -983 +984 @@
   int Nskip;
   off_t Nmax;
-  int PoorImages;
   FitDataSet psfStars;
   FitDataSet kronStars;
@@ -1002 +1002 @@
 int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
 void *setMmos_worker (void *data);
-int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
-
-void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
+int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr);
+
+void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists) {
   info->Nfew = 0;
   info->Nbad = 0;
@@ -1014 +1014 @@
 
   info->Nmax = Nmax;
-  info->PoorImages = PoorImages;
 
   if (allocLists) {
@@ -1069 +1068 @@
 }
 
-int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMmos (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
   off_t i, N, Nmax;
@@ -1082 +1081 @@
   // so do not run setMmos in threaded mode if PLOTSTUFF is set
   if (NTHREADS && !PLOTSTUFF) {
-    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
+    int status = setMmos_threaded (catalog, flatcorr);
     return status;
   }
@@ -1089 +1088 @@
 
   fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
-
-  if (PoorImages) {
-    // XXX use bad stars and measurements for PoorImages? or not?
-    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
-    IMAGE_BAD = 0;
-  }
 
   Nmax = 0;
@@ -1102 +1095 @@
 
   SetMmosInfo info;
-  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
+  SetMmosInfoInit (&info, Nmax, TRUE);
 
   // int savePlotDelay = PLOTDELAY;
@@ -1121 +1114 @@
   fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys);
 
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-    STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-  }
   return (TRUE);
 }
@@ -1134 +1122 @@
   off_t j, NimageReal;
 
-  /* on PoorImages run, skip good images */
-  if (info->PoorImages) {
-    int bad = myMosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
-    if (!bad) return TRUE;
-  }      
-
   FitDataSet *psfStars = &info->psfStars;
   FitDataSet *kronStars = &info->kronStars;
@@ -1147 +1129 @@
   assert (Nmos < Nmosaic);
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int minUbercalDist = 1000;
+
+  // unset this flag at start (set below if zeropoint is calculated)
+  myMosaic->flags &= ~ID_IMAGE_MOSAIC_PHOTCAL;  // unset this flag
+
+  // calculate the statistics for both good and bad exposures, but only set Mmos if the rules allow
+  // we only skip the statistics for nights with too few measurements or exposures
+
+  // if we are fitting TGroup zero points, when we identify the bad TGroups, we proceed
+  // to fit the mosaics which are from the bad TGroups
+
+  // if we are fitting Mosaic zero points, when we identify the bad nights, we proceed
+  // to fit the images which are from the bad Mosaics
+
   int badNight  = (myMosaic[0].flags & ID_IMAGE_NIGHT_POOR);
   int badMosaic = (myMosaic[0].flags & ID_IMAGE_MOSAIC_POOR);
-
-  myMosaic->flags &= ~ID_IMAGE_MOSAIC_PHOTCAL;  // unset this flag
-
-  // if we are fitting TGroup zero points, when we identify the bad TGroups, we proceed
-  // to fit the mosaics which are from the bad TGroups
-
-  // if we are fitting Mosaic zero points, when we identify the bad nights, we proceed
-  // to fit the images which are from the bad Mosaics
-
-  // in BAD_NIGHT mode, we fit ONLY mosaics in bad nights
+  int useMmos   = TRUE;
+
+  // in BAD_NIGHT mode, we fit ONLY mosaics in bad nights (do not fit mosaics from good nights)
   if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT) {
-    if (!badNight) return TRUE;
+    if (!badNight) useMmos = FALSE;
   }
 
   // in BAD_NIGHT_GOOD_MOSAIC mode, we fit ONLY good mosaics in bad nights
   if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC)) {
-    if (!badNight) return TRUE;
-    if (badMosaic) return TRUE;
+    if (!badNight) useMmos = FALSE; // do not fit good nights
+    if (badMosaic) useMmos = FALSE; // do not fit bad mosaics
   }
 
   // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state
   if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) {
-    if (badMosaic) return TRUE;
+    if (badMosaic) useMmos = FALSE; // do not fit bad mosaics
   }
 
@@ -1189 +1180 @@
       // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     }
+    // XXX if KEEP_UBERCAL, calculate stats but do not change Mmos
     if (KEEP_UBERCAL) return TRUE;
   }
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  int minUbercalDist = 1000;
 
   // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE)
@@ -1204 +1192 @@
   int testImage = FALSE;
   // testImage |= (abs(myMosaic[0].start - 1324104046) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1324103823) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003245) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003069) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003125) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003300) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003365) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003191) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003014) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003484) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323003419) < 10);
-  // testImage |= (abs(myMosaic[0].start - 1323002949) < 10);
 
   FILE *fout = NULL;
@@ -1221 +1198 @@
     snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
     fout = fopen (filename, "w");
+    save_test_mosaic_measures (fout, Nmos, catalog, flatcorr);
+    fclose (fout);
   }
 
@@ -1232 +1211 @@
     off_t c = MosaicToCatalog[Nmos][j];
       
-    // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF
-    if (fout) {
-      float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
-      // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-      float Mgrid    = getMgrid (m, c);
-      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      float Mflat    = getMflat (m, c, flatcorr, catalog);
-
-      off_t n = catalog[c].measureT[m].averef;
-      float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-
-      float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
-
-      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
-
-      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad);
-    }
-
     if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
       info->Nbad ++;
@@ -1329 +1290 @@
   /* N_onMosaic[Nmos] is all measurements, N is good measurements */
 
-  if (fout) { fclose (fout); }
-
-  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
-  if (!info->PoorImages) {
-    int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
-    if (mark) {
-      if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name,  Nmos,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[Nmos]); }
-      myMosaic[0].flags |= ID_IMAGE_PHOTOM_FEW;
-      info->Nfew ++;
-      if (testImage) {
-        fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
-      }
-      return TRUE; // skip mosaics with too few good measurements
-    } else {
-      myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
-    }
+  /* skip mosaics with too few good measurements or too many bad measurements */
+  int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
+  if (mark) {
+    myMosaic->flags |= ID_IMAGE_PHOTOM_FEW;
+    myMosaic->McalPSF    = 0.0;
+    myMosaic->McalAPER   = 0.0;
+    myMosaic->dMcal      = NAN;
+    myMosaic->stdev      = NAN;
+    myMosaic->dMmin      = NAN;
+    myMosaic->dMmax      = NAN;
+    myMosaic->McalChiSq  = NAN;
+    myMosaic->nFitPhotom = 0;
+    myMosaic->nValPhotom = N;
+    info->Nfew ++;
+    if (testImage) {
+      fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
+    }
+    return TRUE; // skip mosaics with too few good measurements
+  } else {
+    myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
+  }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEDIAN");
+  liststats (psfStars->alldata->yVector, NULL, NULL, N, &stats);
+  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+
+  // soften the individual errors with 10% of the scatter above
+  for (j = 0; j < N; j++) {
+    double newSigma = hypot(psfStars->alldata->dyVector[j], 0.1*altSigma);
+    psfStars->alldata->dyVector[j] = newSigma;
   }
 
   fit1d_irls (psfStars, N);
   fit1d_irls (brightStars, Nbright);
+  fit1d_irls (kronStars, N);
+
+  if (useMmos) {
+    myMosaic->McalPSF    = psfStars->bSaveArray[0][0];
+    myMosaic->McalAPER   = kronStars->bSaveArray[0][0];
+    myMosaic->flags |= ID_IMAGE_MOSAIC_PHOTCAL;  // set this flag
+  } else {
+    myMosaic->McalPSF    = 0.0;
+    myMosaic->McalAPER   = 0.0;
+  }
 
   // for now, I have no reason to measure these separately for camera-level images
-  myMosaic[0].McalPSF    = psfStars->bSaveArray[0][0];
-  myMosaic[0].dMcal      = psfStars->bSigma[0];     
-  myMosaic[0].stdev      = psfStars->sigma;
-  myMosaic[0].McalChiSq  = psfStars->chisq;         
-  myMosaic[0].nFitPhotom = psfStars->Nmeas;         
-  myMosaic[0].nValPhotom = N;       
-
-  fit1d_irls (kronStars, N);
-  myMosaic->McalAPER   = kronStars->bSaveArray[0][0];
+  myMosaic->dMcal      = psfStars->bSigma[0];       
+  myMosaic->stdev      = psfStars->sigma;
+  myMosaic->McalChiSq  = psfStars->chisq;           
+  myMosaic->nFitPhotom = psfStars->Nmeas;           
+  myMosaic->nValPhotom = N;         
+
+  // bright end scatter
+  myMosaic->dMsys = brightStars->sigma; // stdev of bright star mags
 
   if (testImage) {
     fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic->start, myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom);
+    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq);
   }
 
@@ -1368 +1354 @@
     fprintf (stderr, "Mmos: %6.3f +/- %6.3f %5d %5d | %s\n", myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom, N, image[MosaicToImage[Nmos][0]].name);
     plot_setMcal (psfStars->alldata->yVector, N);
-  }
-
-  // bright end scatter
-  fit1d_irls (brightStars, Nbright);
-  myMosaic->dMsys = brightStars->sigma; // stdev of bright star mags
-
-  if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) {
-    myMosaic[0].McalPSF = 0.0;
-  }
-
-  myMosaic->flags |= ID_IMAGE_MOSAIC_PHOTCAL;  // set this flag
-
-  if (testImage) {
-    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq);
   }
 
@@ -1400 +1372 @@
 }
   
-int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (int j = 0; j < N_onMosaic[Nmos]; j++) {
+      
+    off_t m = MosaicToMeasure[Nmos][j];
+    off_t c = MosaicToCatalog[Nmos][j];
+      
+    // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF
+    float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
+    // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
+    float Mgrid    = getMgrid (m, c);
+    float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
+    float Mflat    = getMflat (m, c, flatcorr, catalog);
+
+    off_t n = catalog[c].measureT[m].averef;
+    float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
+    
+    float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+
+    int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
+    
+    fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad);
+  }
+  return TRUE;
+}
+
+int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
   int i;
@@ -1408 +1408 @@
 
   fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
-
-  if (PoorImages) {
-    // XXX use bad stars and measurements for PoorImages? or not?
-    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
-    IMAGE_BAD = 0;
-  }
 
   off_t Nmax = 0;
@@ -1421 +1415 @@
 
   SetMmosInfo summary;
-  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
+  SetMmosInfoInit (&summary, Nmax, FALSE);
 
   pthread_attr_t attr;
@@ -1445 +1439 @@
     threadinfo[i].flatcorr = flatcorr;
 
-    // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop,
-    // PoorImages) used in the threads to allocate the arrays and set the MaxIterations
-    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
+    // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop) used in
+    // the threads to allocate the arrays and set the MaxIterations
+    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE);
     pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
   }
@@ -1492 +1486 @@
   npass_output ++;
 
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-    STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-  }
   return TRUE;
 }
@@ -1505 +1494 @@
 
   SetMmosInfo results;
-  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages);
+  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE);
 
   while (1) {
@@ -1892 +1881 @@
   int N = 0;
   for (int i = 0; i < Nmosaic; i++) {
-    // if (mosaic[i].flags & IMAGE_BAD) continue;
     if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) continue;
     if (mosaic[i].skipCal) continue;
@@ -1906 +1894 @@
   float MaxChiSq = stats.Upper90; 
 
-  // old version based on Mcal value
-  // MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
-  /// MedOffset = stats.median;
-
   liststats (slist, NULL, NULL, N, &stats);
   float MaxScatter = stats.Upper90;
-  // old: MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
 
   fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f\n", MaxChiSq, MaxScatter);
@@ -1927 +1910 @@
     int mark = FALSE;
 
-    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) {
+    if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) {
       mark = TRUE;
       Nfew ++;

trunk/Ohana/src/relphot/src/StarOps.c

@@ -662 +662 @@
 
 // this function only operations on the PSF magnitudes
+
+// XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers
 
 # define NSIGMA_CLIP 3.0

trunk/Ohana/src/relphot/src/TGroupOps.c

@@ -194 +194 @@
   for (off_t i = 0; i < Nsubset; i++) {
     ImageToTGroup[i] = NULL;
-
-    // ignore non-GPC1 images
-    // if (!isGPC1chip(subset[i].photcode)) {
-    //   Nsimple ++;
-    //   continue;
-    // }
 
     // ignore mosaic images (photcode == 0)
@@ -459 +453 @@
   // this tgroup has been identified as poor.
   tgroup->flags |= ID_IMAGE_NIGHT_POOR;
+  tgroup->McalPSF  = 0.0;
+  tgroup->McalAPER = 0.0;
+  tgroup->dKlam    = 0.0;
 
   // all associated mosaics should be marked as coming from a bad night
@@ -473 +470 @@
     image[im].flags |= ID_IMAGE_NIGHT_POOR;
   }
-
-  tgroup->McalPSF  = 0.0;
-  tgroup->McalAPER = 0.0;
-  tgroup->dKlam    = 0.0;
 }
 
@@ -552 +545 @@
   fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch);
 
-  // XXX print results for testing
-  for (int i = 0; FALSE && (i < NtgroupTimes); i++) {
-    TGroup *tgroup = tgroupTimes[i][0].byCode;
-    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (0) fprintf (stderr, "  TGROUP time %.0f photcode %d, Nmeasure: %d, Nimage: %d\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, (int) tgroup[j].Nmeasure, (int) tgroup[j].Nimage);
-    }
-  }
   return (TRUE);
 }
@@ -574 +560 @@
 
   TGroup *myGroup = ImageToTGroup[idx];
-  TGTimes *parent = (TGTimes *) myGroup->parent;
-  if (!myGroup) {
-    // Image *image = getimage(idx);
-    // fprintf (stderr, "unmatched image %s\n", image[0].name);
-    // skip measurements from simple chips (not matched to a tgroup by definition)
-    return;
-  }
-
-  // test to check we got the right match:
-  {
-    Image *image = getimage(idx);
-    // XXX we are now matching with just tzero.  be careful for cameras with drift
-    // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
-    unsigned int imageStart = image[0].tzero;
-    // XXX to do this check I need a pointer to the group master
-    if ((imageStart < parent->start) || (imageStart > parent->stop)) {
-      fprintf (stderr, "error in image to tgroup match\n");
-      abort();
-    }
-  }
+  if (!myGroup) return;
 
   // this measurement is on one of my tgroups, mark it as mine.
@@ -611 +578 @@
 }
 
-// XXX return the error as well as the value
+// dZpt (if supplied) returns the error on Mgrp
 float getMgrp (off_t meas, int cat, float airmass, float *dZpt) {
 
@@ -638 +605 @@
 }
 
-// XXX remove the lists and replace with StatDataSet structs?
-// XXX move the StatDataSetAlloc up?
 typedef struct {
   int Nfew;
@@ -650 +615 @@
   int Nskip;
   off_t Nmax;
-  int PoorImages;
   FitDataSet psfStars;
   FitDataSet kronStars;
@@ -669 +633 @@
 int   setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr);
 void *setMgrp_worker (void *data);
-int   setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
-
-void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages) {
+int   setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr);
+
+void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists) {
   info->Nfew = 0;
   info->Nbad = 0;
@@ -682 +646 @@
 
   info->Nmax = Nmax;
-  info->PoorImages = PoorImages;
 
   if (allocLists) {
@@ -744 +707 @@
 }
 
-int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMgrp (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
   off_t N;
@@ -755 +718 @@
   // so do not run setMgrp in threaded mode if PLOTSTUFF is set
   if (NTHREADS && !PLOTSTUFF) {
-    int status = setMgrp_threaded (catalog, PoorImages, flatcorr);
+    int status = setMgrp_threaded (catalog, flatcorr);
     return status;
   }
@@ -770 +733 @@
 
   SetMgrpInfo info;
-  SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages);
+  SetMgrpInfoInit (&info, Nmax, TRUE);
 
   for (off_t i = 0; i < NtgroupTimes; i++) {
@@ -776 +739 @@
     for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
       setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr);
-      // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE);
+      // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n",
+      // ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE);
     }
   }
@@ -799 +763 @@
 int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) {
 
-  off_t j, NimageReal;
-
-  // At some point, we identify good and bad tgroups.  We mark bad tgroups and do not fit
-  // a zero point for them. Instead we will fit those exposures or images separately
-  if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_GOOD_NIGHT) {
-    int bad = myTGroup->flags & ID_IMAGE_NIGHT_POOR;
-    if (bad) return TRUE;
-  }
+  off_t j;
+
+  // calculate the statistics for both good and bad nights, but only set Mcal for the good nights
+  // we only skip the statistics for nights with too few measurements or exposures
 
   FitDataSet *psfStars = &info->psfStars;
@@ -815 +775 @@
   assert (Ngrp < NtgroupTimes);
 
-  // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
-  getimages (&NimageReal, NULL);
-
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  int testImage = FALSE;
   TGTimes *tgroup = (TGTimes *) myTGroup->parent;
   double mjdStart = ohana_sec_to_mjd (tgroup->start);
@@ -827 +783 @@
   myTGroup->flags &= ~ID_IMAGE_TGROUP_PHOTCAL; 
   
-  // testImage = TRUE;
-  testImage |= (abs(mjdStart - 57353) < 0.1);
-  testImage |= (abs(mjdStart - 57952) < 0.1);
-  testImage |= (abs(mjdStart - 55843) < 0.1);
-  testImage |= (abs(mjdStart - 56586) < 0.1);
-  testImage |= (abs(mjdStart - 57323) < 0.1);
-  //testImage |= (abs(mjdStart - 57971) < 0.1);
-  //testImage |= (abs(mjdStart - 57974) < 0.1);
-  //testImage |= (abs(mjdStart - 57978) < 0.1);
-  //testImage |= (abs(mjdStart - 57981) < 0.1);
+  // in clean_tgroups (run on each iteration), we identify good and bad tgroups.  We mark
+  // bad tgroups and do not save a zero point for them. Instead we will fit those
+  // exposures or images separately
+  int useMgrp = TRUE;  
+  if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_GOOD_NIGHT) {
+    if (myTGroup->flags & ID_IMAGE_NIGHT_POOR) useMgrp = FALSE;
+  }
+
+  // for testing, supply the MJD of a night or nights here to dump the full list of measurements
+  int testImage = FALSE;
+  // testImage |= (abs(mjdStart - 57353) < 0.1);
 
   FILE *fout = NULL;
@@ -905 +862 @@
     float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON);
 
-    // XXX byCode : here we are only checking for equiv photcodes (not active photcodes)
+    // here we are only checking for equiv photcodes (not active photcodes)
     PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     if (!code) goto skip;
@@ -942 +899 @@
   }
 
-  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
+  /* skip nights with too few good measurements or too many bad measurements */
   int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure);
   if (mark) {
-    if (TRUE) {
-      TGTimes *mygroup = (TGTimes *) myTGroup->parent;
-      if (0) { // XXX EAM
-        fprintf (stderr, "marked tgroup %f,%d, (%d < %d) || (%d < %f*"OFF_T_FMT")\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); 
-      }
-    }
     myTGroup->flags |= ID_IMAGE_PHOTOM_FEW;
+    myTGroup->McalPSF    = 0.0;
+    myTGroup->McalAPER   = 0.0;
+    myTGroup->dMcal      = NAN;
+    myTGroup->stdev      = NAN;
+    myTGroup->dMmin      = NAN;
+    myTGroup->dMmax      = NAN;
+    myTGroup->McalChiSq  = NAN;
+    myTGroup->nFitPhotom = 0;
     myTGroup->nValPhotom = N;
+    myTGroup->dKlam      = (TGROUP_FIT_AIRMASS) ? psfStars->bSaveArray[1][0] : 0;
     info->Nfew ++;
     if (testImage) {
@@ -986 +946 @@
 
   fit1d_irls (psfStars, N);
-
-  // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
-
-  // for now, I have no reason to measure these separately for camera-level images
-  myTGroup->McalPSF    = psfStars->bSaveArray[0][0];
+  fit1d_irls (kronStars, N);
+  fit1d_irls (brightStars, Nbright);
+
+  if (useMgrp) {
+    myTGroup->McalPSF    = psfStars->bSaveArray[0][0];
+    myTGroup->McalAPER   = kronStars->bSaveArray[0][0];   // XXX this does not make sense: I need to apply the airmass slope calculated above first
+    myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL;  // set this flag (unset by default)
+  } else {
+    myTGroup->McalPSF    = 0.0;
+    myTGroup->McalAPER   = 0.0;
+  }
+
+  // record statistics on the fit regardless if it is used or not
   myTGroup->dMcal      = psfStars->bSigma[0];
   myTGroup->stdev      = psfStars->sigma;
@@ -1000 +968 @@
   myTGroup->dKlam      = (TGROUP_FIT_AIRMASS) ? psfStars->bSaveArray[1][0] : 0;
 
-  // XXX this does not make sense: I need to apply the airmass slope calculated above first
-  fit1d_irls (kronStars, N);
-  myTGroup->McalAPER   = kronStars->bSaveArray[0][0];
+  // bright end scatter (systematic error)
+  myTGroup->dMsys = brightStars->sigma;
 
   if (testImage) {
     TGTimes *parent = (TGTimes *) myTGroup->parent;
     fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom);
+    fprintf (stderr, "%f %f  :  %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq);
   }
   if (PLOTSTUFF) {
     fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, N, ohana_sec_to_mjd(tgroup->start));
     plot_setMcal (psfStars->alldata->yVector, N);
-  }
-
-  // bright end scatter
-  // XXX this does not make sense: I need to apply the airmass slope calculated above first
-  // redo this by calculating the corrected bright stars mags 
-  fit1d_irls (brightStars, Nbright);
-  myTGroup->dMsys = brightStars->sigma;
-
-  // keep this??
-  if (myTGroup->McalPSF < -CLOUD_TOLERANCE) {
-    myTGroup->McalPSF = 0.0;
-  }
-
-  myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL;  // set this flag
-
-  if (testImage) {
-    fprintf (stderr, "%f %f  :  %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq);
   }
 
@@ -1070 +1021 @@
 }
 
-int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
   int i;
@@ -1077 +1028 @@
   Image *image = getimages (&N, NULL);
 
-  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
-
   off_t Nmax = 0;
   for (off_t i = 0; i < NtgroupTimes; i++) {
@@ -1088 +1037 @@
 
   SetMgrpInfo summary;
-  SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages);
+  SetMgrpInfoInit (&summary, Nmax, FALSE);
 
   pthread_attr_t attr;
@@ -1112 +1061 @@
     threadinfo[i].flatcorr = flatcorr;
 
-    // we do NOT allocate the arrays here, we only supply basic info (Nmax,
-    // PoorImages) used in the threads to allocate the arrays and set the MaxIterations
-    SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
+    // we do NOT allocate the arrays here, we only supply basic info used in the threads
+    // to allocate the arrays and set the MaxIterations
+    SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE);
     pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]);
   }
@@ -1169 +1118 @@
 
   SetMgrpInfo results;
-  SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
+  SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE); // allocate list, dlist arrays here
 
   while (1) {
@@ -1194 +1143 @@
 }
 
-// XXX should I split these stats out by photcode?
-// XXX double check the IMAGE_BAD filter:
 StatType statsTGroupM (Catalog *catalog) {
   OHANA_UNUSED_PARAM(catalog);
@@ -1209 +1156 @@
   ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
 
-  // is ID_IMAGE_PHOTOM_POOR used for tgroups?
-  int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR;
+  int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW;
 
   int n = 0;
@@ -1316 +1262 @@
   ALLOCATE (slist, double, NtgroupTimes*Nphotcodes);
 
+  // measure the good/bad statistics using nights which are actually calibrated
+  // but they will be applied to all nights below (allowing bad nights to become good)
   int N = 0;
   for (int i = 0; i < NtgroupTimes; i++) {
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & IMAGE_BAD) continue;
+      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
       mlist[N] = tgroup[j].McalChiSq;
       slist[N] = tgroup[j].stdev; // stdev of all measurements
-   // slist[N] = tgroup[j].dMsys; // stdev of bright measurements
       N++;
     }
@@ -1491 +1438 @@
 }
 
+// temporary test output 
 void dump_tgroups (Catalog *catalog, int Npass) {
 
@@ -1606 +1554 @@
   fclose (fout);
 }
+
+/*
+  for testing, put lines like these in relphot_images.c within the iteration loop
+      dump_tgroups (catalog, i + 40);
+      dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog
+      dump_tgroup_imstats (i + 40);
+*/

trunk/Ohana/src/relphot/src/relphot_images.c

@@ -41 +41 @@
   load_images (&db, skylist, &UserPatch, FALSE, USE_ALL_IMAGES);
   MARKTIME("-- load images: %f sec\n", dtime);
-
-  // checkImages("load images");
 
   /* unlock, if we can (else, unlocked below) */
@@ -63 +61 @@
     /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
 
-    // checkImages("load catalog");
-
     /* match measurements with images, mosaics */
     initImageBins  (catalog, Ncatalog, TRUE);
@@ -87 +83 @@
 
     setExclusions (catalog, Ncatalog, TRUE);
-    // checkImages("set Exclu");
 
     global_stats (catalog, Ncatalog, flatcorr, 0);
@@ -118 +113 @@
       SetZptIteration (i);
 
-      dump_tgroups (catalog, i);
-      dump_catalog (catalog, 0, i); // for a test, just dump a specific catalog
-      dump_tgroup_imstats (i);
-
       setMrel  (catalog, Ncatalog, flatcorr); // threaded
-      setMcalTest (catalog, flatcorr);
-
-      dump_tgroups (catalog, i + 20);
-      dump_catalog (catalog, 0, i + 20); // for a test, just dump a specific catalog
-      dump_tgroup_imstats (i + 20);
-
-      setMcal  (catalog, FALSE, flatcorr);
-      setMmos  (catalog, FALSE, flatcorr);
-      setMgrp  (catalog, FALSE, flatcorr);
-      dump_tgroups (catalog, i + 40);
-      dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog
-      dump_tgroup_imstats (i + 40);
+
+      setMcal  (catalog, flatcorr);
+      setMmos  (catalog, flatcorr);
+      setMgrp  (catalog, flatcorr);
 
       setMgrid (catalog, flatcorr);
@@ -166 +149 @@
     if (USE_GRID) dump_grid ();
 
-    /* set Mcal & Mmos for bad images */
-    int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
-    onlyPoorImages = FALSE; // XXX deactivate for now
-    setMcal  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
-    setMmos  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
-    setMgrp  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
@@ -177 +154 @@
     setMcalFromTGroups (); // copy per-tgroup calibrations to the images
 
-    if (SAVE_IMAGE_UPDATES) {
-
-      FITS_DB dbX;
-      off_t *LineNumber;
-      Image *image, *subset;
-      off_t Nimage, Nsubset, i, Nx;
-      char filename[1024];
-
-      gfits_db_init (&dbX);
-      dbX.lockstate = LCK_XCLD;
-      dbX.timeout   = 60.0;
-      dbX.mode      = db.mode;
-      dbX.format    = db.format;
-
-      snprintf (filename, 1024, "%s.bck", ImageCat);
-      if (!gfits_db_lock (&dbX, filename)) {
-        fprintf (stderr, "can't lock backup image image catalog\n");
-        return (FALSE);
-      }
-    
-      // apply the changes from the image subset to the full table:
-
-      // convert database table to internal structure (binary to Image)
-      // 'image' points to the same memory as db->ftable->buffer
-      image = gfits_table_get_Image (&db.ftable, &Nimage, &db.scaledValue, &db.nativeOrder);
-      if (!image) {
-        fprintf (stderr, "ERROR: failed to read images\n");
-        exit (2);
-      }
-      gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
-
-      subset = getimages (&Nsubset, &LineNumber);
-      for (i = 0; i < Nsubset; i++) {
-        if (LineNumber[i] == -1) continue;
-        memcpy (&image[LineNumber[i]], &subset[i], Nx);
-      }
-
-      // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
-      // the entire Image table, but with the updates in place I think this says: re-work the
-      // ftable/vtable usage in this program to be more sensible...
-      gfits_copy_header (&db.header,  &dbX.header);
-      gfits_copy_matrix (&db.matrix,  &dbX.matrix);
-      gfits_copy_header (&db.theader, &dbX.theader);
-      gfits_copy_ftable (&db.ftable,  &dbX.ftable);
-
-      dbX.ftable.header = &dbX.theader;
-      dbX.virtual = FALSE;
-
-      // save Images.dat using the copied structure
-      if (UPDATE_CATFORMAT) {
-        // ensure the db format is updated
-        dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
-        gfits_modify (&dbX.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT);
-
-        char photcodeFile[1024];
-        sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
-        SavePhotcodesFITS (photcodeFile);
-      }
-      dvo_image_save (&dbX, VERBOSE);
-      dvo_image_unlock (&dbX); 
-      gfits_db_free (&dbX);
-      MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
-    }
+    // calculate and save per-chip residuals here:
+    MagResidSave ("mag.resid.fits", catalog);
+
+    save_images_updates (&db);
 
     /* at this point, we have correct cal coeffs in the image/mosaic structures */

trunk/Ohana/src/relphot/src/relphot_parallel_images.c

@@ -112 +112 @@
 
     // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements
-    setMcal  (catalog, FALSE, flatcorr);
-    setMmos  (catalog, FALSE, flatcorr);
+    setMcal  (catalog, flatcorr);
+    setMmos  (catalog, flatcorr);
     // setMgrp  (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups
     MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
@@ -144 +144 @@
   client_logger_message ("all hosts are done the loops\n");
 
-  /* set Mcal & Mmos for bad images (for stack_and_warps, force images to be measured) */
-  int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
-  setMcal  (catalog, onlyPoorImages, flatcorr);
-  setMmos  (catalog, onlyPoorImages, flatcorr);
-  // setMgrp  (catalog, onlyPoorImages, flatcorr); XXX see above
   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 

trunk/Ohana/src/relphot/src/relphot_parallel_regions.c

@@ -1 +1 @@
 # include "relphot.h"
-int save_images_backup (FITS_DB *db);
 
 int relphot_parallel_regions (SkyTable *sky) {
@@ -101 +100 @@
   return TRUE;
 }
-
-int save_images_backup (FITS_DB *db) {
-
-  FITS_DB dbX;
-  char filename[1024];
-
-  INITTIME;
-
-  gfits_db_init (&dbX);
-  dbX.lockstate = LCK_XCLD;
-  dbX.timeout   = 60.0;
-  dbX.mode      = db->mode;
-  dbX.format    = db->format;
-
-  snprintf (filename, 1024, "%s.bck", ImageCat);
-  if (!gfits_db_lock (&dbX, filename)) {
-    fprintf (stderr, "can't lock backup image image catalog\n");
-    return (FALSE);
-  }
-    
-  // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
-  // the entire Image table, but with the updates in place I think this says: re-work the
-  // ftable/vtable usage in this program to be more sensible...
-  gfits_copy_header (&db->header,  &dbX.header);
-  gfits_copy_matrix (&db->matrix,  &dbX.matrix);
-  gfits_copy_header (&db->theader, &dbX.theader);
-  gfits_copy_ftable (&db->ftable,  &dbX.ftable);
-
-  dbX.ftable.header = &dbX.theader;
-
-  dvo_image_save (&dbX, VERBOSE);
-  dvo_image_unlock (&dbX); 
-  MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
-
-  gfits_db_free (&dbX);
-  return TRUE;
-}