Changeset 41462 for trunk

Timestamp:

Jan 4, 2021, 4:33:51 PM (6 years ago)

Author:

eugene

Message:

substantial updates to the zero point statistics analysis

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (3 diffs)
• include/relphot.h (modified) (11 diffs)
• src/GridOps.c (modified) (3 diffs)
• src/ImageOps.c (modified) (5 diffs)
• src/MosaicOps.c (modified) (22 diffs)
• src/StarOps.c (modified) (3 diffs)
• src/TGroupOps.c (modified) (42 diffs)
• src/ZeroPointModes.c (added)
• src/args.c (modified) (3 diffs)
• src/extra.c (modified) (1 diff)
• src/fit1d_irls.c (modified) (3 diffs)
• src/global_stats.c (modified) (5 diffs)
• src/liststats.c (modified) (1 diff)
• src/relphot_images.c (modified) (3 diffs)
• src/relphot_parallel_images.c (modified) (2 diffs)
• src/setMrelCatalog.c (modified) (12 diffs)
• src/test_fitirls.c (modified) (3 diffs)
• test/relphot.tgroups.dvo (modified) (1 diff)

trunk/Ohana/src/relphot/Makefile

@@ -77 +77 @@
 $(SRC)/indexCatalog.$(ARCH).o    \
 $(SRC)/client_logger.$(ARCH).o   \
+$(SRC)/ZeroPointModes.$(ARCH).o  \
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
@@ -124 +125 @@
 $(SRC)/synthetic_zpts.$(ARCH).o  \
 $(SRC)/client_logger.$(ARCH).o   \
+$(SRC)/ZeroPointModes.$(ARCH).o  \
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
@@ -143 +145 @@
 TEST_FITIRLS = \
 $(SRC)/fit1d_irls.$(ARCH).o      \
+$(SRC)/extra.$(ARCH).o  \
 $(SRC)/test_fitirls.$(ARCH).o
 

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

@@ -84 +84 @@
   float secz;
   float ubercalDist;
-  unsigned short nFitPhotom;
+  unsigned int nFitPhotom;
+  unsigned int nValPhotom;
   unsigned int flags;
   char skipCal;               // if TRUE, this mosaic is incomplete and should not be calibrated
+  char inTGroup;
   Coords coords;
 } Mosaic; 
@@ -97 +99 @@
   float dMcal;
   float dMsys;
+  float stdev;
   float McalChiSq;
-  unsigned short nFitPhotom;
+  unsigned int nFitPhotom;
+  unsigned int nValPhotom;
   unsigned int flags;
 
@@ -109 +113 @@
   off_t *measure;
   off_t *catalog;
+
+  off_t Nmosaic;
+  off_t NMOSAIC;
+  Mosaic **mosaic; // pointer to the mosaic structures
 
   void *parent;
@@ -185 +193 @@
   int     Nlist;              // total number of measurements in list (may be more than the number to use)
   int     Nbootstrap;         // number of bootstrap iterations
+  int     MaxIterations;      // number of IRLS iterations
 
   // input data & parameters: 
@@ -387 +396 @@
 int    VERBOSE2;
 int    VERBOSE_IMAGE;
+
 int    MOSAIC_ZEROPT;
 int    TGROUP_ZEROPT;
+
+int    FREEZE_IMAGES;
+int    FREEZE_MOSAICS;
+
 int    TGROUP_FIT_AIRMASS;
 
@@ -414 +428 @@
 double IMAGE_GOOD_FRACTION;
 int    IMAGE_BAD;
-int    FREEZE_IMAGES;
-int    FREEZE_MOSAICS;
 int    CALIBRATE_STACKS_AND_WARPS;
 int    USE_GRID;
@@ -461 +473 @@
 SkyRegion UserPatch;
 char *UserCatalog;
+
+enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, };
+int TGROUP_ZPT_MODE;
+
+enum {MOSAIC_ZPT_MODE_NONE, MOSAIC_ZPT_MODE_BAD_NIGHT, MOSAIC_ZPT_MODE_GOOD_MOSAIC, MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC, MOSAIC_ZPT_MODE_ALL, };
+int MOSAIC_ZPT_MODE;
+
+enum {IMAGE_ZPT_MODE_NONE, IMAGE_ZPT_MODE_BAD_NIGHT, IMAGE_ZPT_MODE_BAD_MOSAIC, IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC, IMAGE_ZPT_MODE_ALL, };
+int IMAGE_ZPT_MODE;
 
 int USE_ALL_IMAGES;
@@ -526 +547 @@
 float         getMgrid            PROTO((off_t meas, int cat));
 float         getMmos             PROTO((off_t meas, int cat));
-float         getMgrp             PROTO((off_t meas, int cat, float airmass));
+float         getMgrp             PROTO((off_t meas, int cat, float airmass, float *dZpt));
 float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source));
 short         getUbercalDist      PROTO((off_t meas, int cat));
@@ -601 +622 @@
 int           reload_images       PROTO((FITS_DB *db));
 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
-void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
+void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
 void          setMcalFromMosaics  PROTO((void));
 void          setMcalFromTGroups  PROTO((void));
 int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
-void          setMgrid            PROTO((Catalog *catalog, 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, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
+int           setMmos             PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
+int           setMgrp             PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
+
 int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
 void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
+
 int           setMave             PROTO((Catalog *catalog, int Ncatalog));
 void          set_ZP              PROTO((double ZERO));
@@ -627 +651 @@
 StatType      statsMosaicX        PROTO((Catalog *catalog));
 StatType      statsMosaicdM       PROTO((Catalog *catalog));
+StatType      statsTGroupM        PROTO((Catalog *catalog));
+StatType      statsTGroupX        PROTO((Catalog *catalog));
+StatType      statsTGroupdM       PROTO((Catalog *catalog));
+
 StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr));
 StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog, int Nsec));
@@ -795 +823 @@
 
 int save_images_backup (FITS_DB *db);
+
+void dump_tgroups (Catalog *catalog, int Npass);
+void dump_catalog (Catalog *catalog, off_t c, int Npass);
+
+void SetZeroPointModes (int nloop);

trunk/Ohana/src/relphot/src/GridOps.c

@@ -415 +415 @@
 
         // skip mosaics marked as BAD
-        Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+        Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
         if (isnan(Mgrp)) {
           Ngrp ++;
@@ -575 +575 @@
         continue;
       }
-      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) {
         Ngrp ++;
@@ -661 +661 @@
       Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) continue;
-      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
       Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);

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

@@ -581 +581 @@
 
 /* determine Mcal values for all images */
-void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+void setMcal (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   off_t i, j, m, c, n;
   int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
 
-  // FREEZE_IMAGES only applies to mosaic data (eg, gpc1)
+  if (IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_NONE) return;
 
   fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
@@ -609 +609 @@
   FitDataSetAlloc (&brightStars, Nmax, 0, 0);
 
+  // until the analysis has converged a bit, do not use the IRLS analysis
+  if (Nloop < 3) {
+    brightStars.MaxIterations = 0;
+    kronStars.MaxIterations = 0;
+    psfStars.MaxIterations = 0;
+  }
+
+
   Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0;
 
@@ -627 +635 @@
     }      
 
-    // only freeze the gpc1 chips
-    // XXX : when do we FREEZE_IMAGES in a mix of GPC1, GPC2, etc?
-    // if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
-    if (FREEZE_IMAGES) continue;
+    // if requested, freeze mosaic chips:
+    if (FREEZE_IMAGES && isMosaicChip (image[i].photcode)) continue;
     if (image[i].photcode == 0) continue; // skip the PHU images
+
+    int badNight  = (image[i].flags & ID_IMAGE_NIGHT_POOR);
+    int badMosaic = (image[i].flags & ID_IMAGE_MOSAIC_POOR);
+
+    // in BAD_NIGHT mode, we fit ONLY images in bad nights
+    if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_NIGHT)) {
+      if (!badNight) continue;
+    }
+
+    // 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) continue;
+    }
+
+    // in BAD_MOSAIC mode, we fit bad mosaics ignoring night state
+    if ((IMAGE_ZPT_MODE == IMAGE_ZPT_MODE_BAD_MOSAIC)) {
+      if (badMosaic) continue;
+    }
 
     // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
@@ -665 +690 @@
           continue;
       }
-      float Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      float Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL); // ignore error for now?
       if (isnan(Mgrp)) {
           Ngrp ++;
@@ -1027 +1052 @@
       Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) continue;
-      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
       Mgrid = getMgrid (m, c);

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

@@ -10 +10 @@
 // relationships between the mosaics and their associated images
 static off_t   *MosaicN_Image; // number of images associated with the given mosaic
-static off_t  **MosaicToImage; // list of imagesa associated with the given mosaic
+static off_t  **MosaicToImage; // list of images associated with the given mosaic
 
 // mosaic index for given image : ImageToMosaic[ImageIndex] = MosaicIndex (ImageIndex : 0 < Nimage)
@@ -219 +219 @@
     mosaic[i].photcode  = 0;
     mosaic[i].skipCal   = FALSE;
+    mosaic[i].inTGroup  = FALSE; // not (yet?) assigned to a TGroup
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -301 +302 @@
     mosaic[j].dMcal     = 0.0;
     mosaic[j].McalChiSq = 0.0;
+
+    mosaic[j].nFitPhotom = 0;
+    mosaic[j].nValPhotom = 0;
 
     mosaic[j].dMsys     = subset[i].flags;
@@ -750 +754 @@
 }
 
+void markBadMosaic (int myMosaic) {
+
+  off_t Nimage;
+  Image *image;
+
+  // this mosaic has been identified as poor.
+  mosaic[myMosaic].flags |= ID_IMAGE_PHOTOM_POOR;
+  mosaic[myMosaic].McalPSF  = 0.0;
+  mosaic[myMosaic].McalAPER = 0.0;
+
+  image = getimages (&Nimage, NULL);
+
+  // all images should be marked as coming from a bad mosaic
+  // these will be fitted independently
+  for (off_t i = 0; i < MosaicN_Image[myMosaic]; i++) {
+    off_t im = MosaicToImage[myMosaic][i];
+    image[im].flags |= ID_IMAGE_MOSAIC_POOR;
+  }
+}
+
+void markGoodMosaic (int myMosaic) {
+
+  off_t Nimage;
+  Image *image;
+
+  // this mosaic has been identified as good.
+  mosaic[myMosaic].flags &= ~ID_IMAGE_PHOTOM_POOR;
+
+  image = getimages (&Nimage, NULL);
+
+  // all images should be marked as NOT coming from a bad mosaic
+  for (off_t i = 0; i < MosaicN_Image[myMosaic]; i++) {
+    off_t im = MosaicToImage[myMosaic][i];
+    image[im].flags &= ~ID_IMAGE_MOSAIC_POOR;
+  }
+}
+
 void initMosaicBins (Catalog *catalog, int Ncatalog, int doMosaicList) {
 
@@ -923 +964 @@
   int Nsys;
   int Nskip;
+  int Nloop;
   off_t Nmax;
   int PoorImages;
@@ -943 +985 @@
 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, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr);
+
+void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) {
   info->Nfew = 0;
   info->Nbad = 0;
@@ -955 +997 @@
 
   info->Nmax = Nmax;
+  info->Nloop = Nloop;
   info->PoorImages = PoorImages;
 
@@ -962 +1005 @@
     FitDataSetAlloc (&info->kronStars,   Nmax, 0, 0); 
     FitDataSetAlloc (&info->brightStars, Nmax, 0, 0); 
+
+    // until the analysis has converged a bit, do not use the IRLS analysis
+    if (Nloop < 3) {
+      info->brightStars.MaxIterations = 0;
+      info->kronStars.MaxIterations = 0;
+      info->psfStars.MaxIterations = 0;
+    }
   }
 }
@@ -1003 +1053 @@
 }
 
-int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMmos (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   off_t i, N, Nmax;
   Image *image;
 
-  if (!MOSAIC_ZEROPT) return (FALSE);
-  if (FREEZE_MOSAICS) return (FALSE);
+  if (!MOSAIC_ZEROPT) return FALSE;
+  if (FREEZE_MOSAICS) return FALSE;
+
+  if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_NONE) return FALSE;
 
   // plots cannot be done in a threaded context (the plot commands collide)
   // 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, Nloop, PoorImages, flatcorr);
     return status;
   }
@@ -1034 +1086 @@
 
   SetMmosInfo info;
-  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
+  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages, Nloop);
 
   // int savePlotDelay = PLOTDELAY;
@@ -1078 +1130 @@
   assert (Nmos >= 0);
   assert (Nmos < Nmosaic);
+
+  int badNight  = (myMosaic[0].flags & ID_IMAGE_NIGHT_POOR);
+  int badMosaic = (myMosaic[0].flags & ID_IMAGE_MOSAIC_POOR);
+
+  // 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
+  if (MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_BAD_NIGHT) {
+    if (!badNight) return TRUE;
+  }
+
+  // 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;
+  }
+
+  // in GOOD_MOSAIC mode, we fit good mosaics ignoring night state
+  if ((MOSAIC_ZPT_MODE == MOSAIC_ZPT_MODE_GOOD_MOSAIC)) {
+    if (badMosaic) return TRUE;
+  }
 
   // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
@@ -1140 +1217 @@
     if (fout) {
       float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
-      // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      // 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);
@@ -1164 +1241 @@
       continue;
     }
-    float Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+    float Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
     if (isnan(Mgrp)) {
       info->Ngrp++;
@@ -1257 +1334 @@
   myMosaic[0].McalChiSq  = psfStars->chisq;         
   myMosaic[0].nFitPhotom = psfStars->Nmeas;         
+  myMosaic[0].nValPhotom = N;       
 
   // fit1d_irls (kronStars, N);    // does it make sense to calculate a zero point based kron mags?
@@ -1296 +1374 @@
 }
   
-int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMmos_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   int i;
@@ -1317 +1395 @@
 
   SetMmosInfo summary;
-  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
+  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop);
 
   pthread_attr_t attr;
@@ -1340 +1418 @@
     threadinfo[i].image    =    image;
     threadinfo[i].flatcorr = flatcorr;
-    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
+
+    // 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, Nloop);
     pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
   }
@@ -1398 +1479 @@
 
   SetMmosInfo results;
-  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
+  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop);
 
   while (1) {
@@ -1709 +1790 @@
       Mcal = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) continue;
-      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
       Mgrid = getMgrid (m, c);
@@ -1767 +1848 @@
 void clean_mosaics () {
 
-  off_t i, N, mark, Nmark, Nscatter, Noffset;
+  off_t i, N, mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq;
   double *mlist, *slist, *dlist;
   double MaxOffset, MedOffset, MaxScatter;
@@ -1800 +1881 @@
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
   
-  Nmark = Nscatter = Noffset = 0;
+  // XXX this needs to be configurable
+  float MaxChiSq = 10.0;
+
+  Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 0;
   for (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;
 
-    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
+    Ntotal ++;
+
+    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW)) {
+      continue;
+    }
     if (mosaic[i].skipCal) continue;
 
@@ -1817 +1905 @@
       Noffset ++;
     }
+    if (mosaic[i].nFitPhotom < 0.7*mosaic[i].nValPhotom) {
+      mark = TRUE;
+      Nrejects ++;
+    }
+    if (mosaic[i].McalChiSq > MaxChiSq) {
+      mark = TRUE;
+      Nchisq ++;
+    }
     if (mark) { 
       Nmark ++;
-      mosaic[i].flags |= ID_IMAGE_PHOTOM_POOR;
+      markBadMosaic(i); // mark the images associated with a bad night
     } else {
-      mosaic[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
-    }
-  }
-
-  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
+      markGoodMosaic(i); // mark the images associated with a good night
+    }
+  }
+
+  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);
 
   free (mlist);

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

@@ -493 +493 @@
         float Mmos  = getMmos  (m, i);
         if (isnan(Mmos)) continue;
-        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
+        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
         if (isnan(Mgrp)) continue;
         float Mgrid = getMgrid (m, i);
@@ -733 +733 @@
           Mmos  = getMmos  (m, i);
           if (isnan(Mmos)) { Nmos ++; continue; }
-          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
+          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
           if (isnan(Mgrp)) { Ngrp ++; continue; }
           Mgrid = getMgrid (m, i);
@@ -1082 +1082 @@
 }
 
+void dump_catalog (Catalog *catalog, off_t c, int Npass) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  FILE *fout = NULL;
+  char filename[64];
+  snprintf (filename, 64, "tcat.%02d.dat", Npass);
+  fout = fopen (filename, "w");
+
+  for (off_t n = 0; n < catalog[c].Naverage; n++) {
+
+    off_t m = catalog[c].averageT[n].measureOffset;
+    for (off_t k = 0; k < catalog[c].averageT[n].Nmeasure; k++, m++) {
+
+      float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);     // image zero point
+      float Mmos     = getMmos  (m, c);                    // mosaic zero point
+      float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
+
+      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
+      // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+      // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+
+      off_t n = catalog[c].measureT[m].averef;
+
+      // nominal magnitude for this measurement (instrumental + Ko(airmass - 1) + Co
+      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 - Mgrp;
+
+      fprintf (fout, "%f %f : %f %f : %f %f %f  : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta);
+    }
+  }
+  fclose (fout);
+}

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

@@ -1 +1 @@
 # include "relphot.h"
+
+void assignMosaicsToTGroups (void);
 
 // tgroupTimes carries the times of the tgroups (initially, just the photometric nights)
@@ -78 +80 @@
       tgroup[j].dMsys     = 0.0;
       tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
+
+      tgroup[j].stdev      = NAN;
+      tgroup[j].nFitPhotom = 0;
+      tgroup[j].nValPhotom = 0;
+
       tgroup[j].flags     = ID_IMAGE_PHOTOM_UBERCAL; // we start by treating the suspected photometric nights as photometric
       tgroup[j].photcode  = photcodes[j][0].code;
@@ -85 +92 @@
       ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE);
       tgroup[j].image[0]  = -1;
+
+      tgroup[j].NMOSAIC   = 100;
+      tgroup[j].Nmosaic   = 0;
+      ALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC);
 
       tgroup[j].Nmeasure  = 0;
@@ -129 +140 @@
     tgroup[j].dMsys     = 0.0;
     tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
+
+    tgroup[j].stdev      = NAN;
+    tgroup[j].nFitPhotom = 0;
+    tgroup[j].nValPhotom = 0;
+
     tgroup[j].flags     = 0;   // tgroups NOT in the original file are not suspected as photometric
     tgroup[j].photcode  = photcodes[j][0].code;
+
     tgroup[j].NIMAGE    = 1000;
     tgroup[j].Nimage    = 0;   // the new TGroup is not suspected to be photometric
     ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE);
     tgroup[j].image[0]  = -1;
+
+    tgroup[j].NMOSAIC   = 100;
+    tgroup[j].Nmosaic   = 0;
+    ALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC);
+
     tgroup[j].Nmeasure  = 0;
     tgroup[j].NMEASURE  = 1000;
@@ -199 +221 @@
   }
 
+  // assign the mosaics to each tgroup -- we do this by image
+  assignMosaicsToTGroups ();
+
   initTGroupsMcal ();
 
@@ -215 +240 @@
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
       free (tgroup[j].image);
+      free (tgroup[j].mosaic);
       free (tgroup[j].measure);
       free (tgroup[j].catalog);
@@ -249 +275 @@
       tgroup[j].dMsys     = 0.0;
       tgroup[j].McalChiSq = 0.0;// NAN or 0.0?
+
+      tgroup[j].stdev      = NAN;
+      tgroup[j].nFitPhotom = 0;
+      tgroup[j].nValPhotom = 0;
 
       if (tgroup[j].Nimage == 0) continue; // no images, ignore the tgroup
@@ -278 +308 @@
 }
 
+void assignMosaicsToTGroups (void) {
+
+  if (!TGROUP_ZEROPT) return;
+
+  // init the tgroup array values
+  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].Nimage == 0) continue; // no images, ignore the tgroup
+
+      // find the mosaics by image (if MOSAIC_ZEROPT is active)
+      for (int im = 0; MOSAIC_ZEROPT && (im < tgroup[j].Nimage); im++) {
+        int imageIdx = tgroup[j].image[im];
+        Mosaic *mosaic = getMosaicForImage (imageIdx);
+
+        if (mosaic->inTGroup) continue;
+
+        // mosaic to tgroup and extend as needed
+        tgroup[j].mosaic[tgroup[j].Nmosaic] = mosaic;
+        mosaic->inTGroup = TRUE;
+
+        // advance array and extend if needed
+        tgroup[j].Nmosaic ++;
+        if (tgroup[j].Nmosaic >= tgroup[j].NMOSAIC) {
+          tgroup[j].NMOSAIC += 100;
+          REALLOCATE (tgroup[j].mosaic, Mosaic *, tgroup[j].NMOSAIC);
+        }
+      }
+    }
+  }
+  return;
+}
+
 TGroup *getTGroupForImage (off_t im) {
 
@@ -345 +409 @@
 
   fprintf (stderr, "*** return Mcal from tgroup.Mcal to image.Mcal ***\n");
+
+  // NOTE: nights which are poor should have had Mcal set to 0, making this a NOP for those nights
 
   // copy the tgroup results to the images.  set the tgroup Mcal to 0.0 since we have moved its
@@ -360 +426 @@
         image[im].dMagSys     = tgroup[i].dMsys;
         image[im].nFitPhotom  = tgroup[i].nFitPhotom;
+        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // XXX copy this down?
+        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // XXX copy this down?
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW);
         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
@@ -371 +439 @@
     }
   }      
+}
+
+void markBadTGroup (TGroup *tgroup) {
+
+  off_t Nimage;
+  Image *image;
+
+  // this tgroup has been identified as poor.
+  tgroup->flags |= ID_IMAGE_NIGHT_POOR;
+
+  // all associated mosaics should be marked as coming from a bad night
+  // these will be fitted independently
+  for (off_t j = 0; j < tgroup->Nmosaic; j++) {
+    tgroup->mosaic[j]->flags |= ID_IMAGE_NIGHT_POOR;
+  }
+
+  // all associated images should be marked as coming from a bad night
+  // these will be fitted independently
+  image = getimages (&Nimage, NULL);
+  for (off_t j = 0; j < tgroup->Nimage; j++) {
+    off_t im = tgroup->image[j];
+    image[im].flags |= ID_IMAGE_NIGHT_POOR;
+  }
+
+  tgroup->McalPSF  = 0.0;
+  tgroup->McalAPER = 0.0;
+  tgroup->dKlam    = 0.0;
+}
+
+void markGoodTGroup (TGroup *tgroup) {
+
+  off_t Nimage;
+  Image *image;
+
+  // this tgroup has been identified as good.
+  tgroup->flags &= ~ID_IMAGE_NIGHT_POOR;
+
+  // all associated mosaics should be marked as NOT coming from a bad night
+  for (off_t j = 0; j < tgroup->Nmosaic; j++) {
+    tgroup->mosaic[j]->flags &= ~ID_IMAGE_NIGHT_POOR;
+  }
+
+  // all associated images should be marked as NOT coming from a bad night
+  image = getimages (&Nimage, NULL);
+  for (off_t j = 0; j < tgroup->Nimage; j++) {
+    off_t im = tgroup->image[j];
+    image[im].flags &= ~ID_IMAGE_NIGHT_POOR;
+  }
 }
 
@@ -488 +604 @@
 }
 
-float getMgrp (off_t meas, int cat, float airmass) {
-
-  float value;
-
-  if (!TGROUP_ZEROPT) return (0);
+// XXX return the error as well as the value
+float getMgrp (off_t meas, int cat, float airmass, float *dZpt) {
+
+  if (dZpt) { *dZpt = 0; }
+
+  if (!TGROUP_ZEROPT) return 0.0;
 
   // unassigned measurements belong to simple chips
   TGroup *myGroup = MeasureToTGroup[cat][meas];
-  if (!myGroup) return (0.0);
-
-  // if (myGroup->flags & IMAGE_BAD) return (NAN);  
-  value = myGroup->McalPSF + myGroup->dKlam*(airmass - 1.0);
-  return (value);
+  if (!myGroup) return 0.0;
+
+  float value = myGroup->McalPSF + myGroup->dKlam*(airmass - 1.0);
+  if (dZpt) { *dZpt = myGroup->dMcal; }
+  return value;
 }
 
@@ -525 +642 @@
   int Nsys;
   int Nskip;
+  int Nloop;
   off_t Nmax;
   int PoorImages;
@@ -545 +663 @@
 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, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr);
+
+void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) {
   info->Nfew = 0;
   info->Nbad = 0;
@@ -558 +676 @@
 
   info->Nmax = Nmax;
+  info->Nloop = Nloop;
   info->PoorImages = PoorImages;
 
@@ -567 +686 @@
     FitDataSetAlloc (&info->brightStars, Nmax, fitOrder, 0); 
 
+    // until the analysis has converged a bit, do not use the IRLS analysis
+    if (Nloop < 100) {
+      info->brightStars.MaxIterations = 0;
+      info->kronStars.MaxIterations = 0;
+      info->psfStars.MaxIterations = 0;
+    }
     if (1) {
       FitDataSetAddPriors (&info->psfStars);
@@ -614 +739 @@
 }
 
-int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMgrp (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   off_t N;
@@ -620 +745 @@
 
   if (!TGROUP_ZEROPT) return (FALSE);
+  if (TGROUP_ZPT_MODE == TGROUP_ZPT_MODE_NONE) return (FALSE);
 
   // plots cannot be done in a threaded context (the plot commands collide)
   // 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, Nloop, PoorImages, flatcorr);
     return status;
   }
 
   image = getimages (&N, NULL);
-
-  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;
@@ -645 +765 @@
 
   SetMgrpInfo info;
-  SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages);
-
-  // int savePlotDelay = PLOTDELAY;
-  // if (PLOTSTUFF) PLOTDELAY = 0.0;
+  SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages, Nloop);
 
   for (off_t i = 0; i < NtgroupTimes; i++) {
@@ -654 +771 @@
     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, 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, 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);
     }
   }
@@ -668 +785 @@
   fprintf (stderr, "%d tgroups 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);
 }
@@ -684 +796 @@
   off_t j, NimageReal;
 
-  /* on PoorImages run, skip good images */
-  if (info->PoorImages) {
-    int bad = myTGroup->flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
-    if (!bad) return TRUE;
-  }      
+  // 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;
+  }
 
   FitDataSet *psfStars = &info->psfStars;
@@ -704 +817 @@
   int testImage = FALSE;
   TGTimes *tgroup = (TGTimes *) myTGroup->parent;
-  testImage |= (abs(tgroup->start - 1245283200) < 10);
+  // testImage |= (abs(tgroup->start - 1245283200) < 10);
+  // testImage = TRUE;
 
   FILE *fout = NULL;
@@ -825 +939 @@
 
   /* too few good measurements or too many bad measurements (skip in PoorImages run) */
-  if (!info->PoorImages) {
-    int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure);
-    if (mark) {
-      if (VERBOSE2) {
-        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);
-      }
-      myTGroup->flags |= ID_IMAGE_PHOTOM_FEW;
-      info->Nfew ++;
-      if (testImage) {
-        fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure));
-      }
-    } else {
-      myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW;
-    }
+  int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure);
+  if (mark) {
+    if (VERBOSE2) {
+      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);
+    }
+    myTGroup->flags |= ID_IMAGE_PHOTOM_FEW;
+    info->Nfew ++;
+    if (testImage) {
+      fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure));
+    }
+  } else {
+    myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW;
   }
 
@@ -852 +964 @@
   myTGroup->McalPSF    = psfStars->bSaveArray[0][0];
   myTGroup->dMcal      = psfStars->bSigma[0];
+  myTGroup->stdev      = psfStars->sigma;
   myTGroup->McalChiSq  = psfStars->chisq;
   myTGroup->nFitPhotom = psfStars->Nmeas;
+  myTGroup->nValPhotom = N;
   myTGroup->dKlam      = psfStars->bSaveArray[1][0];
 
@@ -864 +978 @@
     fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom);
   }
+
+  // fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.2f %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, myTGroup->nValPhotom, myTGroup->McalChiSq, ohana_sec_to_mjd(tgroup->start));
 
   if (PLOTSTUFF) {
@@ -885 +1001 @@
   }
 
-  // FitDataSetFree (&brightStars);
-  // FitDataSetFree (&kronStars);
-  // FitDataSetFree (&psfStars);
-
   return TRUE;
 }
   
-int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
+int setMgrp_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   int i;
@@ -901 +1013 @@
   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;
   for (off_t i = 0; i < NtgroupTimes; i++) {
@@ -916 +1022 @@
 
   SetMgrpInfo summary;
-  SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages);
+  SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop);
 
   pthread_attr_t attr;
@@ -939 +1045 @@
     threadinfo[i].image    =    image;
     threadinfo[i].flatcorr = flatcorr;
-    SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
+
+    // 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
+    SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages, Nloop);
     pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]);
   }
@@ -986 +1095 @@
   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;
 }
@@ -999 +1103 @@
 
   SetMgrpInfo results;
-  SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
+  SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop); // allocate list, dlist arrays here
 
   while (1) {
@@ -1024 +1128 @@
 }
 
-# if (0)
+// XXX should I split these stats out by photcode?
+// XXX double check the IMAGE_BAD filter:
 StatType statsTGroupM (Catalog *catalog) {
   OHANA_UNUSED_PARAM(catalog);
 
-  off_t i;
-  int n;
   double *list, *dlist;
   StatType stats;
@@ -1037 +1140 @@
 
   // XXX do this by time & photcode?
-  ALLOCATE (list, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-
-  n = 0;
+  ALLOCATE (list, double, NtgroupTimes*Nphotcodes);
+  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 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 & NIGHT_BAD) continue;
       list[n] = tgroup[j].McalPSF;
       dlist[n] = 1;
@@ -1062 +1168 @@
   OHANA_UNUSED_PARAM(catalog);
 
-  off_t i, n;
   double *list, *dlist;
   StatType stats;
@@ -1069 +1174 @@
   if (!TGROUP_ZEROPT) return (stats);
 
-  ALLOCATE (list, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-
-  n = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    if (tgroup[i].flags & IMAGE_BAD) continue;
-    list[n] = tgroup[i].dMcal;
-    dlist[n] = 1;
-    n++;
+  ALLOCATE (list, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
+
+  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;
+      list[n] = tgroup[j].dMcal;
+      dlist[n] = 1;
+      n++;
+    }
   }
 
@@ -1088 +1196 @@
 }
 
-StatType statsTGroupN (Catalog *catalog, FlatCorrectionTable *flatcorr) {
-
-  off_t i, j, m, c, n, N;
+StatType statsTGroupX (Catalog *catalog) {
+  OHANA_UNUSED_PARAM(catalog);
+
   double *list, *dlist;
-  float Mcal, Mgrid, Mrel, Mmos;
   StatType stats;
 
@@ -1098 +1205 @@
   if (!TGROUP_ZEROPT) return (stats);
 
-  ALLOCATE (list, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-
-  n = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    if (tgroup[i].flags & IMAGE_BAD)  continue;
-
-    N = 0;
-    for (j = 0; j < tgroup[i].Nmeasure; j++) {
-
-      m = tgroup[i].measure[j];
-      c = tgroup[i].catalog[j];
-
-      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
-      if (isnan(Mcal)) continue;
-      Mmos = getMmos  (m, c);
-      if (isnan(Mmos)) continue;
-      Mgrid = getMgrid (m, c);
-      if (isnan(Mgrid)) continue;
-      Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      if (isnan(Mrel)) continue;
-      N++;
-    }
-    list[n] = N;
-    dlist[n] = 1;
-    n++;
-  }
-  // fprintf (stderr, "Ntgroup: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Ntgroup,  n);
+  ALLOCATE (list, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
+
+  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;
+      list[n] = tgroup[j].McalChiSq;
+      dlist[n] = 1;
+      n++;
+    }
+  }
 
   liststats_setmode (&stats, "MEAN");
@@ -1135 +1227 @@
 }
 
-StatType statsTGroupX (Catalog *catalog) {
-  OHANA_UNUSED_PARAM(catalog);
-
-  off_t i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (!TGROUP_ZEROPT) return (stats);
-
-  ALLOCATE (list, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-
-  n = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    if (tgroup[i].flags & IMAGE_BAD) continue;
-    if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
-
-    list[n]  = tgroup[i].McalChiSq;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats_setmode (&stats, "MEAN");
-
-  liststats (list, dlist, NULL, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-# endif
-
 /* mark tgroup if: abs(Mcal - <Mcal>) too large, dMcal too large */
 void clean_tgroups () {
 
-  off_t i, N, mark, Nmark, Nscatter, Noffset;
+  off_t mark, Ntotal, Nmark, Nscatter, Noffset, Nrejects, Nchisq;
   double *mlist, *slist, *dlist;
   double MaxOffset, MedOffset, MaxScatter;
@@ -1182 +1242 @@
   ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
 
-  for (i = N = 0; i < NtgroupTimes; i++) {
+  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++) {
@@ -1197 +1258 @@
 
   liststats (mlist, dlist, NULL, N, &stats);
-  MaxOffset = MAX (NIGHT_OFFSET, 2*stats.sigma);
+  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, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
+  fprintf (stderr, "TGROUPS: median Mrel: %f, median dMrel: %f, Max Offset: %f, Max Scatter: %f\n", MedOffset, stats.median, MaxOffset, MaxScatter);
   
-  Nmark = Nscatter = Noffset = 0;
-  for (i = 0; i < NtgroupTimes; i++) {
+  // XXX this needs to be configurable
+  float MaxChiSq = 10.0;
+
+  Ntotal = Nmark = Nscatter = Noffset = Nrejects = Nchisq = 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 we are keeping ubercal sacrosanct, then we should not be allowed to break them...
-      if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
       
-      if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
+      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
+      if (tgroup[j].flags & ID_IMAGE_PHOTOM_FEW) {
+        markBadTGroup(&tgroup[j]); // mark the image associated with a bad night
+        continue;
+      }
       
       mark = FALSE;
-      if (tgroup[i].dMcal > MaxScatter) {
+      if (tgroup[j].dMcal > MaxScatter) {
         mark = TRUE;
         Nscatter ++;
       }
-      if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) {
+      if (fabs(tgroup[j].McalPSF - MedOffset) > MaxOffset) {
         mark = TRUE;
         Noffset ++;
       }
+      if (tgroup[j].nFitPhotom < 0.7*tgroup[j].nValPhotom) {
+        mark = TRUE;
+        Nrejects ++;
+      }
+      if (tgroup[j].McalChiSq > MaxChiSq) {
+        mark = TRUE;
+        Nchisq ++;
+      }
       if (mark) { 
         Nmark ++;
-        tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR;
+        markBadTGroup(&tgroup[j]); // mark the tgroup & associated images & mosaics with a bad night
       } else {
-        tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
+        markGoodTGroup(&tgroup[j]); // mark the tgroup & associated images & mosaics with a good night
       }
-      // XXX : unset the ubercal flag if we mark the night as bad?
-    }
-  }
-
-  fprintf (stderr, OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
+    }
+  }
+
+  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);
 
   free (mlist);
@@ -1347 +1424 @@
 }
 
+void dump_tgroups (Catalog *catalog, int Npass) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int Ngrp = 0;
+  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++) {
+
+      FILE *fout = NULL;
+      char filename[64];
+      snprintf (filename, 64, "tgrp.%05d.%02d.dat", Ngrp, Npass);
+      fout = fopen (filename, "w");
+      Ngrp ++;
+
+      TGroup *myTGroup = &tgroup[j];
+
+      for (int k = 0; k < myTGroup->Nmeasure; k++) {
+      
+        off_t m = myTGroup->measure[k];
+        off_t c = myTGroup->catalog[k];
+      
+        float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);     // image zero point
+        float Mmos     = getMmos  (m, c);                    // mosaic zero point
+        float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
+
+        float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
+        // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+        // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+
+        off_t n = catalog[c].measureT[m].averef;
+
+        // nominal magnitude for this measurement (instrumental + Ko(airmass - 1) + Co
+        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 - Mgrp;
+
+        fprintf (fout, "%f %f : %f %f : %f %f %f  : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta);
+      }
+      fclose (fout);
+    }
+  }
+}

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

@@ -303 +303 @@
   }
 
+  // global restriction on fitting image zero points for mosaic cameras
+  FREEZE_IMAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-imfreeze"))) {
+    remove_argument (N, &argc, argv);
+    FREEZE_IMAGES = TRUE;
+  }
+  // global restriction on fitting mosaic zero points
+  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;
+  }
+  USE_MCAL_PSF_FOR_STACK_APER = FALSE;
+  if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) {
+    remove_argument (N, &argc, argv);
+    USE_MCAL_PSF_FOR_STACK_APER = TRUE;
+  }
+
+
+  // default behavior is to fit all chips independently
+  TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE;
+  MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+  IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_ALL;
+
   /* XXX this argument used to do two things: specify the camera name and tell the analysis to
      calculate a common zero point for a single mosaic.  I've moved the MOSAICNAME concept into the
@@ -315 +344 @@
       exit (2);
     }
-  }
-
-  FREEZE_IMAGES = FALSE;
-  if ((N = get_argument (argc, argv, "-imfreeze"))) {
-    remove_argument (N, &argc, argv);
-    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;
-  }
-  USE_MCAL_PSF_FOR_STACK_APER = FALSE;
-  if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) {
-    remove_argument (N, &argc, argv);
-    USE_MCAL_PSF_FOR_STACK_APER = TRUE;
+    IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
+    MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL;
   }
 
@@ -341 +351 @@
   TGROUP_ZEROPT = FALSE;
   char *TGROUP_FILENAME = NULL;
-  if ((N = get_argument (argc, argv, "-tgroup"))) {
+  if ((N = get_argument (argc, argv, "-tgroup")) || (N = get_argument (argc, argv, "-tgroups"))) {
     TGROUP_ZEROPT = TRUE;
     remove_argument (N, &argc, argv);
     TGROUP_FILENAME = strcreate(argv[N]);
     remove_argument (N, &argc, argv);
-    FREEZE_IMAGES = TRUE;
-  }
-  if ((N = get_argument (argc, argv, "-tgroups"))) {
-    if (TGROUP_ZEROPT) {
-      fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n");
-      exit (2);
-    }
-    TGROUP_ZEROPT = TRUE;
-    remove_argument (N, &argc, argv);
-    TGROUP_FILENAME = strcreate(argv[N]);
-    remove_argument (N, &argc, argv);
-    FREEZE_IMAGES = TRUE;
+    // if we fit for TGROUPS, start with the images and mosaics frozen
+    TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_ALL;
+    IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
+    MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+  }
+  if ((N = get_argument (argc, argv, "-tgroup")) || (N = get_argument (argc, argv, "-tgroups"))) {
+    fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n");
+    exit (2);
   }
   TGROUP_FIT_AIRMASS = FALSE;

trunk/Ohana/src/relphot/src/extra.c

@@ -130 +130 @@
 }
 
+double weight_cauchy (double x) {
+  double r = x / 2.385;
+  return (1.0 / (1.0 + SQ(r)));
+}
+
+double VectorFractionInterpolate (double *values, float fraction, int Npts) {
+
+  float F = fraction * Npts;
+  int   N = fraction * Npts;
+
+  if (N < 0        ) return NAN;
+  if (N >= Npts - 2) return NAN;
+
+  // interpolate between N,N+1
+    
+  double S = (F - N) * (values[N+1] - values[N]) + values[N];
+  return S;
+}

trunk/Ohana/src/relphot/src/fit1d_irls.c

@@ -2 +2 @@
 
 // These should probably be tunable:
-# define MAX_ITERATIONS 10
 # define FIT_TOLERANCE 1e-4
 # define FLT_TOLERANCE 1e-6
@@ -35 +34 @@
 
   int converged = FALSE;
-  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
+  for (int iterations = 0; !converged && (iterations < dataset->MaxIterations); iterations++) {
     
     for (int i = 0; i < Npoints; i++) {
@@ -293 +292 @@
   dataset->mterm = 2*order + 1;
   dataset->Nbootstrap = Nbootstrap;
+  dataset->MaxIterations = 10;
 
   // arrays for actual data values

trunk/Ohana/src/relphot/src/global_stats.c

@@ -3 +3 @@
 void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop) {
 
-  StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD;
+  StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD, tgM, tgX, tgD;
 
   liststats_init (&stN);
   liststats_init (&stX);
   liststats_init (&stS);
+
   liststats_init (&imN);
   liststats_init (&imX);
   liststats_init (&imM);
   liststats_init (&imD);
+
   liststats_init (&msN);
   liststats_init (&msX);
@@ -17 +19 @@
   liststats_init (&msD);
 
-  // INITTIME;
+  liststats_init (&tgX);
+  liststats_init (&tgM);
+  liststats_init (&tgD);
 
   fprintf (stderr, "\n");
@@ -29 +33 @@
     int seccode = photcodes[Ns][0].code;
 
-    // MARKTIME("start star stats: %f sec\n", dtime);
     stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr);
-    // MARKTIME("done star N: %f sec\n", dtime);
     stX = statsStarX (catalog, Ncatalog, Nsec);
-    // MARKTIME("done star X: %f sec\n", dtime);
     stS = statsStarS (catalog, Ncatalog, Nsec);
-    // MARKTIME("done star S: %f sec\n", dtime);
   
     fprintf (stderr, "   --- stats for %s ---\n", photcodes[Ns][0].name);
@@ -53 +53 @@
   msX = statsMosaicX (catalog);
   
+  tgM = statsMosaicM (catalog);
+  tgD = statsMosaicdM (catalog);
+  tgX = statsMosaicX (catalog);
+  
   fprintf (stderr, "meas / image:   %7.0f  %7.0f  %7.0f  %7.0f  %7.0f  %6d\n",   imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas);
   fprintf (stderr, "Mcal image:     %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas);
@@ -62 +66 @@
   fprintf (stderr, "dMcal mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
   fprintf (stderr, "chisq mosaic:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);
+
+  fprintf (stderr, "Mcal  tgroup:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   tgM.median, tgM.mean, tgM.sigma, tgM.min, tgM.max, tgM.Nmeas);
+  fprintf (stderr, "dMcal tgroup:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   tgD.median, tgD.mean, tgD.sigma, tgD.min, tgD.max, tgD.Nmeas);
+  fprintf (stderr, "chisq tgroup:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   tgX.median, tgX.mean, tgX.sigma, tgX.min, tgX.max, tgX.Nmeas);
 }
 

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

@@ -417 +417 @@
 }
 
-double weight_cauchy (double x) {
-  double r = x / 2.385;
-  return (1.0 / (1.0 + SQ(r)));
-}
-
-double VectorFractionInterpolate (double *values, float fraction, int Npts) {
-
-  float F = fraction * Npts;
-  int   N = fraction * Npts;
-
-  if (N < 0        ) return NAN;
-  if (N >= Npts - 2) return NAN;
-
-  // interpolate between N,N+1
-    
-  double S = (F - N) * (values[N+1] - values[N]) + values[N];
-  return S;
-}
-
 // we could define the weight to be the only scale factor:
 // \mu      = \sum (value_i * weight_i) / \sum (weight_i)

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

@@ -116 +116 @@
     /* determine fit values */
     for (i = 0; i < NLOOP; i++) {
-      if (PLOTSTUFF) {
-        plot_scatter (catalog, Ncatalog, flatcorr); 
-      }
-      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
-      // checkImages("start loop");
+      dump_tgroups (catalog, i);
+      dump_catalog (catalog, 0, i); // for a test, just dump a specific catalog
 
       setMrel  (catalog, Ncatalog, flatcorr); // threaded
-      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
-      // checkImages("set Mrel");
-
-      if (PLOTSTUFF) {
-        plot_scatter (catalog, Ncatalog, flatcorr); 
-      }
-      setMcal  (catalog, FALSE, flatcorr);
-      // checkImages("set Mcal");
-      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
-
-      setMmos  (catalog, FALSE, flatcorr);
-      setMgrp  (catalog, FALSE, flatcorr);
+      dump_tgroups (catalog, i + 20);
+      dump_catalog (catalog, 0, i + 20); // for a test, just dump a specific catalog
+
+      setMcal  (catalog, i, FALSE, flatcorr);
+      setMmos  (catalog, i, FALSE, flatcorr);
+      setMgrp  (catalog, i, FALSE, flatcorr);
+      dump_tgroups (catalog, i + 40);
+      dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog
 
       setMgrid (catalog, flatcorr);
@@ -147 +140 @@
         plot_chisq (catalog, Ncatalog);
       }
-      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
-      // if (i % 6 == 1) rationalize_images ();
-
-      // NOTE : in the past, I was not iterating enough before cleaning.  make sure we do
-      // at least 8 loops first -- that should get the systematic errors down to the ~1%
-      // level, even in cases where we have an even split between photometric data and
-      // data with 1 mag of extinction.
-
-      if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); 
-      if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
-      if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
-      if ((i > 8) && (i % 8 == 5)) clean_tgroups ();
-      if ((i > 8) && (i % 8 == 5)) clean_images ();
 
       global_stats (catalog, Ncatalog, flatcorr, i);
+
+      SetZeroPointModes (i);
       MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
+      fprintf (stderr, "----- continue loop %d -----\n", i);
     }
 
@@ -178 +161 @@
     /* set Mcal & Mmos for bad images */
     int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
-    setMcal  (catalog, onlyPoorImages, flatcorr);
-    setMmos  (catalog, onlyPoorImages, flatcorr);
-    setMgrp  (catalog, onlyPoorImages, flatcorr);
+    onlyPoorImages = FALSE; // XXX deactivate for now
+    setMcal  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
+    setMmos  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
+    setMgrp  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 

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

@@ -114 +114 @@
 
     // 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, i, FALSE, flatcorr);
+    setMmos  (catalog, i, FALSE, 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);
@@ -151 +151 @@
   /* 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);
+  setMcal  (catalog, 999, onlyPoorImages, flatcorr);
+  setMmos  (catalog, 999, onlyPoorImages, flatcorr);
   // setMgrp  (catalog, onlyPoorImages, flatcorr); XXX see above
   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);

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

@@ -9 +9 @@
 # endif
 
+int magStatsByRanking (StatDataSet *dataset, StatType *stats);
 int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats);
 int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats);
@@ -17 +18 @@
 int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats);
 
-# define MAG_STATS_BY_RANKING magStatsByRankingIRLS
+// # define MAG_STATS_BY_RANKING magStatsByRankingIRLS
+# define MAG_STATS_BY_RANKING magStatsByRanking
 
 # define UBERCAL_WEIGHT 100.0
@@ -56 +58 @@
 }
 
+FILE *TESTFILE = NULL;
+
 int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
@@ -67 +71 @@
   SetMrelInfoReset (results); // reset the counters
 
+  TESTFILE = fopen ("mytest.dat", "w");
+  
   for (j = 0; j < catalog[Nc].Naverage; j++) {
 
@@ -83 +89 @@
     }
   }
+
+  fclose (TESTFILE);
   return (TRUE);
 }
@@ -222 +230 @@
     // ** Choose the calibration (depends on the mode : do I have an image reference or not?) 
     float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
+    float dMgrp = 0.0;
     off_t meas = measureOffset + k;
     if (getImageEntry (meas, cat) < 0) {
@@ -235 +244 @@
       Mmos  = getMmos  (meas, cat);
       if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
-      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass);
+      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
       if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
       Mgrid = getMgrid (meas, cat);
@@ -251 +260 @@
     int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));
 
+    // combine the various errors in quadrature here:
+    float dMsys = hypot (code->photomErrSys, dMgrp);
+
     float Map  = NAN;
     float dMap = NAN;
@@ -258 +270 @@
       // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
       // use PSF flux error instead here:
-      // dMap = MAX (hypot(measure[k].dMap,   code->photomErrSys), MIN_ERROR);
-      dMap = MAX (hypot(measure[k].dM,   code->photomErrSys), MIN_ERROR);
+      // dMap = MAX (hypot(measure[k].dMap,   dMsys), MIN_ERROR);
+      dMap = MAX (hypot(measure[k].dM,   dMsys), MIN_ERROR);
 
       if (CHECK_VALID_MAG(Map, dMap)) {
@@ -277 +289 @@
     float dMkron;
     if (isSetMrelFinal) {
-      dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
+      dMkron = MAX (hypot(measure[k].dMkron, dMsys), MIN_ERROR);
     } else {
-      dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+      dMkron = MAX (hypot(measureT[k].dM, dMsys), MIN_ERROR);
     }
     if (CHECK_VALID_MAG(Mkron, dMkron)) {
@@ -293 +305 @@
 
     float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
-    float dMpsf  = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+    float dMpsf  = MAX (hypot(measureT[k].dM, dMsys), MIN_ERROR);
     if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
       int Npsf = results->psfData[Nsec].Nlist;
@@ -401 +413 @@
       secfilt[Nsec].dMpsfChp = psfstats->error;
       secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
+      fprintf (TESTFILE, "%d %d : %d : %f %f : %f\n", cat, (int) ave, NrankingPSF, psfstats->mean, psfstats->median, psfstats->mean - psfstats->median);
     }
 

trunk/Ohana/src/relphot/src/test_fitirls.c

@@ -5 +5 @@
 int main (int argc, char **argv) {
   
-  int Ntests     =  100;
-  int Npoints    =  100;
-  int Noutliers  =   10;
-  int Nbootstrap =  100;
+  int Ntests      =  100;
+  int Npoints     =  100;
+  int Noutliers   =   10;
+  int Nbootstrap  =  100;
+  int Niterations =   10;
 
   int N;
@@ -29 +30 @@
     remove_argument (N, &argc, argv);
     Nbootstrap = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Niterations"))) {
+    remove_argument (N, &argc, argv);
+    Niterations = atoi (argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -88 +94 @@
   FitDataSet dataset;
   FitDataSetAlloc (&dataset, Npoints + Noutliers, FitSlope, Nbootstrap);
+  dataset.MaxIterations = Niterations;
 
   if (isfinite(ValuePrior) || isfinite(SlopePrior)) {

trunk/Ohana/src/relphot/test/relphot.tgroups.dvo

@@ -451 +451 @@
   set stars_dmag = 0.02 + 0.1 * rnd(tmp)
 end
+
+macro plot.test.files
+  if ($0 != 2)
+    echo "USAGE: plot.test.files (N)"
+    break
+  end
+
+  clear -s
+  resize 2000 1200
+
+  section a0 0.00 0.50 0.25 0.50
+  section a1 0.25 0.50 0.25 0.50
+  section a2 0.50 0.50 0.25 0.50
+  section a3 0.75 0.50 0.25 0.50
+
+  section a4 0.00 0.00 0.25 0.50
+  section a5 0.25 0.00 0.25 0.50
+  section a6 0.50 0.00 0.25 0.50
+  section a7 0.75 0.00 0.25 0.50
+
+  label -fn helvetica 12
+  $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5
+
+  for i 0 8
+    sprintf name tgrp.%05d.%02d.dat $i $1
+    data $name
+    exec ls -l --full-time $name
+    read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12
+  
+    section a$i
+    lim Msys dMag; box $BSTYLE; plot Msys dMag -pt cir -op 0.2 -c blue
+  end
+end