Changeset 40477 for branches/czw_branch/20170908/Ohana/src/relphot

Timestamp:

Jun 27, 2018, 3:20:35 PM (8 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/czw_branch/20170908/Ohana

Files:

• . (modified) (1 prop)
• src/relphot/include/relphot.h (modified) (7 diffs)
• src/relphot/src/BrightCatalog.c (modified) (14 diffs)
• src/relphot/src/GridOps.c (modified) (4 diffs)
• src/relphot/src/ImageMagIO.c (modified) (8 diffs)
• src/relphot/src/ImageOps.c (modified) (23 diffs)
• src/relphot/src/ImageSubset.c (modified) (9 diffs)
• src/relphot/src/MosaicOps.c (modified) (30 diffs)
• src/relphot/src/StarOps.c (modified) (15 diffs)
• src/relphot/src/args.c (modified) (4 diffs)
• src/relphot/src/bcatalog.c (modified) (1 diff)
• src/relphot/src/launch_region_hosts.c (modified) (1 diff)
• src/relphot/src/liststats.c (modified) (5 diffs)
• src/relphot/src/load_images.c (modified) (1 diff)
• src/relphot/src/plot_scatter.c (modified) (2 diffs)
• src/relphot/src/plotstuff.c (modified) (1 diff)
• src/relphot/src/reload_catalogs.c (modified) (2 diffs)
• src/relphot/src/relphot_images.c (modified) (6 diffs)
• src/relphot/src/relphot_objects.c (modified) (1 diff)
• src/relphot/src/relphot_parallel_images.c (modified) (4 diffs)
• src/relphot/src/relphot_parallel_regions.c (modified) (2 diffs)
• src/relphot/src/relphot_synthphot_catalog.c (modified) (1 diff)
• src/relphot/src/select_images.c (modified) (1 diff)
• src/relphot/src/setMrelCatalog.c (modified) (19 diffs)
• src/relphot/src/setMrelFinal.c (modified) (3 diffs)
• src/relphot/src/share_image_mags.c (modified) (2 diffs)
• src/relphot/src/share_mean_mags.c (modified) (2 diffs)
• src/relphot/src/synthetic_mags.c (modified) (1 diff)

branches/czw_branch/20170908/Ohana/src/relphot/include/relphot.h

@@ -77 +77 @@
   unsigned int stop;
   short photcode;
-  float Mcal;
+  float McalPSF;
+  float McalAPER;
   float dMcal;
   float dMsys;
-  unsigned short nFitPhotom;
-  short Xm;
+  float McalChiSq;
   float secz;
   float ubercalDist;
+  unsigned short nFitPhotom;
   unsigned int flags;
   char skipCal;               // if TRUE, this mosaic is incomplete and should not be calibrated
@@ -201 +202 @@
 
 typedef struct {
-  float Mcal;
+  float McalPSF;
+  float McalAPER;
   float dMcal;
   float dMagSys;
-  short Xm;
+  float McalChiSq;
   int nFitPhotom;
   int flags;
@@ -220 +222 @@
 
 typedef struct {
-  float Mcal;
+  float McalPSF;
+  float McalAPER;
   float dMcal;
   unsigned int imageID;
@@ -291 +294 @@
 int    PARALLEL_SERIAL;
 char  *MANUAL_UNIQUER;
+
+int    SKIP_PARALLEL_GROUPS;
 
 int    PARALLEL_REGIONS_MANUAL;
@@ -336 +341 @@
 int    SyntheticPhotometry;
 
+int    USE_MCAL_PSF_FOR_STACK_APER;
+
 char  *PhotcodeList;
 
@@ -404 +411 @@
 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
 
+void clearImages (void);
+void checkImages (char *name);
+int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog);
+
 void makeMosaics (Image *image, off_t Nimage, int mergeMcal);
 Mosaic *getMosaicForImage (off_t im);
@@ -422 +433 @@
 off_t         getImageEntry       PROTO((off_t meas, int cat));
 
-float         getMcal_alt         PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd));
-float         getMcal             PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
+float         getMcal             PROTO((off_t meas, int cat, dvoMagClassType class));
 float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
 float         getMgrid            PROTO((off_t meas, int cat));
 float         getMmos             PROTO((off_t meas, int cat));
-float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
+float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source));
 short         getUbercalDist      PROTO((off_t meas, int cat));
 float         getCenterOffset     PROTO((off_t meas, int cat, Measure *measure, unsigned int *myID));

branches/czw_branch/20170908/Ohana/src/relphot/src/BrightCatalog.c

@@ -66 +66 @@
     GET_COLUMN(D,         "DEC",        double);
     GET_COLUMN(M,         "MAG_SYS",    float);
-    GET_COLUMN(Mcal,      "MAG_CAL",    float);
+    GET_COLUMN(Mkron,     "MAG_KRON",   float);
+    GET_COLUMN(McalPSF,   "MCAL_PSF",   float);
+    GET_COLUMN(McalAPER,  "MCAL_APER",  float);
     GET_COLUMN(Mflat,     "MAG_FLAT",   float);
     GET_COLUMN(dM,        "MAG_ERR",    float);
@@ -89 +91 @@
       measure[i].D         = D[i];
       measure[i].M         = M[i];
-      measure[i].Mcal      = Mcal[i];
+      measure[i].Mkron     = Mkron[i];
+      measure[i].McalPSF   = McalPSF[i];
+      measure[i].McalAPER  = McalAPER[i];
       measure[i].Mflat     = Mflat[i];
       measure[i].dM        = dM[i];
@@ -110 +114 @@
     free (D       );
     free (M       );
-    free (Mcal    );
+    free (Mkron   );
+    free (McalPSF );
+    free (McalAPER);
     free (Mflat   );
     free (dM      );
@@ -186 +192 @@
     // need to create and assign to flat-field correction
     GET_COLUMN(M,      "MAG",      float);
+    GET_COLUMN(Mkron,  "MAG_KRON", float);
     GET_COLUMN(dM,     "MAG_ERR",  float);
     GET_COLUMN(Mchisq, "MAG_CHI",  float);
@@ -199 +206 @@
     ALLOCATE (secfilt, SecFilt, Nrow);
     for (i = 0; i < Nrow; i++) {
-      secfilt[i].M      = M[i];         
-      secfilt[i].dM     = dM[i];
-      secfilt[i].Mchisq = Mchisq[i];
-      secfilt[i].flags  = flags[i];
-      secfilt[i].Ncode  = Ncode[i];
-      secfilt[i].Nused  = Nused[i];
-      secfilt[i].Mmin   = Mmin[i];
-      secfilt[i].Mmax   = Mmax[i];
+      secfilt[i].MpsfChp  = M[i];         
+      secfilt[i].MkronChp = Mkron[i];         
+      secfilt[i].dMpsfChp = dM[i];
+      secfilt[i].Mchisq   = Mchisq[i];
+      secfilt[i].flags    = flags[i];
+      secfilt[i].Ncode    = Ncode[i];
+      secfilt[i].Nused    = Nused[i];
+      secfilt[i].Mmin     = Mmin[i];
+      secfilt[i].Mmax     = Mmax[i];
     }
     fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow);
 
     free (M     );
+    free (Mkron );
     free (dM    );
     free (Mchisq);
@@ -271 +280 @@
     gfits_define_bintable_column (&theader, "D", "DEC",      "dec",                        "degree", 1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_KRON", "magnitude (sys,kron)",        NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MCAL_PSF", "magnitude (cal)",             NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MCAL_APER","magnitude (cal)",             NULL,    1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_FLAT", "magnitude (flat)",            NULL,    1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "magnitude (err)",             NULL,    1.0, 0.0);
@@ -293 +304 @@
     double *D         ; ALLOCATE (D        ,  double, catalog->Nmeasure);
     float  *M         ; ALLOCATE (M        ,  float,  catalog->Nmeasure);
-    float  *Mcal      ; ALLOCATE (Mcal     ,  float,  catalog->Nmeasure);
+    float  *Mkron     ; ALLOCATE (Mkron    ,  float,  catalog->Nmeasure);
+    float  *McalPSF   ; ALLOCATE (McalPSF  ,  float,  catalog->Nmeasure);
+    float  *McalAPER  ; ALLOCATE (McalAPER ,  float,  catalog->Nmeasure);
     float  *Mflat     ; ALLOCATE (Mflat    ,  float,  catalog->Nmeasure);
     float  *dM        ; ALLOCATE (dM       ,  float,  catalog->Nmeasure);
@@ -311 +324 @@
     MeasureTiny *measure = catalog->measure;
     for (i = 0; i < catalog->Nmeasure; i++) {
-      R[i]        = measure[i].R       ;
-      D[i]        = measure[i].D       ;
+      R[i]        = measure[i].R        ;
+      D[i]        = measure[i].D        ;
       M[i]        = measure[i].M        ;
-      Mcal[i]     = measure[i].Mcal     ;
+      Mkron[i]    = measure[i].Mkron    ;
+      McalPSF[i]  = measure[i].McalPSF  ;
+      McalAPER[i] = measure[i].McalAPER ;
       Mflat[i]    = measure[i].Mflat    ;
       dM[i]       = measure[i].dM       ;
@@ -334 +349 @@
     gfits_set_bintable_column (&theader, &ftable, "DEC",        D,         catalog->Nmeasure);
     gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_KRON",   Mkron,     catalog->Nmeasure);
+    gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF",   McalPSF,   catalog->Nmeasure);
+    gfits_set_bintable_column (&theader, &ftable, "MCAL_APER",  McalAPER,  catalog->Nmeasure);
     gfits_set_bintable_column (&theader, &ftable, "MAG_FLAT",   Mflat,     catalog->Nmeasure);
     gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",    dM,        catalog->Nmeasure);
@@ -352 +369 @@
     free (D       );
     free (M       );
-    free (Mcal    );
+    free (Mkron   );
+    free (McalPSF );
+    free (McalAPER);
     free (Mflat   );
     free (dM      );
@@ -436 +455 @@
     gfits_create_table_header (&theader, "BINTABLE", "SECFILT");
 
-    gfits_define_bintable_column (&theader, "E", "MAG",      "",              "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "",              "arcsec", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG",      "",              "mag",   1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_KRON", "",              "mag",   1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "",              "mag",   1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_CHI",  "",              NULL,    1.0, 0.0);
     gfits_define_bintable_column (&theader, "J", "FLAGS",    "",              NULL,    1.0, 0.0);
@@ -454 +474 @@
     // create intermediate storage arrays
     float *M        ; ALLOCATE (M      ,  float, Nsec);
+    float *Mkron    ; ALLOCATE (Mkron  ,  float, Nsec);
     float *dM       ; ALLOCATE (dM     ,  float, Nsec);
     float *Mchisq   ; ALLOCATE (Mchisq ,  float, Nsec);
@@ -465 +486 @@
     SecFilt *secfilt = catalog->secfilt;
     for (i = 0; i < Nsec; i++) {
-      M     [i]       = secfilt[i].M      ;
-      dM    [i]       = secfilt[i].dM     ;
-      Mchisq[i]       = secfilt[i].Mchisq ;
-      flags [i]       = secfilt[i].flags  ;
-      Ncode [i]       = secfilt[i].Ncode  ;
-      Nused [i]       = secfilt[i].Nused  ;
-      Mmin  [i]       = secfilt[i].Mmin   ;
-      Mmax  [i]       = secfilt[i].Mmax   ;
+      M     [i]       = secfilt[i].MpsfChp ;
+      Mkron [i]       = secfilt[i].MkronChp;
+      dM    [i]       = secfilt[i].dMpsfChp;
+      Mchisq[i]       = secfilt[i].Mchisq  ;
+      flags [i]       = secfilt[i].flags   ;
+      Ncode [i]       = secfilt[i].Ncode   ;
+      Nused [i]       = secfilt[i].Nused   ;
+      Mmin  [i]       = secfilt[i].Mmin    ;
+      Mmax  [i]       = secfilt[i].Mmax    ;
     }
 
     // add the columns to the output array
     gfits_set_bintable_column (&theader, &ftable, "MAG",      M     , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_KRON", Mkron , Nsec);
     gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM    , Nsec);
     gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Mchisq, Nsec);
@@ -486 +509 @@
 
     free (M      );
+    free (Mkron  );
     free (dM     );
     free (Mchisq );

branches/czw_branch/20170908/Ohana/src/relphot/src/GridOps.c

@@ -401 +401 @@
 
         // skip images marked as BAD
-        Mcal = getMcal  (m, c, flatcorr, catalog);
+        Mcal = getMcal  (m, c, MAG_CLASS_PSF);
         if (isnan(Mcal)) {
           Ncal ++;
@@ -558 +558 @@
         continue;
       }
-      Mcal = getMcal  (m, c, flatcorr, catalog);
+      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) {
         Ncal ++;
@@ -568 +568 @@
         continue;
       }
-      Mrel  = getMrel  (catalog, m, c);
+      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
       if (isnan(Mrel)) {
         Nrel ++;
@@ -645 +645 @@
         continue;
       }
-      Mcal  = getMcal  (m, c, flatcorr, catalog);
+      Mcal  = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) continue;
       Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) continue;
-      Mrel  = getMrel  (catalog, m, c);
+      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
       if (isnan(Mrel)) continue;
 

branches/czw_branch/20170908/Ohana/src/relphot/src/ImageMagIO.c

@@ -58 +58 @@
   char type[16];
 
-  GET_COLUMN (Mcal,              "MCAL",           float);
+  GET_COLUMN (McalPSF,           "MCAL_PSF",       float);
+  GET_COLUMN (McalAPER,          "MCAL_APER",      float);
   GET_COLUMN (dMcal,             "MCAL_ERR",       float);
   GET_COLUMN (dMagSys,           "MCAL_SYSERR",    float);
-  GET_COLUMN (nFitPhotom,        "NFIT",             int);
-  GET_COLUMN (flags,             "FLAGS",            int);
-  GET_COLUMN (ubercalDist,       "UDIST",            int);
-  GET_COLUMN (imageID,           "ID",               int);
-  GET_COLUMN (Xm,                "CHISQ",          short);
+  GET_COLUMN (McalChiSq,         "MCAL_CHISQ",     float);
+  GET_COLUMN (nFitPhotom,        "NFIT",           int);
+  GET_COLUMN (flags,             "FLAGS",          int);
+  GET_COLUMN (ubercalDist,       "UDIST",          int);
+  GET_COLUMN (imageID,           "ID",             int);
 
   // free the memory associated with the FITS files
@@ -75 +76 @@
   ALLOCATE (image_mags, ImageMag, Nrow);
   for (i = 0; i < Nrow; i++) {
-    image_mags[i].Mcal                      = Mcal       [i];
+    image_mags[i].McalPSF                   = McalPSF    [i];
+    image_mags[i].McalAPER                  = McalAPER   [i];
     image_mags[i].dMcal                     = dMcal      [i];
     image_mags[i].dMagSys                   = dMagSys    [i];
-    image_mags[i].Xm                        = Xm         [i];
+    image_mags[i].McalChiSq                 = McalChiSq  [i];
     image_mags[i].nFitPhotom                = nFitPhotom [i];
     image_mags[i].flags                     = flags      [i];
@@ -86 +88 @@
   fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow);
 
-  free (Mcal       );
+  free (McalPSF    );
+  free (McalAPER   );
   free (dMcal      );
   free (dMagSys    );
-  free (Xm         );
+  free (McalChiSq  );
   free (nFitPhotom );
   free (flags      );
@@ -121 +124 @@
   gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS");
 
-  gfits_define_bintable_column (&theader, "E", "MCAL",           "cal offset",                 "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_PSF",       "PSF cal offset",             "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_APER",      "APER cal offset",            "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MCAL_ERR",       "cal error",                  "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR",    "systematic error",           "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_CHISQ",     "cal chisq",                  "unitless",   1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "NFIT",           "number of fitted stars",     "unitless",   1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "FLAGS",          "analysis flags",             "unitless",   1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "UDIST",          "distance to ubercal images", "images",     1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "ID",             "image ID",                   "unitless",   1.0, 0.0);
-  gfits_define_bintable_column (&theader, "I", "CHISQ",          "cal chisq",                  "unitless",   1.0, FT_BZERO_INT16);
 
   // generate the output array that carries the data
   gfits_create_table (&theader, &ftable);
 
-  float *Mcal        ;
+  float *McalPSF     ;
+  float *McalAPER    ;
   float *dMcal       ;
   float *dMagSys     ;
-  float *Xm          ;
+  float *McalChiSq   ;
   int   *nFitPhotom  ;
   int   *flags       ;
@@ -143 +148 @@
 
   // create intermediate storage arrays
-  ALLOCATE (Mcal        ,         float,          Nimage_mags);
+  ALLOCATE (McalPSF     ,         float,          Nimage_mags);
+  ALLOCATE (McalAPER    ,         float,          Nimage_mags);
   ALLOCATE (dMcal       ,         float,          Nimage_mags);
   ALLOCATE (dMagSys     ,         float,          Nimage_mags);
-  ALLOCATE (Xm          ,         float,          Nimage_mags);
+  ALLOCATE (McalChiSq   ,         float,          Nimage_mags);
   ALLOCATE (nFitPhotom  ,           int,          Nimage_mags);
   ALLOCATE (flags       ,           int,          Nimage_mags);
@@ -154 +160 @@
   // assign the storage arrays
   for (i = 0; i < Nimage_mags; i++) {
-    Mcal       [i]   = image_mags[i].Mcal       ;
+    McalPSF    [i]   = image_mags[i].McalPSF    ;
+    McalAPER   [i]   = image_mags[i].McalAPER   ;
     dMcal      [i]   = image_mags[i].dMcal      ;
     dMagSys    [i]   = image_mags[i].dMagSys    ;
-    Xm         [i]   = image_mags[i].Xm         ;
+    McalChiSq  [i]   = image_mags[i].McalChiSq  ;
     nFitPhotom [i]   = image_mags[i].nFitPhotom ;
     flags      [i]   = image_mags[i].flags      ;
@@ -165 +172 @@
 
   // add the columns to the output array
-  gfits_set_bintable_column (&theader, &ftable, "MCAL",           Mcal       ,         Nimage_mags);
+  gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF",       McalPSF    ,         Nimage_mags);
+  gfits_set_bintable_column (&theader, &ftable, "MCAL_APER",      McalAPER   ,         Nimage_mags);
   gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR",       dMcal      ,         Nimage_mags);
   gfits_set_bintable_column (&theader, &ftable, "MCAL_SYSERR",    dMagSys    ,         Nimage_mags);
-  gfits_set_bintable_column (&theader, &ftable, "CHISQ",          Xm         ,         Nimage_mags);
+  gfits_set_bintable_column (&theader, &ftable, "MCAL_CHISQ",     McalChiSq  ,         Nimage_mags);
   gfits_set_bintable_column (&theader, &ftable, "NFIT",           nFitPhotom ,         Nimage_mags);
   gfits_set_bintable_column (&theader, &ftable, "FLAGS",          flags      ,         Nimage_mags);
@@ -174 +182 @@
   gfits_set_bintable_column (&theader, &ftable, "ID",             imageID    ,         Nimage_mags);
 
-  free (Mcal       );
+  free (McalPSF    );
+  free (McalAPER   );
   free (dMcal      );
   free (dMagSys    );
-  free (Xm         );
+  free (McalChiSq  );
   free (nFitPhotom );
   free (flags      );

branches/czw_branch/20170908/Ohana/src/relphot/src/ImageOps.c

@@ -143 +143 @@
     image[i].photom_map_id = input[i].photom_map_id;
     image[i].flags         = input[i].flags        ;
-    image[i].Mcal          = input[i].Mcal         ;
+    image[i].McalPSF       = input[i].McalPSF      ;
+    image[i].McalAPER      = input[i].McalAPER     ;
     image[i].dMcal         = input[i].dMcal        ;
     image[i].tzero         = input[i].tzero        ;
@@ -179 +180 @@
     subset[i].photom_map_id = image[i].photom_map_id;
     subset[i].flags         = image[i].flags        ;
-    subset[i].Mcal          = image[i].Mcal         ;
+    subset[i].McalPSF       = image[i].McalPSF      ;
+    subset[i].McalAPER      = image[i].McalAPER     ;
     subset[i].dMcal         = image[i].dMcal        ;
     subset[i].tzero         = image[i].tzero        ;
@@ -430 +432 @@
 }
 
-// returns image.Mcal
+// returns image.McalPSF or image.McalAPER
 // NOTE: static flat-field component is included in measure.Mflat
-float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
-  OHANA_UNUSED_PARAM(flatcorr);
-  OHANA_UNUSED_PARAM(catalog);
-
-  off_t i;
-  float value;
+float getMcal (off_t meas, int cat, dvoMagClassType class) {
+
+  off_t i;
 
   i = MeasureToImage[cat][meas];
-  if (i == -1) return (NAN);
-
-  if (image[i].flags & IMAGE_BAD) return (NAN);  
-  value = image[i].Mcal;
-
-  // to do this, I need to pass in the catalog and flatcorr pointers
-  // int flat_id = image[i].photom_map_id;
-  // float offset = 0.0;
-  // if (flat_id) {
-  //   offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
-  // }
-  // value -= offset;
-
-  return (value);
-}
-
-// returns image.Mcal - ff(x,y)
-// NOTE: static flat-field component is included in measure.Mflat
-float getMcal_alt (off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd) {
-  OHANA_UNUSED_PARAM(flatcorr);
-  OHANA_UNUSED_PARAM(Xccd);
-  OHANA_UNUSED_PARAM(Yccd);
-
-  off_t i;
-  float value;
-
-  i = MeasureToImage[cat][meas];
-  if (i == -1) return (NAN);
-
-  if (image[i].flags & IMAGE_BAD) return (NAN);  
-  value = image[i].Mcal;
-
-  // to do this, I need to pass in the catalog and flatcorr pointers
-  // int flat_id = image[i].photom_map_id;
-  // float offset = 0.0;
-  // if (flat_id) {
-  //   offset = FlatCorrectionOffset (flatcorr, flat_id, Xccd, Yccd);
-  // }
-  // value -= offset;
-
-  return (value);
-}
-
-// returns image.Mcal - ff(x,y)
+  if (i == -1) return NAN;
+
+  if (image[i].flags & IMAGE_BAD) return NAN;  
+
+  switch (class) {
+    case MAG_CLASS_PSF:
+      return image[i].McalPSF;
+    case MAG_CLASS_APER:
+    case MAG_CLASS_KRON:
+      return image[i].McalAPER;
+    default:
+      return NAN;
+  }
+  return NAN; // should not be able to reach here
+}
+
 short getUbercalDist (off_t meas, int cat) {
 
@@ -499 +470 @@
 }
 
-// returns image.Mcal - ff(x,y)
 float getCenterOffset (off_t meas, int cat, Measure *measure, unsigned int *myID) {
 
@@ -519 +489 @@
 }
 
-// returns image.Mcal - ff(x,y)
 int MatchImageName (off_t meas, int cat, char *name) {
 
@@ -542 +511 @@
 }
 
-// returns image.Mcal - ff(x,y)
 int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID) {
 
@@ -563 +531 @@
 }
 
-// returns image.Mcal - ff(x,y)
 int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) {
 
@@ -618 +585 @@
 void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
-  off_t i, j, m, c, n, Nmax;
+  off_t i, j, m, c, n;
   int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys;
-  float Msys, Mrel, Mmos, Mgrid, Mflat;
-  // double *list, *dlist, *Mlist, *dMlist;
 
   StatType stats;
@@ -638 +603 @@
   }
 
-  Nmax = 0;
+  off_t Nmax = 0;
   for (i = 0; i < Nimage; i++) {
     Nmax = MAX (Nmax, N_onImage[i]);
   }
-  StatDataSet *refStars = StatDataSetAlloc (1, Nmax);
+
+  StatDataSet *kronStars   = StatDataSetAlloc (1, Nmax);
+  StatDataSet *psfStars    = StatDataSetAlloc (1, Nmax);
   StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
 
@@ -676 +643 @@
     int minUbercalDist = 1000;
     
-    off_t Nref = 0;  // total number of reference stars on the image
-    int Nbright = 0; // number of stars to measure the bright-end scatter
+    off_t Nref = 0;    // number of stars used to measure McalPSF
+    int   Nkron = 0;   // number of stars to measure McalAPER
+    int   Nbright = 0; // number of stars to measure the bright-end scatter
 
     if (N_onImage[i] == 0) {
@@ -693 +661 @@
           continue;
       }
-      Mmos  = getMmos  (m, c);
+      float Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) {
           Nmos ++;
           continue;
       }
-      Mgrid = getMgrid (m, c);
+      float Mgrid = getMgrid (m, c);
       if (isnan(Mgrid)) {
           Ngrid++;
           continue;
       }
-      Mrel  = getMrel  (catalog, m, c);
-      if (isnan(Mrel)) {
+
+      // Mrel* is the average magnitude for this star.  For PS1 stacks, we have too much
+      // PSF variability.  We need to calibrate the PSF magnitudes separately from the
+      // Aperture-like magnitues.  (We have an option to use the kron magnitudes or the
+      // other apertures here).  I basically need to do this analysis separately for each
+      // magnitude type
+
+      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
+      if (isnan(MrelPSF)) {
           Nrel ++;
           continue;
@@ -714 +689 @@
       // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
 
-      Mflat = getMflat (m, c, flatcorr, catalog);
+      float Mflat = getMflat (m, c, flatcorr, catalog);
 
       n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-      if (isnan(Msys)) {
+      float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
+      if (isnan(MsysPSF)) {
         Nsys++;
         continue;
@@ -731 +706 @@
 
     skip:
-      refStars->flxlist[Nref] = Msys - Mrel - Mmos - Mgrid + Mflat;
-      refStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      refStars->wgtlist[Nref] = 1;
-      refStars->msklist[Nref] = 0;
-      if (fabs(refStars->flxlist[Nref]) > 0.03) {
+      psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrid + Mflat;
+      psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+      psfStars->wgtlist[Nref] = 1;
+      psfStars->msklist[Nref] = 0;
+      if (fabs(psfStars->flxlist[Nref]) > 0.03) {
         // fprintf (stderr, "deviant\n");
       }
 
+      float MrelKron = getMrel  (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP);
+      float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON);
+
+      if (isfinite(MrelKron) && isfinite(MsysKron)) {
+        kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrid + Mflat;
+        kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+        kronStars->wgtlist[Nkron] = 1;
+        kronStars->msklist[Nkron] = 0;
+        Nkron ++;
+      }
+
       if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
-        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat, refStars->flxlist[Nref]);
+        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars->flxlist[Nref]);
       }
 
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-        brightStars->flxlist[Nbright] = refStars->flxlist[Nref];
-        brightStars->errlist[Nbright] = refStars->errlist[Nref];
+        brightStars->flxlist[Nbright] = psfStars->flxlist[Nref];
+        brightStars->errlist[Nbright] = psfStars->errlist[Nref];
         brightStars->wgtlist[Nbright] = 1;
         brightStars->msklist[Nbright] = 0;
@@ -770 +756 @@
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
 # if (BASIC_STATS) 
-    liststats (refStars->flxlist, refStars->errlist, NULL, Nref, &stats);
+    liststats (psfStars->flxlist, psfStars->errlist, NULL, Nref, &stats);
 # else
-    liststats_irls (refStars, Nref, &stats);
+    liststats_irls (psfStars, Nref, &stats);
 # endif
-    image[i].Mcal       = stats.mean;
+    image[i].McalPSF    = stats.mean;
     image[i].dMcal      = stats.error;
     image[i].nFitPhotom = Nref;
-    image[i].Xm         = 100.0*log10(stats.chisq);
+    image[i].McalChiSq  = stats.chisq;
     Ncalibrated ++;
+
+    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
+# if (BASIC_STATS) 
+    liststats (kronStars->flxlist, kronStars->errlist, NULL, Nkron, &stats);
+# else
+    liststats_irls (kronStars, Nref, &stats);
+# endif
+    image[i].McalAPER   = stats.mean;
 
     if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
       for (j = 0; j < Nref; j++) {
-        fprintf (stderr, "%1d, %8d : %6.3f  %6.3f  %6.3f  %d\n", (int) i, (int) image[i].imageID, refStars->flxlist[j], refStars->errlist[j], refStars->wgtlist[j], refStars->msklist[j]);
+        fprintf (stderr, "%1d, %8d : %6.3f  %6.3f  %6.3f  %d\n", (int) i, (int) image[i].imageID, kronStars->flxlist[j], kronStars->errlist[j], kronStars->wgtlist[j], kronStars->msklist[j]);
       }
     }
 
     if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
-      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
+      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal);
     }
 
     if (!mark && VERBOSE_IMAGE) {
-      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
+      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal);
     }
 
     if (PLOTSTUFF) {
-      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal);
-      plot_setMcal (refStars->flxlist, Nref, &stats, CLOUD_TOLERANCE);
+      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal);
+      plot_setMcal (psfStars->flxlist, Nref, &stats, CLOUD_TOLERANCE);
     }
 
@@ -803 +797 @@
     image[i].dMagSys = stats.sigma;
 
-    if (image[i].Mcal < -CLOUD_TOLERANCE) {
-      image[i].Mcal = 0.0;
+    if (image[i].McalPSF < -CLOUD_TOLERANCE) {
+      image[i].McalPSF = 0.0;
     }
 
@@ -813 +807 @@
 
   StatDataSetFree (brightStars, 1);
-  StatDataSetFree (refStars, 1);
+  StatDataSetFree (kronStars, 1);
+  StatDataSetFree (psfStars, 1);
 
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);
@@ -848 +843 @@
     if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
-    mlist[N] = image[i].Mcal;
+    mlist[N] = image[i].McalPSF;
     slist[N] = image[i].dMcal;
     dlist[N] = 1;
@@ -873 +868 @@
     mark = FALSE;
     image[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
-    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
+    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset);
     if (mark) { 
       Nmark ++;
@@ -961 +956 @@
     if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
-    Mlist[Nplot] = image[i].Mcal;
+    Mlist[Nplot] = image[i].McalPSF;
     dlist[Nplot] = image[i].dMcal;
     xlist[Nplot] = image[i].secz;
     minAirmass = MIN (image[i].secz, minAirmass);
     maxAirmass = MAX (image[i].secz, maxAirmass);
-    minMcal = MIN (image[i].Mcal, minMcal);
-    maxMcal = MAX (image[i].Mcal, maxMcal);
+    minMcal = MIN (image[i].McalPSF, minMcal);
+    maxMcal = MAX (image[i].McalPSF, maxMcal);
     mindMcal = MIN (image[i].dMcal, mindMcal);
     maxdMcal = MAX (image[i].dMcal, maxdMcal);
@@ -1051 +1046 @@
       c = ImageToCatalog[i][j];
 
-      Mcal  = getMcal  (m, c, flatcorr, catalog);
+      Mcal  = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) continue;
       Mmos  = getMmos  (m, c);
@@ -1092 +1087 @@
     if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
-    list[n] = pow (10.0, 0.01*image[i].Xm);
+    list[n] = image[i].McalChiSq;
     dlist[n] = 1;
     n++;
@@ -1125 +1120 @@
     if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
-    list[n] = image[i].Mcal;
+    list[n] = image[i].McalPSF;
     dlist[n] = 1;
     n++;
@@ -1170 +1165 @@
   return (stats);
 }
+
+void clearImages (void) {
+  image = NULL;
+}
+
+void checkImages (char *name) {
+
+  off_t i;
+  fprintf (stderr, "--- %s ---\n", name);
+  
+  OhanaMemblock *ref = (OhanaMemblock *) image - 1;
+  fprintf (stderr, "file: %s\n", ref->file);
+
+  fprintf (stderr, "0x%08llx \n", (unsigned long long) image);
+
+  for (i = 0; i < Nimage; i++) {
+    fprintf (stderr, "%s : 0x%08x : %6.3f  %6.3f\n", image[i].name, image[i].flags, image[i].McalPSF, image[i].McalAPER);
+  }
+}

branches/czw_branch/20170908/Ohana/src/relphot/src/ImageSubset.c

@@ -63 +63 @@
   char type[16];
 
-  GET_COLUMN (Mcal,      "MCAL",         float);
+  GET_COLUMN (McalPSF,   "MCAL_PSF",     float);
+  GET_COLUMN (McalAPER,  "MCAL_APER",    float);
   GET_COLUMN (dMcal,     "MCAL_ERR",     float);
   GET_COLUMN (imageID,   "IMAGE_ID",     int);
@@ -79 +80 @@
   ALLOCATE (image, ImageSubset, Nrow);
   for (i = 0; i < Nrow; i++) {
-    image[i].Mcal          = Mcal[i];
+    image[i].McalPSF       = McalPSF[i];
+    image[i].McalAPER      = McalAPER[i];
     image[i].dMcal         = dMcal[i];
     image[i].imageID       = imageID[i];
@@ -93 +95 @@
   fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow);
 
-  free (Mcal);
+  free (McalPSF);
+  free (McalAPER);
   free (dMcal);
   free (imageID);
@@ -136 +139 @@
   gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET");
 
-  gfits_define_bintable_column (&theader, "E", "MCAL",         "zero point offset", "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_PSF",     "zero point offset", "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_APER",    "zero point offset", "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MCAL_ERR",     "zero point error",  "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "IMAGE_ID",     "image ID",          NULL,         1.0, FT_BZERO_INT32);
@@ -151 +155 @@
   gfits_create_table (&theader, &ftable);
 
-  float *Mcal, *dMcal;
+  float *McalPSF, *McalAPER, *dMcal;
   unsigned int *imageID, *map, *flags, *tzero;
   int *tessID, *projID, *skycellID;
@@ -158 +162 @@
 
   // create intermediate storage arrays
-  ALLOCATE (Mcal,      float,          Nimage);
+  ALLOCATE (McalPSF,   float,          Nimage);
+  ALLOCATE (McalAPER,  float,          Nimage);
   ALLOCATE (dMcal,     float,          Nimage);
   ALLOCATE (imageID,   unsigned int,   Nimage);
@@ -172 +177 @@
   // assign the storage arrays
   for (i = 0; i < Nimage; i++) {
-    Mcal[i]      = image[i].Mcal;
+    McalPSF[i]   = image[i].McalPSF;
+    McalAPER[i]  = image[i].McalAPER;
     dMcal[i]     = image[i].dMcal;
     imageID[i]   = image[i].imageID;
@@ -186 +192 @@
 
   // add the columns to the output array
-  gfits_set_bintable_column (&theader, &ftable, "MCAL",         Mcal,      Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "MCAL_PSF",     McalPSF,   Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "MCAL_APER",    McalAPER,  Nimage);
   gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR",     dMcal,     Nimage);
   gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID",     imageID,   Nimage);
@@ -198 +205 @@
   gfits_set_bintable_column (&theader, &ftable, "UBERCAL_DIST", ucdist,    Nimage);
 
-  free (Mcal);
+  free (McalPSF);
+  free (McalAPER);
   free (dMcal);
   free (imageID);

branches/czw_branch/20170908/Ohana/src/relphot/src/MosaicOps.c

@@ -111 +111 @@
     
     /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
-    mosaic[Nmosaic].start  = start;
-    mosaic[Nmosaic].stop   = stop;
-    mosaic[Nmosaic].Mcal   = 0.0; // note : mosaic stores only offsets relative to the original image values
-    mosaic[Nmosaic].dMcal  = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].dMsys  = 0.0;
-    mosaic[Nmosaic].Xm = 0.0;
-    mosaic[Nmosaic].flags  = image[i].flags;
-    mosaic[Nmosaic].secz   = image[i].secz;
-    mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    mosaic[Nmosaic].start     = start;
+    mosaic[Nmosaic].stop      = stop;
+    mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
+    mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[Nmosaic].dMsys     = 0.0;
+    mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
+    mosaic[Nmosaic].flags     = image[i].flags;
+    mosaic[Nmosaic].secz      = image[i].secz;
+    mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
 
     // XXX do we need to do something about flag consistency across a mosaic?
@@ -208 +209 @@
   for (i = 0; i < Nmosaic; i++) {
     /* a new mosaic, define ranges */
-    mosaic[i].start    = startTimesMosaic[i];
-    mosaic[i].stop     = 0;
-    mosaic[i].Mcal     = 0.0;
-    mosaic[i].dMcal    = 0.0;
-    mosaic[i].dMsys    = 0.0;
-    mosaic[i].Xm       = 0.0;
-    mosaic[i].flags    = 0;
-    mosaic[i].secz     = NAN;
-    mosaic[i].photcode = 0;
-    mosaic[i].skipCal  = FALSE;
+    mosaic[i].start     = startTimesMosaic[i];
+    mosaic[i].stop      = 0;
+    mosaic[i].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[i].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
+    mosaic[i].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[i].dMsys     = 0.0;
+    mosaic[i].McalChiSq = 0.0;// NAN or 0.0?
+    mosaic[i].flags     = 0;
+    mosaic[i].secz      = NAN;
+    mosaic[i].photcode  = 0;
+    mosaic[i].skipCal   = FALSE;
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -294 +296 @@
       abort();
     }
-    mosaic[j].stop     = stop;
-    mosaic[j].Mcal     = 0.0;
-    mosaic[j].dMcal    = 0.0;
-    mosaic[j].Xm       = 0.0;
-    mosaic[j].dMsys    = subset[i].flags;
-    mosaic[j].flags    = subset[i].flags;
-    mosaic[j].secz     = subset[i].secz;
-    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
+    mosaic[j].stop      = stop;
+
+    mosaic[j].McalPSF   = 0.0;
+    mosaic[j].McalAPER  = 0.0;
+    mosaic[j].dMcal     = 0.0;
+    mosaic[j].McalChiSq = 0.0;
+
+    mosaic[j].dMsys     = subset[i].flags;
+    mosaic[j].flags     = subset[i].flags;
+    mosaic[j].secz      = subset[i].secz;
+    mosaic[j].photcode  = GetPhotcodeEquivCodebyCode (subset[i].photcode);
   }
 
@@ -377 +382 @@
   // init the mosaic array values
   for (i = 0; i < Nmosaic; i++) {
-    mosaic[i].start    = startTimesMosaic[i];
-    mosaic[i].stop     = 0;
-    mosaic[i].Mcal     = 0.0;
-    mosaic[i].dMcal    = 0.0;
-    mosaic[i].dMsys    = 0.0;
-    mosaic[i].Xm       = 0.0;
-    mosaic[i].flags    = 0;
-    mosaic[i].secz     = NAN;
-    mosaic[i].photcode = 0;
-    mosaic[i].skipCal  = FALSE;
+    mosaic[i].start     = startTimesMosaic[i];
+    mosaic[i].stop      = 0;
+    mosaic[i].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[i].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
+    mosaic[i].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
+    mosaic[i].dMsys     = 0.0;
+    mosaic[i].McalChiSq = 0.0;// NAN or 0.0?
+    mosaic[i].flags     = 0;
+    mosaic[i].secz      = NAN;
+    mosaic[i].photcode  = 0;
+    mosaic[i].skipCal   = FALSE;
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -438 +444 @@
     }
     mosaic[j].stop     = stop;
-    mosaic[j].Mcal     = 0.0;
-    mosaic[j].dMcal    = 0.0;
-    mosaic[j].Xm       = 0.0;
-    mosaic[j].dMsys    = image[i].flags;
-    mosaic[j].flags    = image[i].flags;
-    mosaic[j].secz     = image[i].secz;
-    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    mosaic[j].McalPSF   = 0.0;
+    mosaic[j].McalAPER  = 0.0;
+    mosaic[j].dMcal     = 0.0;
+    mosaic[j].McalChiSq = 0.0;
+    mosaic[j].dMsys     = image[i].flags;
+    mosaic[j].flags     = image[i].flags;
+    mosaic[j].secz      = image[i].secz;
+    mosaic[j].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
   }
   MARKTIME("assign images to mosaic: %f sec\n", dtime);
@@ -594 +601 @@
     mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
 
-    mosaic[i].Mcal   = 0.0;
-    mosaic[i].dMcal  = 0.0;
-    mosaic[i].Xm     = 0.0;
+    mosaic[i].McalPSF   = 0.0; 
+    mosaic[i].McalAPER  = 0.0; 
+    mosaic[i].dMcal     = 0.0; 
+    mosaic[i].dMsys     = 0.0;
+    mosaic[i].McalChiSq = 0.0; 
   }
   return;
@@ -612 +621 @@
   double dS, dX, dY;
   double R, D, Rmin, Rmax, Dmin, Dmax;
-  double Mcal, dMcal, Xm;
+  double McalPSF, McalAPER, dMcal, McalChiSq;
 
   fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
@@ -621 +630 @@
     Dmax = Rmax = -360.0;
     dS = 0.0;
-    Mcal = dMcal = Xm = 0;
+    McalPSF = McalAPER = dMcal = McalChiSq = 0;
     for (j = 0; j < MosaicN_Image[i]; j++) {
       m = MosaicToImage[i][j];
@@ -678 +687 @@
       // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
 
-      Mcal  += image[m].Mcal;
-      dMcal += image[m].dMcal;
-      Xm    += image[m].Xm;
-
-      image[m].Mcal   = 0.0;
-      image[m].dMcal  = NAN;
-      image[m].Xm     = NAN_S_SHORT;
+      McalPSF   += image[m].McalPSF;
+      McalAPER  += image[m].McalAPER;
+      dMcal     += image[m].dMcal;
+      McalChiSq += image[m].McalChiSq;
+
+      image[m].McalPSF   = 0.0;
+      image[m].McalAPER  = 0.0;
+      image[m].dMcal     = NAN;
+      image[m].McalChiSq = NAN;
     }
     dS /= MosaicN_Image[i];
@@ -696 +707 @@
     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
 
-    mosaic[i].Mcal   = Mcal / MosaicN_Image[i];
-    mosaic[i].dMcal  = dMcal / MosaicN_Image[i];
-    mosaic[i].Xm     = Xm / MosaicN_Image[i];
+    mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
+    mosaic[i].McalAPER  = McalAPER / MosaicN_Image[i];
+    mosaic[i].dMcal     = dMcal / MosaicN_Image[i];
+    mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
   }
   if (!USE_GRID) return;
@@ -724 +736 @@
     for (j = 0; j < MosaicN_Image[i]; j++) {
       im = MosaicToImage[i][j];
-      image[im].Mcal += mosaic[i].Mcal;
-      image[im].dMcal = mosaic[i].dMcal;
-      image[im].Xm = mosaic[i].Xm;
+      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);
-    }
-    mosaic[i].Mcal = 0.0;
+      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);
+    }
+    mosaic[i].McalPSF  = 0.0;
+    mosaic[i].McalAPER = 0.0;
   }      
 }
@@ -886 +900 @@
 
   if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
-  value = mosaic[i].Mcal;
+  value = mosaic[i].McalPSF;
   return (value);
 }
@@ -911 +925 @@
   off_t Nmax;
   int PoorImages;
-  double *list;
-  double *dlist;
-  double *Mlist;
-  double *dMlist;
+  double *psfMagList;
+  double *psfErrList;
+  double *brightMagList;
+  double *brightErrList;
+  double *kronMagList;
+  double *kronErrList;
 } SetMmosInfo;
 
@@ -945 +961 @@
 
   if (allocLists) {
-    ALLOCATE (info->list, double, Nmax);
-    ALLOCATE (info->dlist, double, Nmax);
-    ALLOCATE (info->Mlist, double, Nmax);
-    ALLOCATE (info->dMlist, double, Nmax);
+    ALLOCATE (info->psfMagList, double, Nmax);
+    ALLOCATE (info->psfErrList, double, Nmax);
+    ALLOCATE (info->kronMagList, double, Nmax);
+    ALLOCATE (info->kronErrList, double, Nmax);
+    ALLOCATE (info->brightMagList, double, Nmax);
+    ALLOCATE (info->brightErrList, double, Nmax);
   }
 }
 
 void SetMmosInfoFree (SetMmosInfo *info) {
-  free (info->list);
-  free (info->dlist);
-  free (info->Mlist);
-  free (info->dMlist);
+  free (info->psfMagList   );
+  free (info->psfErrList   );
+  free (info->kronMagList  );
+  free (info->kronErrList  );
+  free (info->brightMagList);
+  free (info->brightErrList);
 }
 
@@ -1057 +1077 @@
   liststats_setmode (&stats, "INNER_WTMEAN");
 
-  double *list   = info->list;
-  double *dlist  = info->dlist;
-  double *Mlist  = info->Mlist;
-  double *dMlist = info->dMlist;
+  double *psfMagList    = info->psfMagList;
+  double *psfErrList    = info->psfErrList;
+  double *kronMagList   = info->kronMagList;
+  double *kronErrList   = info->kronErrList;
+  double *brightMagList = info->brightMagList;
+  double *brightErrList = info->brightErrList;
 
   assert (Nmos >= 0);
@@ -1125 +1147 @@
   int N = 0;
   for (j = 0; j < N_onMosaic[Nmos]; j++) {
-    float Msys, Mrel, Mcal, Mgrid, Mflat;
       
     off_t m = MosaicToMeasure[Nmos][j];
     off_t c = MosaicToCatalog[Nmos][j];
       
+    // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF
     if (fout) {
-      Mcal  = getMcal  (m, c, flatcorr, catalog);
-      Mgrid = getMgrid (m, c);
-      Mrel  = getMrel  (catalog, m, c);
-      Mflat = getMflat (m, c, flatcorr, catalog);
+      float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
+      float Mgrid    = getMgrid (m, c);
+      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
+      float Mflat    = getMflat (m, c, flatcorr, catalog);
 
       off_t n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-
-      float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
+      float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
+
+      float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
 
       int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
 
-      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad);
+      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad);
     }
 
@@ -1150 +1172 @@
       continue;
     }
-    Mcal  = getMcal  (m, c, flatcorr, catalog);
+    float Mcal = getMcal  (m, c, MAG_CLASS_PSF);
     if (isnan(Mcal)) {
       info->Ncal++;
       continue;
     }
-    Mgrid = getMgrid (m, c);
+    float Mgrid = getMgrid (m, c);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
       continue;
     }
-    Mrel  = getMrel  (catalog, m, c);
-    if (isnan(Mrel)) {
+
+    // Mrel* is the average magnitude for this star.  For PS1 stacks, we have too much
+    // PSF variability.  We need to calibrate the PSF magnitudes separately from the
+    // Aperture-like magnitues.  (We have an option to use the kron magnitudes or the
+    // other apertures here).  I basically need to do this analysis separately for each
+    // magnitude type
+    
+    float MrelPSF = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
+    if (isnan(MrelPSF)) {
       info->Nrel ++;
       continue;
     }
+    float MrelKron = getMrel  (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP);
       
     // image.Mcal is not supposed to include the flat-field correction, so we need to
@@ -1171 +1201 @@
     // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
 
-    Mflat = getMflat (m, c, flatcorr, catalog);
+    float Mflat = getMflat (m, c, flatcorr, catalog);
 
     off_t n = catalog[c].measureT[m].averef;
-    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-    if (isnan(Msys)) {
+    float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
+    if (isnan(MsysPSF)) {
       info->Nsys++;
       continue;
     }
+    float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON);
 
     PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
@@ -1193 +1224 @@
     assert (Nbright >= 0);
 
-    list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
-    dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+    psfMagList[N]  = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    psfErrList[N]  = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+    kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
+    kronErrList[N] = psfErrList[N];
     if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-      Mlist[Nbright] = list[N];
-      dMlist[Nbright] = dlist[N];
+      brightMagList[Nbright] = psfMagList[N];
+      brightErrList[Nbright] = psfErrList[N];
       Nbright ++;
     }
@@ -1223 +1256 @@
   }
 
-  liststats (list, dlist, NULL, N, &stats);
+  liststats (psfMagList, psfErrList, NULL, N, &stats);
   if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
-  myMosaic[0].Mcal  = stats.mean;
-  myMosaic[0].dMcal = stats.error;
+  // for now, I have no reason to measure these separately for camera-level images
+  myMosaic[0].McalPSF    = stats.mean;
+  myMosaic[0].dMcal      = stats.error;
+  myMosaic[0].McalChiSq  = stats.chisq;
   myMosaic[0].nFitPhotom = N;
-  myMosaic[0].Xm    = 100.0*log10(stats.chisq);
+
+  liststats (kronMagList, kronErrList, NULL, N, &stats);
+  myMosaic[0].McalAPER   = stats.mean;
 
   if (testImage) {
@@ -1237 +1274 @@
   if (PLOTSTUFF) {
     fprintf (stderr, "Mmos: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[MosaicToImage[Nmos][0]].name);
-    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+    plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE);
   }
 
   // bright end scatter
-  liststats (Mlist, dMlist, NULL, Nbright, &stats);
+  liststats (brightMagList, brightErrList, NULL, Nbright, &stats);
   myMosaic[0].dMsys = stats.sigma;
 
-  if (myMosaic[0].Mcal < -CLOUD_TOLERANCE) {
-    myMosaic[0].Mcal = 0.0;
+  if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) {
+    myMosaic[0].McalPSF = 0.0;
   }
 
   if (testImage) {
-    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, pow (10.0, 0.01*myMosaic[0].Xm));
+    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].McalPSF, myMosaic[0].dMsys, myMosaic[0].McalChiSq);
   }
 
@@ -1427 +1464 @@
     imageOffset[i] = 0.0;
 
-    if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
+    if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) {
+      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
     } 
-    if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) {
-      imageOffset[i] = -mosaic[i].Mcal; 
-      // NOTE the negative sign: down below, we are going to add in the negative of Mcal
+    if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) {
+      imageOffset[i] = -mosaic[i].McalPSF; 
+      // NOTE the negative sign: down below, we are going to add in the negative of McalPSF
       // to this image, and the propagated mean values for other images
-      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
+      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
     } 
-    if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
+    if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) {
+      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
     } 
   }
@@ -1540 +1577 @@
     }
 
-    // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal)
+    // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF)
     if (imageOffset[i] > 0.0) continue;
 
@@ -1564 +1601 @@
   // find the images / mosaics with negative clouds and save their offset
   for (i = 0; i < Nmosaic; i++) {
-    mosaic[i].Mcal += imageOffset[i];
+    mosaic[i].McalPSF += imageOffset[i];
   }
 
@@ -1605 +1642 @@
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
 
-    list[n] = mosaic[i].Mcal;
+    list[n] = mosaic[i].McalPSF;
     dlist[n] = 1;
     n++;
@@ -1677 +1714 @@
       c = MosaicToCatalog[i][j];
 
-      Mcal = getMcal  (m, c, flatcorr, catalog);
+      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) continue;
       Mgrid = getMgrid (m, c);
       if (isnan(Mgrid)) continue;
-      Mrel = getMrel  (catalog, m, c);
+      Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
       if (isnan(Mrel)) continue;
       N++;
@@ -1719 +1756 @@
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
 
-    list[n]  = pow (10.0, 0.01*mosaic[i].Xm);
+    list[n]  = mosaic[i].McalChiSq;
     dlist[n] = 1;
     n++;
@@ -1751 +1788 @@
     if (mosaic[i].flags & IMAGE_BAD) continue;
     if (mosaic[i].skipCal) continue;
-    mlist[N] = mosaic[i].Mcal;
+    mlist[N] = mosaic[i].McalPSF;
     slist[N] = mosaic[i].dMcal;
     dlist[N] = 1;
@@ -1781 +1818 @@
       Nscatter ++;
     }
-    if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) {
+    if (fabs(mosaic[i].McalPSF - MedOffset) > MaxOffset) {
       mark = TRUE;
       Noffset ++;
@@ -1858 +1895 @@
 
   for (i = 0; i < Nmosaic; i++) {
-    Mlist[i] = mosaic[i].Mcal;
+    Mlist[i] = mosaic[i].McalPSF;
     dlist[i] = mosaic[i].dMcal;
     xlist[i] = mosaic[i].secz;

branches/czw_branch/20170908/Ohana/src/relphot/src/StarOps.c

@@ -34 +34 @@
 }  
 
-float getMrel (Catalog *catalog, off_t meas, int cat) {
-
-  int Nsec, Nsecfilt, photcode;
-  int ave;
-  float value;
-
-  ave = catalog[cat].measureT[meas].averef;
-  photcode = catalog[cat].measureT[meas].photcode;
+// 
+float getMrel (Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source) {
+
+  int ave = catalog[cat].measureT[meas].averef;
+  int photcode = catalog[cat].measureT[meas].photcode;
 
   int ecode = GetPhotcodeEquivCodebyCode (photcode);
-  Nsec = GetPhotcodeNsec(ecode);
-  Nsecfilt = GetPhotcodeNsecfilt ();
+  int Nsec = GetPhotcodeNsec(ecode);
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int entry = Nsecfilt*ave+Nsec;
+
+  SecFilt *secfilt = &catalog[cat].secfilt[entry];
 
   // is this star OK?
-  if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN);  
- 
-  value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
-  return (value);
+  if (secfilt->flags & STAR_BAD) return (NAN);  
+
+  switch (class) {
+    case MAG_CLASS_PSF:
+      switch (source) {
+        case MAG_SRC_CHP:
+          return secfilt->MpsfChp;
+        case MAG_SRC_WRP:
+          return secfilt->MpsfWrp;
+        case MAG_SRC_STK:
+          return secfilt->MpsfStk;
+        default:
+          return NAN;
+      }
+      break;
+    case MAG_CLASS_KRON:
+      switch (source) {
+        case MAG_SRC_CHP:
+          return secfilt->MkronChp;
+        case MAG_SRC_WRP:
+          return secfilt->MkronWrp;
+        case MAG_SRC_STK:
+          return secfilt->MkronStk;
+        default:
+          return NAN;
+      }
+      break;
+    case MAG_CLASS_APER:
+      switch (source) {
+        case MAG_SRC_CHP:
+          return secfilt->MapChp;
+        case MAG_SRC_WRP:
+          return secfilt->MapWrp;
+        case MAG_SRC_STK:
+          return secfilt->MapStk;
+        default:
+          return NAN;
+      }
+      break;
+    case MAG_CLASS_NONE:
+    case MAG_CLASS_DEV: // DeVaucouleur Model (only for galphot)
+    case MAG_CLASS_EXP: // Exponential Model (only for galphot)
+    case MAG_CLASS_SER: // Sersic Model (only for galphot)
+      return NAN;
+      break;
+  }
+  return NAN; // should not be able to reach here
 }
 
@@ -429 +473 @@
   int i;
   off_t j, k, m;
-  float Mcal, Mmos, Mgrid;
 
   MEAS_BAD = ID_MEAS_NOCAL;
@@ -440 +483 @@
       for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
         if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-        Mcal  = getMcal  (m, i, flatcorr, catalog);
-        if (isnan(Mcal)) continue;
-        Mmos  = getMmos  (m, i);
+        float McalPSF  = getMcal  (m, i, MAG_CLASS_PSF);
+        float McalAPER = getMcal  (m, i, MAG_CLASS_APER);
+        if (isnan(McalPSF)) continue;
+        float Mmos  = getMmos  (m, i);
         if (isnan(Mmos)) continue;
-        Mgrid = getMgrid (m, i);
+        float Mgrid = getMgrid (m, i);
         if (isnan(Mgrid)) continue;
+
+        // Note that this operation is setting measure->McalAPER to image->McalAPER only
+        // for the stacks.  Elsewhere (setMrelCatalog) we are using image->McalPSF for all
+        // types of measurements EXCEPT stacks.
+
+        // we need to use McalAPER for stacks (and only stacks)
+        int useStackAper = !USE_MCAL_PSF_FOR_STACK_APER && isGPC1stack(catalog[i].measureT[m].photcode);
 
         // note that measurements for which the image is not selected will not be modified
@@ -451 +502 @@
 
         // set the output calibration
-        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
+        catalog[i].measure[m].McalPSF  = McalPSF  + Mmos + Mgrid;
+        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Mmos + Mgrid : McalPSF  + Mmos + Mgrid;
 
         if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-          myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag");
+          myAssert (isfinite(catalog[i].measure[m].McalPSF), "oops, broke an ubercal mag");
         }
       }
@@ -501 +553 @@
 }
 
+int dumpMags (FILE *fout, Catalog *catalog, int Ncatalog) {
+
+  int i, n;
+  off_t j;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; (j < catalog[i].Naverage) && (j < 2); j++) {
+      fprintf (fout, "%08x %08x %10.6f %10.6f : \n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D); 
+      for (n = 0; n < 5; n++) {
+        fprintf (fout, "secf  %5d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n", n, catalog[i].secfilt[j*Nsecfilt + n].MpsfChp, catalog[i].secfilt[j*Nsecfilt + n].MkronChp, catalog[i].secfilt[j*Nsecfilt + n].MpsfStk, catalog[i].secfilt[j*Nsecfilt + n].MkronStk, catalog[i].secfilt[j*Nsecfilt + n].MpsfWrp, catalog[i].secfilt[j*Nsecfilt + n].MkronWrp);
+      }
+      int m = catalog[i].averageT[j].measureOffset;
+      for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) {
+        fprintf (fout, "meas %5d %5d | %6.3f %6.3f | %6.3f %6.3f\n", m+n, catalog[i].measureT[m+n].photcode, catalog[i].measureT[m+n].M, catalog[i].measureT[m+n].Mkron, catalog[i].measureT[m+n].McalPSF, catalog[i].measureT[m+n].McalAPER); 
+      }
+    }
+  }
+  return TRUE;
+}
+
 void clean_stars (Catalog *catalog, int Ncatalog) {
 
@@ -513 +587 @@
   if (VERBOSE) fprintf (stderr, "marking poor stars\n");
 
-  /* find Xm median -> ChiSq lim must be > median */
+  /* find Mchisq median -> ChiSq lim must be > median */
   for (i = Ntot = 0; i < Ncatalog; i++) {
     Ntot += catalog[i].Naverage; 
@@ -535 +609 @@
         if (isnan(Mchisq)) continue;
         xlist[Ntot] = Mchisq;
-        slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
         dlist[Ntot] = 1;
         Ntot ++;
@@ -551 +625 @@
     for (i = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
-        dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
         float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
         mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
@@ -578 +652 @@
 // NSIGMA_REJECT (5) sigma of the mean
 
+// this function only operations on the PSF magnitudes
+
 # define NSIGMA_CLIP 3.0
 # define NSIGMA_REJECT 5.0
@@ -644 +720 @@
           
           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-          Mcal  = getMcal  (m, i, flatcorr, catalog);
+          Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
           if (isnan(Mcal)) { Ncal ++; continue; }
           Mmos  = getMmos  (m, i);
@@ -695 +771 @@
           
           // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-          Mcal  = getMcal  (m, i, flatcorr, catalog);
+          Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
           if (isnan(Mcal)) continue;
           Mmos  = getMmos  (m, i);
@@ -776 +852 @@
         int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
         if (ecode != seccode) { continue;}
-        Mcal = getMcal  (m, i, flatcorr, catalog);
+        Mcal = getMcal  (m, i, MAG_CLASS_PSF);
         if (isnan(Mcal)) { continue;}
         Mmos = getMmos  (m, i);
@@ -865 +941 @@
       if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
 
-      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
       if (isnan(dM)) continue;
       list[n] = dM;
@@ -906 +982 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-        dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
         bin = dMrel / 0.00025;
         bin = MAX (0, MIN (NBIN-1, bin));
@@ -946 +1022 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
-        xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
+        xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
         value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
         if (isnan((double)(value))) continue;

branches/czw_branch/20170908/Ohana/src/relphot/src/args.c

@@ -136 +136 @@
     remove_argument (N, &argc, argv);
   }
+  SKIP_PARALLEL_GROUPS = 0;
+  if ((N = get_argument (argc, argv, "-skip-parallel-groups"))) {
+    remove_argument (N, &argc, argv);
+    SKIP_PARALLEL_GROUPS = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   // elements needed for parallel regions / parallel images
@@ -325 +331 @@
     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;
   }
 
@@ -544 +555 @@
   PARALLEL_MANUAL = FALSE;
   PARALLEL_SERIAL = FALSE;
+  SKIP_PARALLEL_GROUPS = 0;
 
   HOST_ID = 0;
@@ -694 +706 @@
     remove_argument (N, &argc, argv);
     KEEP_UBERCAL = FALSE;
+  }
+
+  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;
   }
 

branches/czw_branch/20170908/Ohana/src/relphot/src/bcatalog.c

@@ -153 +153 @@
         if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
           if (!KEEP_UBERCAL) {
-            subcatalog[0].measureT[Nmeasure].Mcal = 0.0;
+            subcatalog[0].measureT[Nmeasure].McalPSF  = 0.0;
+            subcatalog[0].measureT[Nmeasure].McalAPER = 0.0;
             subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
           } 
         } else {
           if (getImageEntry (Nmeasure, Ncat) >= 0) {
-            subcatalog[0].measureT[Nmeasure].Mcal = 0.0;
+            subcatalog[0].measureT[Nmeasure].McalPSF  = 0.0;
+            subcatalog[0].measureT[Nmeasure].McalAPER = 0.0;
           }
         }

branches/czw_branch/20170908/Ohana/src/relphot/src/launch_region_hosts.c

@@ -58 +58 @@
     strextend (&command, "-cloud-limit %f", CLOUD_TOLERANCE);
 
-    if (VERBOSE)             strextend (&command, "-v");
-    if (VERBOSE2)            strextend (&command, "-vv");
-    if (RESET)               strextend (&command, "-reset");
-    if (RESET_ZEROPTS)       strextend (&command, "-reset-zpts");
-    if (!KEEP_UBERCAL)       strextend (&command, "-reset-ubercal");
-    if (DophotSelect)        strextend (&command, "-dophot %d", DophotValue);
-    if (ImagSelect)          strextend (&command, "-instmag %f %f", ImagMin, ImagMax);
-    if (MaxDensityUse)       strextend (&command, "-max-density %f", MaxDensityValue);
-    if (SyntheticPhotometry) strextend (&command, "-synthphot");
-    if (USE_BASIC_CHECK)     strextend (&command, "-basic-image-search");
+    if (VERBOSE)                    strextend (&command, "-v");
+    if (VERBOSE2)                   strextend (&command, "-vv");
+    if (RESET)                      strextend (&command, "-reset");
+    if (RESET_ZEROPTS)              strextend (&command, "-reset-zpts");
+    if (!KEEP_UBERCAL)              strextend (&command, "-reset-ubercal");
+    if (DophotSelect)               strextend (&command, "-dophot %d", DophotValue);
+    if (ImagSelect)                 strextend (&command, "-instmag %f %f", ImagMin, ImagMax);
+    if (MaxDensityUse)              strextend (&command, "-max-density %f", MaxDensityValue);
+    if (SyntheticPhotometry)        strextend (&command, "-synthphot");
+    if (USE_BASIC_CHECK)            strextend (&command, "-basic-image-search");
 
-    if (UPDATE)              strextend (&command, "-update");
-    if (MOSAIC_ZEROPT)       strextend (&command, "-mosaic");
-    if (FREEZE_IMAGES)       strextend (&command, "-imfreeze");
-    if (FREEZE_MOSAICS)      strextend (&command, "-mosfreeze");
-    if (PARALLEL)            strextend (&command, "-parallel");
-    if (PARALLEL_MANUAL)     strextend (&command, "-parallel-manual");
-    if (PARALLEL_SERIAL)     strextend (&command, "-parallel-serial");
+    if (UPDATE)                     strextend (&command, "-update");
+    if (MOSAIC_ZEROPT)              strextend (&command, "-mosaic");
+    if (FREEZE_IMAGES)              strextend (&command, "-imfreeze");
+    if (FREEZE_MOSAICS)             strextend (&command, "-mosfreeze");
+    if (CALIBRATE_STACKS_AND_WARPS) strextend (&command, "-only-stacks-and-warps");
+    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
+
+    if (PARALLEL)                   strextend (&command, "-parallel");
+    if (PARALLEL_MANUAL)            strextend (&command, "-parallel-manual");
+    if (PARALLEL_SERIAL)            strextend (&command, "-parallel-serial");
 
     // XXX deprecate this if we are happy with the new version

branches/czw_branch/20170908/Ohana/src/relphot/src/liststats.c

@@ -194 +194 @@
 # define WEIGHT_THRESHOLD 0.3
 
-int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts);
+int fit_least_squares (double *fit, double *err, double *y, double *dy, double *wgt, double *wt, int Npts);
 double VectorFractionInterpolate (double *values, float fraction, int Npts);
 double weight_cauchy (double x);
@@ -253 +253 @@
     
     double oldValue = value;
-    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->wtvals, Npoints)) {
+    if (!fit_least_squares (&value, NULL, dataset->flxlist, dataset->errlist, dataset->wgtlist, dataset->wtvals, Npoints)) {
       value = oldValue;
       break;
@@ -330 +330 @@
     }
 
-    if (!fit_least_squares (&value, dataset->ysample, dataset->dysample, dataset->wtsample, NULL, Nkeep)) continue;
+    if (!fit_least_squares (&value, NULL, dataset->ysample, dataset->dysample, dataset->wtsample, NULL, Nkeep)) continue;
 
     dataset->bvalue[Nboot] = value;
@@ -342 +342 @@
   stats->error = (Shi - Slo) / 2.0;
 
+  // bootstrap can sometimes yield an excessively-optimistic result for the error.  Do not let
+  // the reported error be smaller than the formal error 
+  double errvalue;
+  if (fit_least_squares (&value, &errvalue, dataset->ykeep, dataset->dykeep, dataset->wtkeep, NULL, Nkeep)) {
+    stats->error = MAX (stats->error, errvalue);
+  }
+
   return TRUE;
 }
 
 // wgt is externally-supplied weight, wt is optional
-int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts) { 
+int fit_least_squares (double *fit, double *err, double *y, double *dy, double *wgt, double *wt, int Npts) { 
 
   int i;
@@ -364 +371 @@
   }
   if (S0 == 0.0) return FALSE;
-  *fit = S1 / S0;
+  *fit = S1  / S0;
+  if (err) { *err = 1.0 / S0; }
   return TRUE;
 }

branches/czw_branch/20170908/Ohana/src/relphot/src/load_images.c

@@ -52 +52 @@
       free (subset);
       subset = image;
+    }
+  }
+
+  // reset image values as needed.  always allow 'few' images to succeed, if possible (new
+  // 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);
+    }
+    if (RESET) {
+      if (RESET_ZEROPTS) {
+        if (!KEEP_UBERCAL || !(subset[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) {
+          subset[i].McalPSF  = 0.0;
+          subset[i].McalAPER = 0.0;
+          subset[i].dMcal    = NAN;
+          subset[i].flags   &= ~ID_IMAGE_PHOTOM_UBERCAL;
+        }
+      }
+      subset[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
+      subset[i].ubercalDist = 1000;
     }
   }

branches/czw_branch/20170908/Ohana/src/relphot/src/plot_scatter.c

@@ -40 +40 @@
 
                 if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-                Mcal = getMcal  (m, i, flatcorr, catalog);
+                Mcal = getMcal  (m, i, MAG_CLASS_PSF);
                 if (isnan(Mcal)) continue;
                 Mmos = getMmos  (m, i);
@@ -47 +47 @@
                 if (isnan(Mgrid)) continue;
 
-                Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
+                Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
                 if (isnan(Mrel)) continue;
 

branches/czw_branch/20170908/Ohana/src/relphot/src/plotstuff.c

@@ -179 +179 @@
 
   KapaInitGraph (graphdata);
-  graphdata[0].style = 2;
+  graphdata[0].style = KAPA_PLOT_POINTS; /* points */
   graphdata[0].ptype = 2;
   graphdata[0].ltype = 0;

branches/czw_branch/20170908/Ohana/src/relphot/src/reload_catalogs.c

@@ -180 +180 @@
   for (i = 0; i < Ngroups; i++) {
     // update only a group of unique machines at a time
+    if (i < SKIP_PARALLEL_GROUPS) continue;
     reload_catalog_parallel_group (&groups[i], sky, imageFile);
   }
@@ -247 +248 @@
     if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); }
     if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
+    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
 
     if (!(STAGES & STAGE_CHIP))  { strextend (&command, "-skip-chip"); }

branches/czw_branch/20170908/Ohana/src/relphot/src/relphot_images.c

@@ -42 +42 @@
   MARKTIME("-- load images: %f sec\n", dtime);
 
+  // checkImages("load images");
+
   /* unlock, if we can (else, unlocked below) */
   if (!UPDATE) dvo_image_unlock (&db); 
 
   // load the flat correction table (if defined)
+  // we no longer need to do this: the flatcorrection is now saved in measure.Mflat
 # if (0)
   char flatcorrfile[256];
@@ -60 +63 @@
     /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
 
+    // checkImages("load catalog");
+
     /* match measurements with images, mosaics */
     initImageBins  (catalog, Ncatalog, TRUE);
@@ -77 +82 @@
 
     setExclusions (catalog, Ncatalog, TRUE);
+    // checkImages("set Exclu");
 
     global_stats (catalog, Ncatalog, flatcorr, 0);
@@ -108 +114 @@
         plot_scatter (catalog, Ncatalog, flatcorr); 
       }
+      if (VERBOSE2) dumpMags (stderr, catalog, Ncatalog);
+      // checkImages("start loop");
+
       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);
       setMgrid (catalog, flatcorr);
@@ -213 +228 @@
         // ensure the db format is updated
         dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
+        gfits_modify (&dbX.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT);
+
         char photcodeFile[1024];
         sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
@@ -257 +274 @@
     // ensure the db format is updated
     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
+    gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT);
   }
   if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {

branches/czw_branch/20170908/Ohana/src/relphot/src/relphot_objects.c

@@ -251 +251 @@
     if (USE_BASIC_CHECK)   { strextend (&command, "-basic-image-search"); }
     if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
+    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
 
     if (!(STAGES & STAGE_CHIP))  { strextend (&command, "-skip-chip"); }

branches/czw_branch/20170908/Ohana/src/relphot/src/relphot_parallel_images.c

@@ -36 +36 @@
   client_logger_message ("loaded images\n");
 
-  // once we have read this table, we should remove it for repeat runs
-  // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection?
-
   makeMosaics (image, Nimage, TRUE);
 
@@ -44 +41 @@
 
   // load the flat correction table (if defined)
+  // we no longer need to do this: the flatcorrection is now saved in measure.Mflat
 # if (0)
   char flatcorrfile[256];
@@ -150 +148 @@
   client_logger_message ("all hosts are done the loops\n");
 
-  /* set Mcal & Mmos for bad images */
-  setMcal  (catalog, TRUE, flatcorr);
-  setMmos  (catalog, TRUE, flatcorr);
+  /* 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);
   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
@@ -171 +170 @@
   freeGridBins (Ncatalog);
   freeImages((char *)image);
+  free (image);
 
   SkyListFree(skylist);

branches/czw_branch/20170908/Ohana/src/relphot/src/relphot_parallel_regions.c

@@ -51 +51 @@
     // ensure the db format is updated
     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
+    gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT);
   }
 
@@ -84 +85 @@
 
   // save the changes to the image parameters
-  dvo_image_save (&db, VERBOSE); // this function modifies the db.ftable.buffer: do not free stored Image table
+  dvo_image_save (&db, VERBOSE);
+
+  // dvo_image_save frees db.ftable.buffer (== image) and replaces it: do not free stored Image table
+  clearImages ();
+
   dvo_image_unlock (&db); 
   MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime);

branches/czw_branch/20170908/Ohana/src/relphot/src/relphot_synthphot_catalog.c

@@ -92 +92 @@
       float *value = (float *) zpts->matrix[Nsec].buffer;
       float ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
-      secfilt[Nsec].M    = measure[measSYN[i]].M + ZP;
-      secfilt[Nsec].dM   = 0.6;
-      secfilt[Nsec].Mchisq = NAN;
-      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
+      secfilt[Nsec].MpsfChp  = measure[measSYN[i]].M + ZP;
+      secfilt[Nsec].dMpsfChp = 0.6;
+      secfilt[Nsec].Mchisq   = NAN;
+      secfilt[Nsec].flags   |= ID_SECF_USE_SYNTH;
     }
   }

branches/czw_branch/20170908/Ohana/src/relphot/src/select_images.c

@@ -251 +251 @@
     image[nimage] = timage[i]; 
     inSubset[i] = TRUE;
-    /* always allow 'few' images to succeed, if possible (new images / detections may have
-     * been added) */
-    if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { 
-      image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
-    }
-    if (RESET) {
-      if (RESET_ZEROPTS) {
-        if (!KEEP_UBERCAL || !(image[nimage].flags & ID_IMAGE_PHOTOM_UBERCAL)) {
-          image[nimage].Mcal = 0.0;
-          image[nimage].dMcal = NAN;
-          image[nimage].flags &= ~ID_IMAGE_PHOTOM_UBERCAL;
-        }
-      }
-      image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
-      image[nimage].ubercalDist = 1000;
-    }
     line_number[nimage] = i;
     nimage ++;

branches/czw_branch/20170908/Ohana/src/relphot/src/setMrelCatalog.c

@@ -2 +2 @@
 
 # if (0)
-# define TEST_OBJ_ID 0x0000000e
-# define TEST_CAT_ID 0x000076ee
+# define TEST_OBJ_ID 0x000149b0 
+# define TEST_CAT_ID 0x00000001
 # else
 # define TEST_OBJ_ID 0
@@ -102 +102 @@
 // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
 // 4) some reference photcode of some kind can be specified as fixed and have a high weight
+
+// Although I calculate McalAPER for exposures, I am only using McalPSF for chips.  Note
+// in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for chips
+// and warps, but not stacks
 
 // set mean of chip measurements (selected by photcode range for now):
@@ -224 +228 @@
       // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
       Mmos = Mgrid = 0;
-      Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
+      Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value
     } else {
-      // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y)
-      // NOTE: getMcal_alt does not include measure.Mflat
-      Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
+      // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat
+      Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
       if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
       Mmos  = getMmos  (meas, cat);
@@ -268 +271 @@
     }
 
-    float Mkron, dMkron;
+    // NOTE: we need to calculate the averge chip Kron on each pass to be able to calibrate the stacks
+    float Mkron = PhotCatTiny (&measureT[k], MAG_CLASS_KRON);
+    float dMkron;
     if (isSetMrelFinal) {
-      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
       dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
-      if (CHECK_VALID_MAG(Mkron, dMkron)) {
-        int Nkron = results->kronData[Nsec].Nlist;
-        results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
-        results->kronData[Nsec].errlist[Nkron] = dMkron;
-        results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
-        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
-        results->kronData[Nsec].measSeq[Nkron] = k;
-        results->kronData[Nsec].msklist[Nkron] = 0;
-        results->kronData[Nsec].Nlist ++;
-      }
+    } else {
+      dMkron = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+    }
+    if (CHECK_VALID_MAG(Mkron, dMkron)) {
+      int Nkron = results->kronData[Nsec].Nlist;
+      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
+      results->kronData[Nsec].errlist[Nkron] = dMkron;
+      results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
+      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
+      results->kronData[Nsec].measSeq[Nkron] = k;
+      results->kronData[Nsec].msklist[Nkron] = 0;
+      results->kronData[Nsec].Nlist ++;
     }
 
@@ -377 +383 @@
         }
       }
-      secfilt[Nsec].M      = Mpsf + ZP;
-      secfilt[Nsec].dM     = 0.6;
-      secfilt[Nsec].Mchisq = 0.0;
-      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
+      secfilt[Nsec].MpsfChp  = Mpsf + ZP;
+      secfilt[Nsec].dMpsfChp = 0.6;
+      secfilt[Nsec].Mchisq   = 0.0;
+      secfilt[Nsec].flags   |= ID_SECF_USE_SYNTH;
       continue;
     } 
@@ -386 +392 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *psfstats = &results->psfstats;
-    int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
-    if (Nranking < Nminmeas) { 
+    int NrankingPSF = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
+    if (NrankingPSF < Nminmeas) { 
       secfilt[Nsec].flags |= ID_OBJ_FEW;
     } else {
-      secfilt[Nsec].M    = psfstats->mean;
-      secfilt[Nsec].dM   = psfstats->error;
-      secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
-    }
-    int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
+      secfilt[Nsec].MpsfChp  = psfstats->mean;
+      secfilt[Nsec].dMpsfChp = psfstats->error;
+      secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
+    }
 
     // when running -averages, we have no information about the images, so we cannot set this
@@ -401 +406 @@
     }
 
+    StatType *kronstats = &results->kronstats;
+    int NrankingKRON = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
+    if (NrankingKRON) {
+      secfilt[Nsec].MkronChp  = kronstats->mean; 
+      secfilt[Nsec].dMkronChp = kronstats->error; 
+      secfilt[Nsec].sMkronChp = kronstats->sigma; 
+      secfilt[Nsec].NusedKron = NrankingKRON;
+    }
+
     if (isSetMrelFinal) {
       if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
@@ -406 +420 @@
       }
 
-      // mark the measurements matching this ranking 
-      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
-
-      if (Nranking) {
-        secfilt[Nsec].Mstdev = psfstats->sigma;
-        secfilt[Nsec].Nused  = psfstats->Nmeas;
-        secfilt[Nsec].Mmax   = psfstats->max;
-        secfilt[Nsec].Mmin   = psfstats->min;
+      if (NrankingPSF) {
+        secfilt[Nsec].sMpsfChp = psfstats->sigma;
+        secfilt[Nsec].Nused    = psfstats->Nmeas;
+        secfilt[Nsec].Mmax     = psfstats->max;
+        secfilt[Nsec].Mmin     = psfstats->min;
       }
       secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
@@ -421 +432 @@
       MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax);
 
+      // mark the measurements matching this ranking 
+      int minRankPSF = (NrankingPSF > 0) ? results->psfData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
+
+      int minRankKron = (NrankingKRON > 0) ? results->kronData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
+
       StatType *apstats = &results->apstats;
-      Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
-      if (Nranking) {
-        secfilt[Nsec].Map     = apstats->mean; 
-        secfilt[Nsec].dMap    = apstats->error; 
-        secfilt[Nsec].sMap    = apstats->sigma; 
-        secfilt[Nsec].NusedAp = Nranking;
-      }
-      int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
+      int NrankingAPER = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
+      if (NrankingAPER) {
+        secfilt[Nsec].MapChp  = apstats->mean; 
+        secfilt[Nsec].dMapChp = apstats->error; 
+        secfilt[Nsec].sMapChp = apstats->sigma; 
+        secfilt[Nsec].NusedAp = NrankingAPER;
+      }
+      int minRankAper = (NrankingAPER > 0) ? results->aperData[Nsec].ranking[0] : 10;
       markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
-
-      StatType *kronstats = &results->kronstats;
-      Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
-      if (Nranking) {
-        secfilt[Nsec].Mkron     = kronstats->mean; 
-        secfilt[Nsec].dMkron    = kronstats->error; 
-        secfilt[Nsec].sMkron    = kronstats->sigma; 
-        secfilt[Nsec].NusedKron = Nranking;
-      }
-      int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
-      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
 
       // does this object appear extended in > 50% of measurements?
@@ -518 +525 @@
 // 2) select the BEST detections per filter (regardless of PRIMARY)
 // 3) apply the zero point and AB->Jy transformations
+
+// I calculate McalAPER and McalPSF independently for stacks.  I use McalAPER for Mkron
+// and Map, and McalPSF for Mpsf.  Note in StarOps.c:setMcalOutput I am setting
+// measure->McalAPER to image->McalPSF for chips and warps, but not stacks
+
 int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) {
 
@@ -528 +540 @@
   SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];
 
-  off_t k;
-
-  float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0;
+  float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0;
 
   // set the primary projection cell and skycell for this coordinate
@@ -546 +556 @@
   int haveStackObject = FALSE;
 
-  int Ns;
-  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+  for (int Ns = 0; Ns < Nphotcodes; Ns++) {
 
     int thisCode = photcodes[Ns][0].code;
@@ -556 +565 @@
     int haveStack = FALSE;
 
-    float psfQFbest = 0.0;
-
-    off_t stackBestMeasure = -1;
-    off_t stackPrimaryMeasure = -1;
-
     int isBad = FALSE;
     int isSuspect = FALSE;
 
-    int Nstack = 0; // number for this photcode
-    int NstackDet = 0; // number for this photcode
-
+    int NstackMeas = 0; // number of stack measurements for this photcode
+    int NstackDet  = 0; // number of stack detections for this photcode (not forced)
+
+    int   primaryLevelMax = 0;
+    off_t primaryEntryMax = -1;
+    float primaryValueMax = 0.0;
+
+    int   bestLevelMax = 0;
+    off_t bestEntryMax = -1;
+    float bestValueMax = 0.0;
+
+    int Nprimary = 0;
+
+    // reset all stack-related values for this secfilt:
     secfilt[Nsec].stackBestOff = -1;
     secfilt[Nsec].stackPrmryOff = -1;
-
-    off_t meas = measureOffset;
-    for (k = 0; k < Nmeasure; k++, meas++) {
+    secfilt[Nsec].FpsfStk   = NAN;
+    secfilt[Nsec].dFpsfStk  = NAN;
+    secfilt[Nsec].FkronStk  = NAN;
+    secfilt[Nsec].dFkronStk = NAN;
+    secfilt[Nsec].FapStk    = NAN;
+    secfilt[Nsec].dFapStk   = NAN;
+    secfilt[Nsec].MpsfStk   = NAN;
+    secfilt[Nsec].MkronStk  = NAN;
+    secfilt[Nsec].MapStk    = NAN;
+    secfilt[Nsec].Nstack    = 0;
+    secfilt[Nsec].NstackDet = 0;
+    secfilt[Nsec].flags    &= ~ID_SECF_STACK_FLAGS; // reset the stack flags
+
+    for (off_t k = 0; k < Nmeasure; k++) {
 
       // only examine gpc1 stack data
@@ -581 +607 @@
       if (code->equiv != thisCode) { continue; }
 
-      Nstack ++;
-      if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++;
+      NstackMeas ++;
+      if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++; // 0x04 = PM_SOURCE_MODE2_MATCHED
 
       // clear this bit for all measurements
@@ -590 +616 @@
       haveStack = TRUE;
       haveStackObject = TRUE;
+
+      // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
+
+      int isPrimary = FALSE;
       
-      // ** find the PRIMARY measurement
-
-      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the 
-      // primary skycell for this position
-      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
-        stackPrimaryMeasure = k;
+      // ** is this a PRIMARY measurement? (there may be more than one)
+
+      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement
+      // is from the primary skycell for this position. (note that MatchImageSkycellID
+      // requires the entry in the full Measure table, not just this object)
+      if (MatchImageSkycellID (measureOffset + k, cat, tessID, projectID, skycellID)) {
         measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
-        secfilt[Nsec].stackPrmryOff = meas;
-        myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
-      }
-
-      // ** now choose the BEST measurements (may also be PRIMARY)
-
-      // ensure that we at least have a single best measure 
-      if (stackBestMeasure == -1) stackBestMeasure = k;
-
-      // choose the best psfQFperf value for the BEST measurement
-      if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > psfQFbest)) {
-        psfQFbest = measure[k].psfQFperf;
-        stackBestMeasure = k;
-      }
-      // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary.
-      if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.98)) {
-        psfQFbest = 1000; // force this to be the best entry
-        stackBestMeasure = k;
-      }
-
-      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
-
+        isPrimary = TRUE;
+        Nprimary ++;
+      }
+
+      // soften the error floor (can dM be 0.0?)
+      // XXX EAM : this was int SNvalue -- isfinite(SNvalue) would always be finite
+      float SNvalue = isfinite(measure[k].dM) ? 1.0 / hypot (measure[k].dM, MIN_ERROR) : NAN;
+      int psfQFperfAboveLimit = isfinite(SNvalue) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.95);
+
+      // ** determine the BEST level
+      int bestLevel = 0;
+      if (!isPrimary && !psfQFperfAboveLimit) bestLevel = 1;
+      if ( isPrimary && !psfQFperfAboveLimit) bestLevel = 2;
+      if (!isPrimary &&  psfQFperfAboveLimit) bestLevel = 3;
+      if ( isPrimary &&  psfQFperfAboveLimit) bestLevel = 4;
+
+      // here is the rule for choosing the best value:
+      float bestValue = (bestLevel > 2) ? SNvalue : measure[k].psfQFperf;
+
+      // if we have not reached this level before, set the new level
+      if (bestLevelMax < bestLevel) {
+        bestValueMax = bestValue;
+        bestLevelMax = bestLevel;
+        bestEntryMax = k;
+      }
+      // if we have reached this level before, set the new value if we beat the old one
+      if ((bestLevelMax == bestLevel) && (bestValueMax < bestValue)) {
+        bestValueMax = bestValue;
+        bestEntryMax = k;
+      }
+      myAssert (bestEntryMax > -1, "this should not happen");
+
+      // ** determine the PRIMARY level
+      int primaryLevel = 0;
+      if ( isPrimary && !psfQFperfAboveLimit) primaryLevel = 1;
+      if ( isPrimary &&  psfQFperfAboveLimit) primaryLevel = 2;
+
+      // here is the rule for choosing the PRIMARY value:
+      float primaryValue = 0.0;
+      if (primaryLevel == 1) { primaryValue = measure[k].psfQFperf; }
+      if (primaryLevel == 2) { primaryValue = SNvalue; }
+
+      // if we have not reached this level before, set the new level
+      if (isPrimary && (primaryLevelMax < primaryLevel)) {
+        primaryValueMax = primaryValue;
+        primaryLevelMax = primaryLevel;
+        primaryEntryMax = k;
+      }
+      // if we have reached this level before, set the new value if we beat the old one
+      if (isPrimary && (primaryLevelMax == primaryLevel) && (primaryValueMax < primaryValue)) {
+        primaryValueMax = primaryValue;
+        primaryEntryMax = k;
+      }
+    }
+
+    // if we do not have a stack measurement for this photcode, skip everything below
+    if (!haveStack) continue;
+  
+    // now we have a BEST and a PRIMARY entry (may be the same entry)
+
+    off_t meas = bestEntryMax;
+    off_t measSeq = meas + measureOffset;
+
+    myAssert (meas >= 0, "this should not happen");
+
+    // measurements without an image are either external reference photometry or
+    // data for which the associated image has not been loaded (probably because of
+    // overlaps).  we only want measurements associated with stack images in this loop
+    
+    // match measurement to its image (this is just a check, right?)
+    if (getImageEntry (measSeq, cat) < 0) {
       // measurements without an image are either external reference photometry or
       // data for which the associated image has not been loaded (probably because of
-      // overlaps).  we only want measurements associated with stack images in this loop
-
-      // match measurement to its image (this is just a check, right?)
-      if (getImageEntry (meas, cat) < 0) {
-        // measurements without an image are either external reference photometry or
-        // data for which the associated image has not been loaded (probably because of
-        // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
-        Mmos = Mgrid = 0;
-        Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
-      } else {
-        Mcal  = getMcal_alt  (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd);
-        if (isnan(Mcal))  SKIP_THIS_MEAS_STACK(Ncal);
-        Mmos  = getMmos  (meas, cat);
-        if (isnan(Mmos))  SKIP_THIS_MEAS_STACK(Nmos);
-        Mgrid = getMgrid (meas, cat);
-        if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
-      }
-
-      // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
-      Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
-      if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);
-
-      // data quality assessment
-      isBad |= (measure[k].photFlags & code->photomBadMask);
-      isBad |= (measure[k].psfQF < 0.85);
-      isBad |= isnan(measure[k].psfQF);
-      isBad |= measure[k].dM > 0.2; // S/N < 5.0
-
-      isSuspect |= (measure[k].photFlags & code->photomPoorMask);
-      isSuspect |= (measure[k].psfQFperf < 0.85);
-    }
-
-    if (!haveStack) continue;
-
-    if (!isSuspect && !isBad) {
-      NstackGood ++;
-    }
-    if (isSuspect && !isBad) {
-      NstackSuspect ++;
-    }
-
-    // measurements which are bad will not have a valid stack entry and are skipped
-    k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure;
-    if (k < 0) continue;
-
-    // we are now populating stackDetectID not stack Image ID in secfilt
-    // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
-    // ID = measure[k].extID; // for the stack, this is the stackDetectID
-
-    // get the zero point for the selected image
-    float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid);
+      // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
+      McalPSF  = measure[meas].McalPSF; // check that this is zero for loaded REF value
+      McalAPER = McalPSF; // check that this is zero for loaded REF value
+      Mmos     = 0.0;
+      Mgrid    = 0.0;
+    } else {
+      McalPSF   = getMcal  (measSeq, cat, MAG_CLASS_PSF);
+      McalAPER  = USE_MCAL_PSF_FOR_STACK_APER ? getMcal  (measSeq, cat, MAG_CLASS_PSF) : getMcal  (measSeq, cat, MAG_CLASS_KRON);
+      Mmos      = getMmos  (measSeq, cat);
+      Mgrid     = getMgrid (measSeq, cat);
+    }
+    
+    // ** Here is the math to relate mag,zp to flux
 
     // flux_cgs : erg sec^1 cm^-2 Hz^-1
@@ -703 +744 @@
     // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
 
+    // get the zero point for the selected image
+    // Use a different zero point for the PSF-like and APERTURE-like magnitudes
+    float zpPSF  = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF  + Mmos + Mgrid);
+    float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid);
+
     // zpFactor to go from instrumental flux to Janskies
-    float zpFactor = pow(10.0, -0.4*zp + 3.56);
+    float zpFactorPSF  = pow(10.0, -0.4*zpPSF  + 3.56);
+    float zpFactorAPER = pow(10.0, -0.4*zpAPER + 3.56);
 
     // need to put in AB mag factor to get to Janskies (or uJy?)
-    secfilt[Nsec].FpsfStk   = zpFactor * measure[k].FluxPSF;  
-    secfilt[Nsec].dFpsfStk  = zpFactor * measure[k].dFluxPSF; 
-    secfilt[Nsec].FkronStk  = zpFactor * measure[k].FluxKron; 
-    secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron;
-    secfilt[Nsec].FapStk    = zpFactor * measure[k].FluxAp;
+    secfilt[Nsec].FpsfStk   = zpFactorPSF  * measure[meas].FluxPSF;  
+    secfilt[Nsec].dFpsfStk  = zpFactorPSF  * measure[meas].dFluxPSF; 
+    secfilt[Nsec].FkronStk  = zpFactorAPER * measure[meas].FluxKron; 
+    secfilt[Nsec].dFkronStk = zpFactorAPER * measure[meas].dFluxKron;
+    secfilt[Nsec].FapStk    = zpFactorAPER * measure[meas].FluxAp;
 
     // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
     // use PSF flux error instead here:
-    // secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
-    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxPSF;
+    // secfilt[Nsec].dFapStk   = zpFactorAPER * measure[meas].dFluxAp;
+    secfilt[Nsec].dFapStk   = zpFactorAPER * measure[meas].dFluxPSF;
 
     // Jy to AB mags
-    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
-    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
-    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;
-
-    secfilt[Nsec].stackBestOff = k + measureOffset;
+    secfilt[Nsec].MpsfStk   = (measure[meas].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk)  : NAN;
+    secfilt[Nsec].MkronStk  = (measure[meas].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
+    secfilt[Nsec].MapStk    = (measure[meas].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk)   : NAN;
+
+    // record the measurement which gave the best value
+    secfilt[Nsec].stackBestOff = (bestEntryMax == -1) ? -1 : bestEntryMax + measureOffset;
     myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range");
 
-    secfilt[Nsec].Nstack    = Nstack;
+    // record the selected primary measurement
+    secfilt[Nsec].stackPrmryOff = (primaryEntryMax == -1) ? -1 : primaryEntryMax + measureOffset;
+    myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
+
+    // this is the selected measurement used by secfilt[]
+    measure[meas].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
+    if (Nprimary) {
+      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
+      if (Nprimary > 1) {
+        secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY_MULTIPLE;
+      }
+    }
+
+    // stack measurement used for 'best' was a detection (not forced from the other bands)
+    if ((measure[bestEntryMax].photFlags2 & 0x00000004) == 0) {
+      secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET;
+    }
+
+    // stack measurement used for 'primary' was a detection (not forced from the other bands)
+    if ((primaryEntryMax >= 0) && ((measure[primaryEntryMax].photFlags2 & 0x00000004) == 0)) {
+      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET;
+    }
+
+    secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK;
+
+    // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
+    // float Finst = PhotFluxInst (&measure[meas], MAG_CLASS_PSF);
+    
+    // data quality assessment
+    isBad |= (measure[meas].photFlags & photcodes[Ns][0].photomBadMask);
+    isBad |= (measure[meas].psfQF < 0.85);
+    isBad |= isnan(measure[meas].psfQF);
+    isBad |= measure[meas].dM > 0.2; // S/N < 5.0
+
+    isSuspect |= (measure[meas].photFlags & photcodes[Ns][0].photomPoorMask);
+    isSuspect |= (measure[meas].psfQFperf < 0.85);
+
+    if (!isSuspect && !isBad) {
+      NstackGood ++;
+    }
+    if (isSuspect && !isBad) {
+      NstackSuspect ++;
+    }
+
+    secfilt[Nsec].Nstack    = NstackMeas;
     secfilt[Nsec].NstackDet = NstackDet;
 
-    // this is the measurement used by secfilt[]
-    measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
-    if (k == stackPrimaryMeasure) {
-      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
-    }
-
-    // stack measurement used for 'best' was a detection (not forced from the other bands)
-    if ((measure[k].photFlags2 & 0x00000004) == 0) {
-      secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET;
-    }
-    // stack measurement used for 'primary' was a detection (not forced from the other bands)
-    if ((stackPrimaryMeasure >= 0) && ((measure[stackPrimaryMeasure].photFlags2 & 0x00000004) == 0)) {
-      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET;
-    }
-
-    secfilt[Nsec].flags |= ID_SECF_HAS_PS1_STACK;
   } // Nsecfilt loop
 
@@ -755 +831 @@
 
   int PrimaryIsBest = TRUE;
-  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+  for (int Ns = 0; Ns < Nphotcodes; Ns++) {
     if (!(secfilt[Ns].flags & ID_SECF_HAS_PS1_STACK)) continue; // no stack detection in PS1, nothing is best
     if (secfilt[Ns].flags & ID_SECF_STACK_PRIMARY) continue;    // primary stack detection is best
@@ -780 +856 @@
 // * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
 // analysis is done on flux, not mags (as the faintest objects will be nearly insignificant)
+
+// Although I calculate McalAPER for exposures, I am only using McalPSF for warps..
+// Note in StarOps.c:setMcalOutput I am setting measure->McalAPER to image->McalPSF for
+// chips and warps, but not stacks
+
 int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
   OHANA_UNUSED_PARAM(flatcorr);
@@ -857 +938 @@
       // data for which the associated image has not been loaded (probably because of
       // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
-      Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
+      Mcal = measure[k].McalPSF; // check that this is zero for loaded REF value
     } else {
-      // use getMcal not getMcal_alt?
-      Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
-      // Mcal  = getMcal (meas, cat);
+      Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
       if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
     }

branches/czw_branch/20170908/Ohana/src/relphot/src/setMrelFinal.c

@@ -65 +65 @@
           if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
             if (!KEEP_UBERCAL) {
-              catalog[0].measure[m].Mcal = 0.0;
+              catalog[0].measure[m].McalPSF  = 0.0;
+              catalog[0].measure[m].McalAPER = 0.0;
               catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
             } 
           } else {
             if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
-              catalog[0].measure[m].Mcal = 0.0;
+              catalog[0].measure[m].McalPSF  = 0.0;
+              catalog[0].measure[m].McalAPER = 0.0;
             }
           }
@@ -105 +107 @@
   ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure);
   setMeasureRank (catalog);
-  setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.M = <measure.M + measure.Mflat - image.Mcal>
+  setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.MpsfChp = <measure.M + measure.Mflat - image.Mcal>
   setMcalOutput (catalog, 1, flatcorr); // sets measure.Mcal = image.Mcal
 
@@ -226 +228 @@
         off_t Nim = getImageEntry (m, 0);
         if (Nim > -1) {
-          if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
+          if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip;
           if (isnan(getMmos (m, 0))) goto skip;
         }

branches/czw_branch/20170908/Ohana/src/relphot/src/share_image_mags.c

@@ -74 +74 @@
       continue;
     }
-    images[seq].Mcal        = image_mags[i].Mcal;
+    images[seq].McalPSF     = image_mags[i].McalPSF;
+    images[seq].McalAPER    = image_mags[i].McalAPER;
     images[seq].dMcal       = image_mags[i].dMcal;
     images[seq].dMagSys     = image_mags[i].dMagSys;
-    images[seq].Xm          = image_mags[i].Xm;
+    images[seq].McalChiSq   = image_mags[i].McalChiSq;
     images[seq].nFitPhotom  = image_mags[i].nFitPhotom;
     images[seq].flags       = image_mags[i].flags;
@@ -91 +92 @@
 int set_image_mags (ImageMag *image_mags, Image *image) {
 
-  image_mags->Mcal        = image->Mcal;
+  image_mags->McalPSF     = image->McalPSF;
+  image_mags->McalAPER    = image->McalAPER;
   image_mags->dMcal       = image->dMcal;
   image_mags->dMagSys     = image->dMagSys;
-  image_mags->Xm          = image->Xm;
+  image_mags->McalChiSq   = image->McalChiSq;
   image_mags->nFitPhotom  = image->nFitPhotom;
   image_mags->flags       = image->flags;

branches/czw_branch/20170908/Ohana/src/relphot/src/share_mean_mags.c

@@ -111 +111 @@
     if (Nsec < 0) continue;
 
-    catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].M = meanmags[i].M;
+    catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].MpsfChp = meanmags[i].M;
   }
   free (meanmags);
@@ -122 +122 @@
 int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec) {
 
-  meanmags->M  = secfilt->M;
-  meanmags->dM = secfilt->dM;
+  meanmags->M  = secfilt->MpsfChp;
+  meanmags->dM = secfilt->dMpsfChp;
   meanmags->Mchisq = secfilt->Mchisq;
   meanmags->Nsec = Nsec; // key to secfilt entry

branches/czw_branch/20170908/Ohana/src/relphot/src/synthetic_mags.c

@@ -37 +37 @@
   // (r - i > 0.5) : w = r + 0.268 - 0.435 (r-i) - 0.078(r-i)^2
 
-  float Mr = secfilt[NSr].M;
-  float Mi = secfilt[NSi].M;
+  float Mr = secfilt[NSr].MpsfChp;
+  float Mi = secfilt[NSi].MpsfChp;
 
   if (!isfinite(Mr)) return FALSE;