Changeset 41555

Timestamp:

Apr 24, 2021, 2:28:09 PM (5 years ago)

Author:

eugene

Message:

add function to measure the residual stats for each chip, fix resets and flags on resets; fix update of image mcal from tgroup and mosaic; improved determination of bad nights and bad mosaics; allow tgroups to fit only zero points without airmass terms; clean up the test detection dumps; soften the errors to avoid overweighting extreme outliers with high S/N; more careful sequencing of fitting modes; add manual iteration option

Location:

trunk/Ohana/src/relphot

Files:

• include/relphot.h (modified) (6 diffs)
• src/ImageOps.c (modified) (9 diffs)
• src/MosaicOps.c (modified) (12 diffs)
• src/TGroupOps.c (modified) (21 diffs)
• src/ZeroPointModes.c (modified) (9 diffs)
• src/args.c (modified) (2 diffs)
• src/liststats.c (modified) (6 diffs)
• src/load_images.c (modified) (2 diffs)
• src/relphot_images.c (modified) (3 diffs)
• src/setMrelCatalog.c (modified) (4 diffs)

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

@@ -88 +88 @@
   float dMcal;
   float dMsys;
+  float stdev;
+  float dMmin;
+  float dMmax; 
   float McalChiSq;
   float secz;
@@ -107 +110 @@
   float dMsys;
   float stdev;
+  float dMmin;
+  float dMmax; 
   float McalChiSq;
   unsigned int nFitPhotom;
@@ -155 +160 @@
   double min;
   double max;
+  double Upper90;
   double Upper80;
   double Lower20;
+  double Lower10;
   double total;
   int    Nmeas;
@@ -415 +422 @@
 int    TEST_IMAGE2;
 
+int    MANUAL_ITERATION;
 int    NLOOP;
 int    NGRID;
@@ -643 +651 @@
 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
 
+void          setMcalTest         PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+
 int           setMave             PROTO((Catalog *catalog, int Ncatalog));
 void          set_ZP              PROTO((double ZERO));
@@ -835 +845 @@
 void dump_tgroups (Catalog *catalog, int Npass);
 void dump_catalog (Catalog *catalog, off_t c, int Npass);
+void dump_tgroup_imstats (int Npass);
 
 void SetZptIteration (int current);

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

@@ -1 +1 @@
 # include "relphot.h"
-# define BASIC_STATS 0
 
 // the MeasureToImage, ImageToCatalog, and ImageToMeasure arrays are a substantial part of the memory footprint.
@@ -631 +630 @@
     int isMosaic  = isMosaicChip(image[i].photcode);
 
+    // unset at start by default (set if actually used below)
+    image[i].flags &= ~ID_IMAGE_IMAGE_PHOTCAL;
+
     // if requested, freeze mosaic chips:
     if (FREEZE_IMAGES && isMosaic) continue;
@@ -719 +721 @@
       }
 
+      float Moff = Mmos + Mgrp + Mgrid - Mflat;
+
       PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
       if (!code) goto skip;
@@ -727 +731 @@
 
     skip:
-      psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat;
+      psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff;
       psfStars.alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
 
@@ -734 +738 @@
 
       if (isfinite(MrelKron) && isfinite(MsysKron)) {
-        kronStars.alldata-> yVector[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat;
+        kronStars.alldata-> yVector[Nkron] = MsysKron - MrelKron - Moff;
         kronStars.alldata->dyVector[Nkron] = psfStars.alldata->dyVector[Nref];
         Nkron ++;
@@ -765 +769 @@
       }      
     }
+    if (mark) continue;
 
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
@@ -776 +781 @@
     // 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?
-    fit1d_irls (&kronStars, Nkron);
-    image[i].McalAPER   = kronStars.bSaveArray[0][0];
+    // XXX EAM : temporarily use this to track calibration quality
+    if (0) {
+      fit1d_irls (&kronStars, Nkron);
+      image[i].McalAPER   = kronStars.bSaveArray[0][0];
+    }
 
     if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
@@ -795 +803 @@
       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
@@ -814 +824 @@
   FitDataSetFree (&psfStars);
 
+  return;
+}
+
+/* determine McalTEST values for all images -- this is not used to set measure.Mcal, but only to test */
+void setMcalTest (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+
+  off_t i, j, m, c, n;
+  int Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  off_t Nmax = 0;
+  for (i = 0; i < Nimage; i++) {
+    Nmax = MAX (Nmax, N_onImage[i]);
+  }
+
+  // we are making a 0-order fit and not doing bootstrap analysis:
+  FitDataSet psfStars;
+  FitDataSetAlloc (&psfStars, Nmax, 0, 0);
+
+  // for testing, need to allow this to be optional
+  if (FALSE && UseStandardOLS(ZPT_IMAGES)) {
+    psfStars.MaxIterations = 0;
+  }
+  psfStars.MaxIterations = 0;
+
+  Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0;
+
+  int Ncalibrated = 0;
+  for (i = 0; i < Nimage; i++) {
+    if (image[i].photcode == 0) continue; // skip the PHU images
+
+    off_t Nref = 0;    // number of stars used to measure McalPSF
+
+    if (N_onImage[i] == 0) {
+      // fprintf (stderr, "image missing detections? %s %d "OFF_T_FMT"\n", image[i].name, image[i].nstar, N_onImage[i]);
+      continue;
+    }
+
+    for (j = 0; j < N_onImage[i]; j++) {
+      
+      m = ImageToMeasure[i][j];
+      c = ImageToCatalog[i][j];
+      
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
+          Nbad++;
+          continue;
+      }
+      float Mmos  = getMmos  (m, c);
+      if (isnan(Mmos)) {
+          Nmos ++;
+          continue;
+      }
+      float Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL); // ignore error for now?
+      if (isnan(Mgrp)) {
+          Ngrp ++;
+          continue;
+      }
+      float Mgrid = getMgrid (m, c);
+      if (isnan(Mgrid)) {
+          Ngrid++;
+          continue;
+      }
+
+      // MrelPSF is the average magnitude for this star.  
+      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
+      if (isnan(MrelPSF)) {
+          Nrel ++;
+          continue;
+      }
+      
+      // image.Mcal is not supposed to include the flat-field correction, so we need to
+      // apply that offset as well here for this image (in other words, each detection is
+      // being compared to the model, excluding the zero point, Mcal.  The model includes
+      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
+
+      float Mflat = getMflat (m, c, flatcorr, catalog);
+
+      // get the PSF magnitude for thie measurement, with airmass slope applied
+      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);
+      if (isnan(MsysPSF)) {
+        Nsys++;
+        continue;
+      }
+
+      float Moff = Mmos + Mgrp + Mgrid - Mflat;
+
+      psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff;
+      psfStars.alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+      Nref++;
+    }
+
+    if (Nref < 5) {
+      // fprintf (stderr, "very few detections? %s %d "OFF_T_FMT"\n", image[i].name, image[i].nstar, Nref);
+      // continue;
+      image[i].McalAPER   = NAN;
+      image[i].dMcal      = NAN;
+      image[i].nFitPhotom = Nref;
+      image[i].McalChiSq  = NAN;
+      continue;
+    }
+
+    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
+    fit1d_irls (&psfStars, Nref);
+    image[i].McalAPER   = psfStars.bSaveArray[0][0];
+    image[i].dMcal      = psfStars.bSigma[0];
+    image[i].nFitPhotom = psfStars.Nmeas;
+    image[i].McalChiSq  = psfStars.chisq;
+    Ncalibrated ++;
+  }
+
+  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);
+
+  FitDataSetFree (&psfStars);
   return;
 }

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

@@ -685 +685 @@
       // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
 
+      if (!isfinite(image[m].McalPSF)) {
+        image[m].McalPSF = 0.0;
+        image[m].McalAPER = 0.0;
+        image[m].dMcal = 0.0;
+        image[m].McalChiSq = 0.0;
+        fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
+      }
+
       McalPSF   += image[m].McalPSF;
       McalAPER  += image[m].McalAPER;
@@ -694 +702 @@
       image[m].dMcal     = NAN;
       image[m].McalChiSq = NAN;
+
     }
     dS /= MosaicN_Image[i];
@@ -734 +743 @@
     for (j = 0; j < MosaicN_Image[i]; j++) {
       im = MosaicToImage[i][j];
-      image[im].McalPSF    += mosaic[i].McalPSF;
-      image[im].McalAPER   += mosaic[i].McalAPER;
-      image[im].dMcal       = mosaic[i].dMcal;
-      image[im].McalChiSq   = mosaic[i].McalChiSq;
-      image[im].ubercalDist = mosaic[i].ubercalDist;
-      image[im].dMagSys     = mosaic[i].dMsys;
-      image[im].nFitPhotom  = mosaic[i].nFitPhotom;
+      if (mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL) {
+        image[im].McalPSF     = mosaic[i].McalPSF;
+        image[im].McalAPER    = mosaic[i].McalAPER;
+        image[im].dMcal       = mosaic[i].dMcal;
+        image[im].McalChiSq   = mosaic[i].McalChiSq;
+        image[im].ubercalDist = mosaic[i].ubercalDist;
+        image[im].dMagSys     = mosaic[i].dMsys;
+        image[im].nFitPhotom  = mosaic[i].nFitPhotom;
+      }
       image[im].flags      |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW);
       image[im].flags      |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
       image[im].flags      |= (mosaic[i].flags & ID_IMAGE_MOSAIC_POOR);
+      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_TGROUP_PHOTCAL);
+      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL);
+      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_IMAGE_PHOTCAL);
     }
     mosaic[i].McalPSF  = 0.0;
@@ -756 +770 @@
 
   // this mosaic has been identified as poor.
-  mosaic[myMosaic].flags |= ID_IMAGE_PHOTOM_POOR;
+  mosaic[myMosaic].flags |= ID_IMAGE_MOSAIC_POOR;
   mosaic[myMosaic].McalPSF  = 0.0;
   mosaic[myMosaic].McalAPER = 0.0;
@@ -776 +790 @@
 
   // this mosaic has been identified as good.
-  mosaic[myMosaic].flags &= ~ID_IMAGE_PHOTOM_POOR;
+  mosaic[myMosaic].flags &= ~ID_IMAGE_MOSAIC_POOR;
 
   image = getimages (&Nimage, NULL);
@@ -1297 +1311 @@
     assert (Nbright >= 0);
 
-    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    float Moff =  Mcal + Mgrp + Mgrid - Mflat;
+
+    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff;
     psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
 
-    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
+    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Moff;
     kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N];
 
@@ -1337 +1353 @@
   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;         
@@ -1863 +1880 @@
 void clean_mosaics () {
 
-  off_t i, N, mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq;
-  double *mlist, *slist, *dlist;
-  double MaxOffset, MedOffset, MaxScatter;
+  double *mlist, *slist;
 
   if (!MOSAIC_ZEROPT) return;
@@ -1874 +1889 @@
   ALLOCATE (mlist, double, Nmosaic);
   ALLOCATE (slist, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  for (i = N = 0; i < Nmosaic; i++) {
+
+  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;
-    mlist[N] = mosaic[i].McalPSF;
-    slist[N] = mosaic[i].dMcal;
-    dlist[N] = 1;
+    mlist[N] = mosaic[i].McalChiSq;
+    slist[N] = mosaic[i].stdev;
     N++;
   }
@@ -1889 +1903 @@
   liststats_setmode (&stats, "MEAN");
 
-  liststats (mlist, dlist, NULL, N, &stats);
-  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
-  MedOffset = stats.median;
-
-  liststats (slist, dlist, NULL, N, &stats);
-  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
-  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
+  liststats (mlist, NULL, NULL, N, &stats);
+  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);
   
-  // XXX this needs to be configurable
-  float MaxChiSq = 10.0;
-
-  Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0;
-  for (i = 0; i < Nmosaic; i++) {
+  int Ntotal = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0;
+  for (int i = 0; i < Nmosaic; i++) {
     // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1907 +1923 @@
     Ntotal ++;
 
+    if (mosaic[i].skipCal) continue;
+
+    int mark = FALSE;
+
     if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) {
-      continue;
-    }
-    if (mosaic[i].skipCal) continue;
-
-    mark = FALSE;
-    if (mosaic[i].dMcal > MaxScatter) {
+      mark = TRUE;
+      Nfew ++;
+    }
+    if (mosaic[i].stdev > MaxScatter) {
       mark = TRUE;
       Nscatter ++;
-    }
-    if (fabs(mosaic[i].McalPSF - MedOffset) > MaxOffset) {
-      mark = TRUE;
-      Noffset ++;
-    }
-    if (mosaic[i].nFitPhotom < 0.7*mosaic[i].nValPhotom) {
-      mark = TRUE;
-      Nrejects ++;
     }
     if (mosaic[i].McalChiSq > MaxChiSq) {
@@ -1937 +1947 @@
   }
 
-  fprintf (stderr, OFF_T_FMT" of "OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Ntotal, Nscatter, Noffset);
+  fprintf (stderr, "%d of %d mosaics marked poor (%d scatter, %d chisq, %d few)\n",  Nmark, Ntotal, Nscatter, Nchisq, Nfew);
 
   free (mlist);
   free (slist);
-  free (dlist);
 }
 

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

@@ -2 +2 @@
 
 void assignMosaicsToTGroups (void);
+int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr);
 
 // tgroupTimes carries the times of the tgroups (initially, just the photometric nights)
@@ -82 +83 @@
       tgroup[j].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
       tgroup[j].dMsys     = 0.0;
+      tgroup[j].dMmin     = NAN;
+      tgroup[j].dMmax     = NAN;
       tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
 
@@ -142 +145 @@
     tgroup[j].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
     tgroup[j].dMsys     = 0.0;
+    tgroup[j].dMmin     = NAN;
+    tgroup[j].dMmax     = NAN;
     tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
 
@@ -421 +426 @@
         off_t im = tgroup[i].image[j];
         double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF;
-        image[im].McalPSF    += Mgrp;
-        image[im].McalAPER   += Mgrp;
-        image[im].dMcal       = tgroup[i].dMcal;
-        image[im].McalChiSq   = tgroup[i].McalChiSq;
-        image[im].dMagSys     = tgroup[i].dMsys;
-        image[im].nFitPhotom  = tgroup[i].nFitPhotom;
+        if (tgroup[i].flags & ID_IMAGE_TGROUP_PHOTCAL) {
+          image[im].McalPSF     = Mgrp;
+          image[im].McalAPER    = Mgrp;
+          image[im].dMcal       = tgroup[i].dMcal;
+          image[im].McalChiSq   = tgroup[i].McalChiSq;
+          image[im].dMagSys     = tgroup[i].dMsys;
+          image[im].nFitPhotom  = tgroup[i].nFitPhotom;
+        }
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // probably not necessary : it is set in markBadTGroup
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // probably not necessary : it is set in markBadTGroup
@@ -432 +439 @@
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL);
+        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_TGROUP_PHOTCAL);
+        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_MOSAIC_PHOTCAL);
+        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_IMAGE_PHOTCAL);
         
         // fprintf (stderr, "TG to IMAGE: %f +/- %f -> %f (%d)\n", tgroup[i].McalPSF, tgroup[i].dMcal, image[im].McalPSF, tgroup[i].nFitPhotom);
@@ -687 +697 @@
       info->psfStars.MaxIterations = 0;
     }
-    if (1) {
-      FitDataSetAddPriors (&info->psfStars);
+    FitDataSetAddPriors (&info->psfStars); // why only psfStars?
+    if (TGROUP_FIT_AIRMASS) {
       info->psfStars.bPriorValue[1] = 0.0;  // note that we are fitting relative to the nominal slope
       info->psfStars.bPriorSigma[1] = 0.04; // XXX this prior sigma needs to be user-configured
@@ -812 +822 @@
   int testImage = FALSE;
   TGTimes *tgroup = (TGTimes *) myTGroup->parent;
-  // testImage |= (abs(tgroup->start - 1245283200) < 10);
+  double mjdStart = ohana_sec_to_mjd (tgroup->start);
+
+  // unset this flag at start (set below if zeropoint is calculated)
+  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);
 
   FILE *fout = NULL;
@@ -820 +843 @@
     snprintf (filename, 64, "test.%05d.%02d.dat", (int) Ngrp, npass_output);
     fout = fopen (filename, "w");
+    fprintf (fout, "# tgroup %s (%.0f)\n", ohana_sec_to_date(tgroup->start), mjdStart);
+    save_test_night_measures (fout, myTGroup, catalog, flatcorr);
+    fclose (fout);
   }
 
@@ -831 +857 @@
     off_t c = myTGroup->catalog[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);     // image zero point
-      float Mmos     = getMmos  (m, c);                    // mosaic zero point
-      // XXX these two should both be 0.0 if we are fitting tgroups : skip those calls above?
-
-      float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
-      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-      float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
-
-      off_t n = catalog[c].measureT[m].averef;
-
-      // magnitude for this measurement
-      float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-
-      // for historical reasons, Mflat is defined with the wrong sign
-      float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat;
-
-      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
-
-      fprintf (fout, "%f %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, catalog[c].measureT[m].airmass, delta, isBad);
-    }
-
     if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
       info->Nbad ++;
@@ -916 +918 @@
     assert (Nbright >= 0);
 
-    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
-    psfStars->alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0);
+    float Moff =  Mcal + Mmos + Mgrid - Mflat;
+
+    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff;
     psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
 
-    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
-    kronStars->alldata-> xVector[N] = psfStars->alldata-> xVector[N];
+    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Moff;
     kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N];
+
+    if (TGROUP_FIT_AIRMASS) {
+      psfStars->alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0);
+      kronStars->alldata-> xVector[N] = psfStars->alldata-> xVector[N];
+    }
 
     if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
       brightStars->alldata-> yVector[Nbright] = psfStars->alldata-> yVector[N];
-      brightStars->alldata-> xVector[Nbright] = psfStars->alldata-> xVector[N];
       brightStars->alldata->dyVector[Nbright] = psfStars->alldata->dyVector[N];
+      if (TGROUP_FIT_AIRMASS) {
+        brightStars->alldata-> xVector[Nbright] = psfStars->alldata-> xVector[N];
+      }
       Nbright ++;
     }
     N++;
   }
-  if (fout) { fclose (fout); }
 
   /* too few good measurements or too many bad measurements (skip in PoorImages run) */
   int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure);
   if (mark) {
-    if (VERBOSE2) {
+    if (TRUE) {
       TGTimes *mygroup = (TGTimes *) myTGroup->parent;
-      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);
+      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->nValPhotom = N;
     info->Nfew ++;
     if (testImage) {
@@ -949 +960 @@
   } else {
     myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW;
+  }
+
+  /* we have a non-trivial fraction of measurements which are extreme outliers, and
+     measurements with very small errors.  To avoid over-weighting measurements with very
+     small errors, let's add something in quadarature.  
+  */
+
+  // use liststats to find the 20-pct, median, 80-pct points
+  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)
+
+  TGTimes *mygroup = (TGTimes *) myTGroup->parent;
+  fprintf (stderr, "TGroup Stats %f,%d : %d %d %6.3f %6.3f\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, (int) N, (int) myTGroup->Nmeasure, stats.median, altSigma);
+
+  // 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;
   }
 
@@ -962 +993 @@
   myTGroup->dMcal      = psfStars->bSigma[0];
   myTGroup->stdev      = psfStars->sigma;
+  myTGroup->dMmin      = psfStars->min;
+  myTGroup->dMmax      = psfStars->max;
   myTGroup->McalChiSq  = psfStars->chisq;
   myTGroup->nFitPhotom = psfStars->Nmeas;
   myTGroup->nValPhotom = N;
-  myTGroup->dKlam      = psfStars->bSaveArray[1][0];
+  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
@@ -991 +1024 @@
   }
 
+  myTGroup->flags |= ID_IMAGE_TGROUP_PHOTCAL;  // set this flag
+
   if (testImage) {
     fprintf (stderr, "%f %f  :  %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq);
@@ -998 +1033 @@
 }
   
+int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (int j = 0; j < myTGroup->Nmeasure; j++) {
+      
+    off_t m = myTGroup->measure[j];
+    off_t c = myTGroup->catalog[j];
+      
+    // XXX these two should both be 0.0 if we are fitting tgroups
+    float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);     // image zero point
+    float Mmos     = getMmos  (m, c);                    // mosaic zero point
+
+    float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+    float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
+    float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+
+    off_t n = catalog[c].measureT[m].averef;
+
+    // magnitude for this measurement
+    float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
+
+    // for historical reasons, Mflat is defined with the wrong sign
+    float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat - myTGroup->McalPSF;
+
+    int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
+
+    double mjdMeas = ohana_sec_to_mjd (catalog[c].measureT[m].t);
+
+    fprintf (fout, "%f %f : %f %f %f %f %f  : %f %f %d : %f\n",
+             catalog[c].averageT[n].R, catalog[c].averageT[n].D,
+             MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, catalog[c].measureT[m].airmass, delta, isBad, mjdMeas);
+
+  }
+  return TRUE;
+}
+
 int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
@@ -1221 +1293 @@
 }
 
-/* mark tgroup if: abs(Mcal - <Mcal>) too large, dMcal too large */
+/* new rules for bad tgroups (bad nights) 
+ * too few exposures (Nmosaic < X)
+ * dMsys or stdev > X
+ * X%-ile of dMsys or stdev
+ * McalChiSq > X
+ * X%-ile of McalChiSq
+
+ * I can define rules for these but they would be arbitrary
+ * I think I should do an analysis of the observed distribution
+
+ */
+
 void clean_tgroups () {
 
-  off_t mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq;
-  double *mlist, *slist, *dlist;
-  double MaxOffset, MedOffset, MaxScatter;
+  double *mlist, *slist;
 
   if (!TGROUP_ZEROPT) return;
@@ -1234 +1315 @@
   ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes);
   ALLOCATE (slist, double, NtgroupTimes*Nphotcodes);
-  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
 
   int N = 0;
@@ -1241 +1321 @@
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
       if (tgroup[j].flags & IMAGE_BAD) continue;
-      mlist[N] = tgroup[j].McalPSF;
-      slist[N] = tgroup[j].dMcal;
-      dlist[N] = 1;
+      mlist[N] = tgroup[j].McalChiSq;
+      slist[N] = tgroup[j].stdev; // stdev of all measurements
+   // slist[N] = tgroup[j].dMsys; // stdev of bright measurements
       N++;
     }
@@ -1251 +1331 @@
   liststats_setmode (&stats, "MEAN");
 
-  liststats (mlist, dlist, NULL, N, &stats);
-  MaxOffset = MAX (NIGHT_OFFSET, 3*stats.sigma);
-  MedOffset = stats.median;
-
-  liststats (slist, dlist, NULL, N, &stats);
-  MaxScatter = MAX (NIGHT_SCATTER, 2*stats.median);
-  fprintf (stderr, "TGROUPS: median Mrel: %f, median dMrel: %f, Max Offset: %f, Max Scatter: %f\n", MedOffset, stats.median, MaxOffset, MaxScatter);
+  liststats (mlist, NULL, NULL, N, &stats);
+  float MaxChiSq = stats.Upper90; 
+
+  liststats (slist, NULL, NULL, N, &stats);
+  float MaxScatter = stats.Upper90;
+
+  fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f\n", MaxChiSq, MaxScatter);
   
-  // XXX this needs to be configurable
-  float MaxChiSq = 10.0;
-
-  Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0;
+  int Ntotal = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0;
   for (int i = 0; i < NtgroupTimes; i++) {
     TGroup *tgroup = tgroupTimes[i][0].byCode;
@@ -1269 +1346 @@
       Ntotal ++;
 
-      // ignore the nights with too few measurements
-      // XXX what to do with nights with too few measurements?
-      // XXX probably free the images to be fitted
+      int mark = FALSE;
+
+      // too few measurements (set in setMgrp_tgroup)
       if (tgroup[j].flags & ID_IMAGE_PHOTOM_FEW) {
-        markBadTGroup(&tgroup[j]); // mark the image associated with a bad night
-        continue;
+        mark = TRUE;
+        NfewNights ++;
       }
-      
-      mark = FALSE;
-      if (tgroup[j].dMcal > MaxScatter) {
+      // too few exposures (configure)
+      if (tgroup[j].Nmosaic < 4) {
+        mark = TRUE;
+        NfewExp ++;
+      }
+      // scatter too large
+      if (tgroup[j].stdev > MaxScatter) {
         mark = TRUE;
         Nscatter ++;
       }
-      if (fabs(tgroup[j].McalPSF - MedOffset) > MaxOffset) {
-        mark = TRUE;
-        Noffset ++;
-      }
-      if (tgroup[j].nFitPhotom < 0.7*tgroup[j].nValPhotom) {
-        mark = TRUE;
-        Nrejects ++;
-      }
+      // chisq too large
       if (tgroup[j].McalChiSq > MaxChiSq) {
         mark = TRUE;
@@ -1303 +1377 @@
   }
 
-  fprintf (stderr, OFF_T_FMT" of "OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset, "OFF_T_FMT" rejects, "OFF_T_FMT" chisq)\n",  Nmark, Ntotal, Nscatter, Noffset, Nrejects, Nchisq);
+  fprintf (stderr, "%d of %d tgroups marked poor (%d scatter, %d few nights, %d few exposures, %d chisq)\n",  Nmark, Ntotal, Nscatter, NfewNights, NfewExp, Nchisq);
 
   free (mlist);
   free (slist);
-  free (dlist);
 }
 
@@ -1465 +1538 @@
   }
 }
+
+void dump_tgroup_imstats (int Npass) {
+
+  FILE *fout = NULL;
+  char filename[64];
+
+  snprintf (filename, 64, "imstats.grp.%02d.dat", Npass);
+  fout = fopen (filename, "w");
+
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      fprintf (fout, "%.0f %5d | %6ld %6d %6d %4ld %4ld | 0x%08x %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n",
+               ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode,
+               tgroup[j].Nmeasure, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, tgroup[j].Nimage, tgroup[j].Nmosaic,
+               tgroup[j].flags, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].dMsys, tgroup[j].stdev, tgroup[j].dMmin, tgroup[j].dMmax, tgroup[j].McalChiSq);
+    }
+  }
+  fclose (fout);
+
+  // ************ mosaic stats
+
+  snprintf (filename, 64, "imstats.mos.%02d.dat", Npass);
+  fout = fopen (filename, "w");
+
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+
+      for (off_t k = 0; k < tgroup->Nmosaic; k++) {
+        Mosaic *mosaic = tgroup[j].mosaic[k];
+
+        fprintf (fout, "%4ld %3ld %3ld %.0f %5d %5d %5d | %12.5f %5.3f 0x%08x | %6.3f %6.3f %6.3f %6.3f\n",
+                 i, j, k, ohana_sec_to_mjd(tgroupTimes[i][0].start), mosaic->photcode, mosaic->nFitPhotom, mosaic->nValPhotom,
+                 ohana_sec_to_mjd(mosaic->start), mosaic->secz, mosaic->flags,
+                 mosaic->McalPSF, mosaic->dMcal, mosaic->dMsys, mosaic->McalChiSq);
+      }
+    }
+  }
+  fclose (fout);
+
+  // ************ image stats
+
+  snprintf (filename, 64, "imstats.chp.%02d.dat", Npass);
+  fout = fopen (filename, "w");
+
+  off_t Nimage = 0;
+  Image *image = getimages (&Nimage, NULL);
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+
+      for (off_t k = 0; k < tgroup->Nimage; k++) {
+        off_t im = tgroup[j].image[k];
+
+        fprintf (fout, "%4ld %3ld %3ld %.0f %5d | %12.5f %5.3f 0x%08x | %6.3f %6.3f %6.3f %6.3f %6.3f\n",
+                 i, j, k, ohana_sec_to_mjd(tgroupTimes[i][0].start), image[im].photcode,
+                 ohana_sec_to_mjd(image[im].tzero), image[im].secz, image[im].flags,
+                 image[im].McalPSF, image[im].McalAPER, image[im].dMcal, image[im].dMagSys, image[im].McalChiSq);
+      }
+    }
+  }
+
+  fclose (fout);
+}

trunk/Ohana/src/relphot/src/ZeroPointModes.c

@@ -4 +4 @@
 
 void SetZptIteration (int current) {
+  // XXX add 'manual' option here: ask for desired zpt iteration
+  if (MANUAL_ITERATION) {
+    int newValue;
+    fprintf (stdout, "Set Iteration (current = %d): ", CurrentLoop);
+    fscanf (stdin, "%d", &newValue);
+    fprintf (stdout, "Using iteration %d\n", newValue);
+    CurrentLoop = newValue;
+    return;
+  }
+
+  // automatic iteration:
   CurrentLoop = current;
 }
@@ -14 +25 @@
 
   if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
-    if ((CurrentLoop > 4) && (CurrentLoop <= 8)) {
+    if (CurrentLoop < 6) {
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+    }
+    if ((CurrentLoop >= 6) && (CurrentLoop < 12)) {
       // after iterating a few times on the TGroup zero points:
       // * identify the photometric nights
@@ -23 +39 @@
       clean_tgroups(); // do this on every pass to update the status of nights 
     }
-    if ((CurrentLoop > 8) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop >= 12) && (CurrentLoop < 18)) {
       // after iterating a few times on the TGroup & Mosaic zero points:
       // * identify the good / bad Mosaics
@@ -33 +49 @@
       clean_mosaics(); // do this on every pass 
     }
+    if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) {
+      // after iterating a few times on the TGroup & Mosaic zero points:
+      // * identify the good / bad Mosaics
+      // * fit the chips from the bad Mosaics
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights
+      clean_tgroups(); // do this on every pass 
+      clean_mosaics(); // do this on every pass 
+      clean_images(); // do this on every pass 
+    }
   }
 
   if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
+    if (CurrentLoop <= 4) {
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+    }
     if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
       // after iterating a few times on the TGroup zero points:
@@ -48 +80 @@
 
   if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
+    if (CurrentLoop <= 4) {
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE;
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
+    }
     if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
       // after iterating a few times on the Mosaic zero points:
@@ -60 +97 @@
 
   if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
-    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
-      // after iterating a few times on the Image zero points:
-      // * ???
-      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE;
-      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
-      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_ALL;
-    }
+    // if we are not fitting TGROUP or MOSAIC, just fit all images
+    TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE;
+    MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+    IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_ALL;
   }
 }
@@ -83 +117 @@
         return FALSE;
       case ZPT_MOSAIC:
+        if (CurrentLoop < 9) return TRUE;
+        return FALSE;
+      case ZPT_IMAGES:
+      default:
+        if (CurrentLoop < 15) return TRUE;
+        return FALSE;
+    }
+  }
+
+  if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
+    switch (mode) {
+      case ZPT_STARS: 
+        if (CurrentLoop < 3) return TRUE;
+        return FALSE;
+      case ZPT_TGROUP: 
+      case ZPT_MOSAIC:
+        if (CurrentLoop < 3) return TRUE;
+        return FALSE;
+      case ZPT_IMAGES:
+      default:
         if (CurrentLoop < 7) return TRUE;
         return FALSE;
-      case ZPT_IMAGES:
-      default:
-        if (CurrentLoop < 11) return TRUE;
-        return FALSE;
-    }
-  }
-
-  if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
+    }
+  }
+
+  if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
     switch (mode) {
       case ZPT_STARS: 
@@ -108 +158 @@
   }
 
-  if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
-    switch (mode) {
-      case ZPT_STARS: 
-        if (CurrentLoop < 3) return TRUE;
-        return FALSE;
-      case ZPT_TGROUP: 
-      case ZPT_MOSAIC:
-        if (CurrentLoop < 3) return TRUE;
-        return FALSE;
-      case ZPT_IMAGES:
-      default:
-        if (CurrentLoop < 7) return TRUE;
-        return FALSE;
-    }
-  }
-
   if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
     switch (mode) {
@@ -139 +173 @@
   return FALSE;
 }
+
+/*  TGROUP & MOSAIC:
+    0 - TGROUP ONLY & OLS
+    1 - TGROUP ONLY & OLS
+    2 - TGROUP ONLY & OLS
+
+    3 - TGROUP ONLY & IRLS
+    4 - TGROUP ONLY & IRLS
+    5 - TGROUP ONLY & IRLS
+
+    6 - TGROUP + MOSAIC & OLS MOSAIC
+    7 - TGROUP + MOSAIC & OLS MOSAIC
+    8 - TGROUP + MOSAIC & OLS MOSAIC
+
+    9 - TGROUP + MOSAIC & IRLS MOSAIC
+   10 - TGROUP + MOSAIC & IRLS MOSAIC
+   11 - TGROUP + MOSAIC & IRLS MOSAIC
+
+   12 - TGROUP + MOSAIC + IMAGE & OLS IMAGE
+   13 - TGROUP + MOSAIC + IMAGE & OLS IMAGE
+   14 - TGROUP + MOSAIC + IMAGE & OLS IMAGE
+
+   15 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE
+   16 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE
+   17 - TGROUP + MOSAIC + IMAGE & IRLS IMAGE
+
+*/
+

trunk/Ohana/src/relphot/src/args.c

@@ -192 +192 @@
   }
 
+  MANUAL_ITERATION = FALSE;
+  if ((N = get_argument (argc, argv, "-manual-iteration"))) {
+    remove_argument (N, &argc, argv);
+    MANUAL_ITERATION = TRUE;
+  }
+
   NGRID = 8;
   if ((N = get_argument (argc, argv, "-ngrid"))) {
@@ -325 +331 @@
     USE_MCAL_PSF_FOR_STACK_APER = TRUE;
   }
-
 
   // default behavior is to fit all chips independently

trunk/Ohana/src/relphot/src/liststats.c

@@ -26 +26 @@
   stats->Upper80 = NAN;
   stats->Lower20 = NAN;
+  stats->Upper90 = NAN;
+  stats->Lower10 = NAN;
   stats->total   = NAN;
   stats->Nmeas = 0;
@@ -49 +51 @@
   if (N < 1) return (FALSE);
 
+  int N90 = MIN (N-1, 0.9*N);
   int N80 = MIN (N-1, 0.8*N);
   int N20 = MAX (0, 0.2*N);
+  int N10 = MAX (0, 0.1*N);
 
   if (weight) {
@@ -68 +72 @@
   stats[0].Upper80 = value[N80];
   stats[0].Lower20 = value[N20];
+  stats[0].Upper90 = value[N90];
+  stats[0].Lower10 = value[N10];
 
   switch (stats->statmode) {
@@ -116 +122 @@
   WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
   WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
+  if (!dvalue) WeightedMean = FALSE;  // warn the user?
 
   /* calculating the per-star offset based on the desired weighting scheme */
@@ -137 +144 @@
         M  += value[i] * weight[i];
         W  +=            weight[i];
-        dM += SQ (weight[i] * dvalue[i]);
+        dM += dvalue ? SQ (weight[i] * dvalue[i]) : 0.0;
         R  += weight[i];
         Nm += 1.0;  
       } 
       Mo  = M / W;
-      dMo = sqrt (dM) / R;
+      dMo = dvalue ? sqrt (dM) / R : NAN;
     }
   } else {
@@ -159 +166 @@
       for (i = ks; i < ke; i++) {
         M   += value[i];
-        dM  += SQ (dvalue[i]);
+        dM  += dvalue ? SQ (dvalue[i]) : 0.0;
         Nm  += 1.0;  
       } 
       Mo = M / Nm;
-      dMo = sqrt (dM) / Nm;
+      dMo = dvalue ? sqrt (dM) / Nm : NAN;
     }
   }

trunk/Ohana/src/relphot/src/load_images.c

@@ -58 +58 @@
   // images / detections may have been added
   for (off_t i = 0; i < Nsubset; i++) {
-    if (subset[i].flags & ID_IMAGE_PHOTOM_FEW) { 
-      subset[i].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
-    }
+    // reset these two regardless (we will re-determine)
+    subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW);
     if (RESET) {
       if (RESET_ZEROPTS) {
@@ -70 +69 @@
         }
       }
-      subset[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
+      subset[i].flags &= ~(ID_IMAGE_MOSAIC_POOR | ID_IMAGE_NIGHT_POOR);
+      subset[i].flags &= ~(ID_IMAGE_IMAGE_PHOTCAL | ID_IMAGE_MOSAIC_PHOTCAL | ID_IMAGE_TGROUP_PHOTCAL);
       subset[i].ubercalDist = 1000;
     }

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

@@ -117 +117 @@
     for (i = 0; i < NLOOP; i++) {
       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);
@@ -129 +134 @@
       dump_tgroups (catalog, i + 40);
       dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog
+      dump_tgroup_imstats (i + 40);
 
       setMgrid (catalog, flatcorr);
@@ -143 +149 @@
 
       global_stats (catalog, Ncatalog, flatcorr, i);
-
       SetZeroPointModes ();
+
       MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
       fprintf (stderr, "----- continue loop %d -----\n", i);

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -272 +272 @@
     float Map  = NAN;
     float dMap = NAN;
+
+    // XXX I am not convinced the sign on Mflat here is correct
+    float Moff =  Mflat + Mcal + Mmos + Mgrp + Mgrid;
+
     if (isSetMrelFinal) {
       Map = PhotCat (&measure[k], MAG_CLASS_APER);
@@ -282 +286 @@
       if (CHECK_VALID_MAG(Map, dMap)) {
         int Nap = results->aperData[Nsec].Nlist;
-        results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
+        results->aperData[Nsec].flxlist[Nap] = Map - Moff; // this is consistent with PhotRel
         results->aperData[Nsec].errlist[Nap] = dMap;
         results->aperData[Nsec].wgtlist[Nap] = modifiedWeight;
@@ -302 +306 @@
     if (CHECK_VALID_MAG(Mkron, dMkron)) {
       int Nkron = results->kronData[Nsec].Nlist;
-      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
+      results->kronData[Nsec].flxlist[Nkron] = Mkron - Moff; // this is consistent with PhotRel
       results->kronData[Nsec].errlist[Nkron] = dMkron;
       results->kronData[Nsec].wgtlist[Nkron] = modifiedWeight;
@@ -315 +319 @@
     if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
       int Npsf = results->psfData[Nsec].Nlist;
-      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
+      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Moff; // this is consistent with PhotRel
       results->psfData[Nsec].errlist[Npsf] = dMpsf;
       results->psfData[Nsec].wgtlist[Npsf] = modifiedWeight;