Changeset 36976

Timestamp:

Jul 7, 2014, 10:57:27 AM (12 years ago)

Author:

eugene

Message:

change secfilt.Xm to secfilt.Mchisq; fix zero point for AB mag / Jy conversion; a little more clarity on the rejections in bcatalog; clearly distinguish measurements from single images (eg, stacks and warps) and mosaiced chips; mosaic identification done by photcode not mosaic name; handle mixed calibration of both simple and mosaiced chips

Location:

branches/eam_branches/ipp-20140610/Ohana/src/relphot

Files:

• include/relphot.h (modified) (2 diffs)
• src/BrightCatalog.c (modified) (4 diffs)
• src/ImageOps.c (modified) (14 diffs)
• src/MeanMagIO.c (modified) (6 diffs)
• src/MosaicOps.c (modified) (24 diffs)
• src/StarOps.c (modified) (8 diffs)
• src/args.c (modified) (1 diff)
• src/assign_images.c (modified) (1 diff)
• src/bcatalog.c (modified) (6 diffs)
• src/initialize.c (modified) (4 diffs)
• src/setMrelCatalog.c (modified) (6 diffs)
• src/share_mean_mags.c (modified) (1 diff)

branches/eam_branches/ipp-20140610/Ohana/src/relphot/include/relphot.h

@@ -48 +48 @@
   float dMsys;
   unsigned short nFitPhotom;
-  short Xm;
+  float Mchisq;
   float secz;
   float ubercalDist;
@@ -119 +119 @@
   float M;
   float dM;
-  float Xm;
+  float Mchisq;
   int Nsec;
   unsigned int objID;

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/BrightCatalog.c

@@ -182 +182 @@
  
     // need to create and assign to flat-field correction
-    GET_COLUMN(M,     "MAG",      float);
-    GET_COLUMN(dM,    "MAG_ERR",  float);
-    GET_COLUMN(Xm,    "MAG_CHI",  float);
-    GET_COLUMN(flags, "FLAGS",    int);
-    GET_COLUMN(Ncode, "NCODE",    short);
-    GET_COLUMN(Nused, "NUSED",    short);
-    GET_COLUMN(Mmin,  "MAG_MIN",  float);
-    GET_COLUMN(Mmax,  "MAG_MAX",  float);
+    GET_COLUMN(M,      "MAG",      float);
+    GET_COLUMN(dM,     "MAG_ERR",  float);
+    GET_COLUMN(Mchisq, "MAG_CHI",  float);
+    GET_COLUMN(flags,  "FLAGS",    int);
+    GET_COLUMN(Ncode,  "NCODE",    short);
+    GET_COLUMN(Nused,  "NUSED",    short);
+    GET_COLUMN(Mmin,   "MAG_MIN",  float);
+    GET_COLUMN(Mmax,   "MAG_MAX",  float);
     gfits_free_header (&theader);
     gfits_free_table  (&ftable);
@@ -196 +196 @@
     ALLOCATE (secfilt, SecFilt, Nrow);
     for (i = 0; i < Nrow; i++) {
-      secfilt[i].M     = M[i];         
-      secfilt[i].dM    = dM[i];
-      secfilt[i].Xm    = Xm[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].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];
     }
     fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow);
 
-    free (M    );
-    free (dM   );
-    free (Xm   );
-    free (flags);
-    free (Ncode);
-    free (Nused);
-    free (Mmin );
-    free (Mmax );
+    free (M     );
+    free (dM    );
+    free (Mchisq);
+    free (flags );
+    free (Ncode );
+    free (Nused );
+    free (Mmin  );
+    free (Mmax  );
     catalog->secfilt = secfilt;
     // assert Nsecfilt * Naverage = Nrow?
@@ -447 +447 @@
     float *M        ; ALLOCATE (M      ,  float, Nsec);
     float *dM       ; ALLOCATE (dM     ,  float, Nsec);
-    float *Xm       ; ALLOCATE (Xm     ,  float, Nsec);
+    float *Mchisq   ; ALLOCATE (Mchisq ,  float, Nsec);
     int   *flags    ; ALLOCATE (flags  ,  int,   Nsec);
     short *Ncode    ; ALLOCATE (Ncode  ,  short, Nsec);
@@ -457 +457 @@
     SecFilt *secfilt = catalog->secfilt;
     for (i = 0; i < Nsec; i++) {
-      M     [i]       = secfilt[i]. M      ;
-      dM    [i]       = secfilt[i]. dM     ;
-      Xm    [i]       = secfilt[i]. Xm     ;
-      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].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   ;
     }
 
     // add the columns to the output array
-    gfits_set_bintable_column (&theader, &ftable, "MAG",      M    , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM   , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Xm   , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "FLAGS",    flags, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "NCODE",    Ncode, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "NUSED",    Nused, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_MIN",  Mmin , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_MAX",  Mmax , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG",      M     , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM    , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Mchisq, Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "FLAGS",    flags , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "NCODE",    Ncode , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "NUSED",    Nused , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_MIN",  Mmin  , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_MAX",  Mmax  , Nsec);
 
     free (M      );
     free (dM     );
-    free (Xm     );
+    free (Mchisq );
     free (flags  );
     free (Ncode  );

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/ImageOps.c

@@ -701 +701 @@
   double MaxOffset, MaxScatter, MedOffset;
 
-  if (FREEZE_IMAGES) return;
+  // if (FREEZE_IMAGES) return;
 
   if (VERBOSE) fprintf (stderr, "marking poor images\n");
@@ -712 +712 @@
   for (i = N = 0; i < Nimage; i++) {
     if (image[i].flags & IMAGE_BAD) continue;
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
+
     mlist[N] = image[i].Mcal;
     slist[N] = image[i].dMcal;
@@ -803 +807 @@
   Graphdata graphdata;
 
-  if (FREEZE_IMAGES) return;
+  // if (FREEZE_IMAGES) return;
 
   ALLOCATE (xlist, double, Nimage);
@@ -816 +820 @@
   float mindMcal    = +100.0;
   float maxdMcal    = -100.0;
-  for (i = 0; i < Nimage; i++) {
-    Mlist[i] = image[i].Mcal;
-    dlist[i] = image[i].dMcal;
-    xlist[i] = image[i].secz;
+
+  int Nplot = 0;
+
+  for (i = 0; i < Nimage; i++) {
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
+
+    Mlist[Nplot] = image[i].Mcal;
+    dlist[Nplot] = image[i].dMcal;
+    xlist[Nplot] = image[i].secz;
     minAirmass = MIN (image[i].secz, minAirmass);
     maxAirmass = MAX (image[i].secz, maxAirmass);
@@ -826 +837 @@
     mindMcal = MIN (image[i].dMcal, mindMcal);
     maxdMcal = MAX (image[i].dMcal, maxdMcal);
+    
+    Nplot ++;
   }
 
@@ -858 +871 @@
   for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
   bzero (Mlist, NBIN*sizeof(double));
-  for (i = 0; i < Nimage; i++) {
+
+  for (i = 0; i < Nimage; i++) {
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
     bin = image[i].dMcal / 0.00025;
     bin = MAX (0, MIN (NBIN - 1, bin));
@@ -883 +899 @@
   StatType stats;
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+
+  // we no longer blindly apply FREEZE_IMAGES to all images, only to mosaics
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -891 +909 @@
   for (i = 0; i < Nimage; i++) {
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
 
     N = 0;
@@ -926 +947 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -935 +956 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
 
     list[n] = pow (10.0, 0.01*image[i].Xm);
@@ -956 +980 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -965 +989 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
 
     list[n] = image[i].Mcal;
@@ -986 +1013 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -995 +1022 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (FREEZE_IMAGES && isGPC1cam) continue;
 
     list[n] = image[i].dMcal;

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MeanMagIO.c

@@ -60 +60 @@
   GET_COLUMN (M,         "MAG",          float);
   GET_COLUMN (dM,        "MAG_ERR",      float);
-  GET_COLUMN (Xm,        "MAG_CHISQ",    float);
+  GET_COLUMN (Mchisq,    "MAG_CHISQ",    float);
   GET_COLUMN (Nsec,      "NSEC",         int);
   GET_COLUMN (objID,     "OBJ_ID",       int);
@@ -73 +73 @@
   ALLOCATE (meanmags, MeanMag, Nrow);
   for (i = 0; i < Nrow; i++) {
-    meanmags[i].M              = M    [i];
-    meanmags[i].dM             = dM   [i];
-    meanmags[i].Xm             = Xm   [i];
-    meanmags[i].Nsec           = Nsec [i];
-    meanmags[i].objID          = objID[i];
-    meanmags[i].catID          = catID[i];
+    meanmags[i].M              = M      [i];
+    meanmags[i].dM             = dM     [i];
+    meanmags[i].Mchisq         = Mchisq [i];
+    meanmags[i].Nsec           = Nsec   [i];
+    meanmags[i].objID          = objID  [i];
+    meanmags[i].catID          = catID  [i];
   }
   fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow);
 
-  free (M    );
-  free (dM   );
-  free (Xm   );
-  free (Nsec );
-  free (objID);
-  free (catID);
+  free (M     );
+  free (dM    );
+  free (Mchisq);
+  free (Nsec  );
+  free (objID );
+  free (catID );
 
   *nmeanmags = Nrow;
@@ -125 +125 @@
   gfits_create_table (&theader, &ftable);
 
-  float *M, *dM, *Xm;
+  float *M, *dM, *Mchisq;
   int *Nsec;
   unsigned int *objID, *catID;
@@ -132 +132 @@
   ALLOCATE (M,         float,          Nmeanmags);
   ALLOCATE (dM,        float,          Nmeanmags);
-  ALLOCATE (Xm,        float,          Nmeanmags);
+  ALLOCATE (Mchisq,    float,          Nmeanmags);
   ALLOCATE (Nsec,      int,            Nmeanmags);
   ALLOCATE (objID,     unsigned int,   Nmeanmags);
@@ -139 +139 @@
   // assign the storage arrays
   for (i = 0; i < Nmeanmags; i++) {
-    M    [i]   = meanmags[i].M    ;
-    dM   [i]   = meanmags[i].dM   ;
-    Xm   [i]   = meanmags[i].Xm   ;
-    Nsec [i]   = meanmags[i].Nsec ;
-    objID[i]   = meanmags[i].objID;
-    catID[i]   = meanmags[i].catID;
+    M     [i]   = meanmags[i].M     ;
+    dM    [i]   = meanmags[i].dM    ;
+    Mchisq[i]   = meanmags[i].Mchisq;
+    Nsec  [i]   = meanmags[i].Nsec  ;
+    objID [i]   = meanmags[i].objID ;
+    catID [i]   = meanmags[i].catID ;
   }
 
@@ -150 +150 @@
   gfits_set_bintable_column (&theader, &ftable, "MAG",       M,         Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",   dM,        Nmeanmags);
-  gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm,        Nmeanmags);
+  gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq,    Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "NSEC",      Nsec,      Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "OBJ_ID",    objID,     Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "CAT_ID",    catID,     Nmeanmags);
 
-  free (M    );
-  free (dM   );
-  free (Xm   );
-  free (Nsec );
-  free (objID);
-  free (catID);
+  free (M     );
+  free (dM    );
+  free (Mchisq);
+  free (Nsec  );
+  free (objID );
+  free (catID );
 
   FILE *f = fopen (filename, "w");

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/MosaicOps.c

@@ -96 +96 @@
     
     /* 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].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].Mchisq = 0.0;
     mosaic[Nmosaic].flags  = image[i].flags;
-    mosaic[Nmosaic].secz  = image[i].secz;
+    mosaic[Nmosaic].secz   = image[i].secz;
     mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
 
@@ -172 +172 @@
   // generate a list of all subset image start times
   ALLOCATE (startTimes, unsigned int, Nsubset);
+  int Nmoschip = 0;
   for (i = 0; i < Nsubset; i++) {
-    startTimes[i] = subset[i].tzero;
-  }
-  sort_times (startTimes, Nsubset);
+    int isGPC1cam = ((subset[i].photcode >= 10000) && (subset[i].photcode <= 10500));
+    if (!isGPC1cam) {
+      continue;
+    }
+    startTimes[Nmoschip] = subset[i].tzero;
+    Nmoschip ++;
+  }
+  sort_times (startTimes, Nmoschip);
   MARKTIME("create array of all image obstimes: %f sec\n", dtime);
   
@@ -184 +190 @@
 
   // generate a list of the unique start times (these define the mosaics)
-  for (i = 0; i < Nsubset; i++) {
+  for (i = 0; i < Nmoschip; i++) {
     if (startTimes[i] < startTimesMosaic[Nmosaic]) {
       fprintf (stderr, "error?\n");
@@ -210 +216 @@
   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].start    = startTimesMosaic[i];
+    mosaic[i].stop     = 0;
+    mosaic[i].Mcal     = 0.0;
+    mosaic[i].dMcal    = 0.0;
+    mosaic[i].dMsys    = 0.0;
+    mosaic[i].Mchisq   = 0.0;
+    mosaic[i].flags    = 0;
+    mosaic[i].secz     = NAN;
     mosaic[i].photcode = 0;
-    mosaic[i].skipCal = FALSE;
+    mosaic[i].skipCal  = FALSE;
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -238 +244 @@
     }
     
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    if (!pname) continue;
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
+    // this is gpc1-specific.  needs a more general solution
+    int isGPC1cam = ((image[i].photcode >= 10000) && (image[i].photcode <= 10500));
+    if (!isGPC1cam) {
+      continue;
+    }
+
+    if (0) {
+      // XXX this version fails if we have stack and fforce data: GPC1.g.SkyChip
+      /* select valid mosaic images by photcode */
+      pname = GetPhotcodeNamebyCode (image[i].photcode);
+      if (!pname) continue;
+      status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
+      if (status) continue;
+    }
 
     /* set image time range */
@@ -262 +277 @@
 
   // assign each image to a mosaic
+  int Nsimple = 0;
   for (i = 0; i < Nsubset; i++) {
     ImageToMosaic[i] = -1;
 
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (subset[i].photcode);
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
+    // this is gpc1-specific.  needs a more general solution
+    int isGPC1cam = ((subset[i].photcode >= 10000) && (subset[i].photcode <= 10500));
+    if (!isGPC1cam) {
+      Nsimple ++;
+      continue;
+    }
+
+    if (0) {
+      // XXX this version fails if we have stack and fforce data: GPC1.g.SkyChip
+      /* select valid mosaic images by photcode */
+      pname = GetPhotcodeNamebyCode (subset[i].photcode);
+      status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
+      if (status) continue;
+    }
 
     start = subset[i].tzero;
@@ -302 +328 @@
     mosaic[j].Mcal  = 0.0;
     mosaic[j].dMcal = 0.0;
-    mosaic[j].Xm    = 0.0;
+    mosaic[j].Mchisq = 0.0;
     mosaic[j].dMsys = subset[i].flags;
     mosaic[j].flags = subset[i].flags;
@@ -316 +342 @@
   initMosaicGrid (subset, Nsubset);
 
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic);
+  fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
   return;
 }
@@ -384 +410 @@
     mosaic[i].dMcal = 0.0;
     mosaic[i].dMsys = 0.0;
-    mosaic[i].Xm    = 0.0;
+    mosaic[i].Mchisq    = 0.0;
     mosaic[i].flags = 0;
     mosaic[i].secz  = NAN;
@@ -443 +469 @@
     mosaic[j].Mcal  = 0.0;
     mosaic[j].dMcal = 0.0;
-    mosaic[j].Xm    = 0.0;
+    mosaic[j].Mchisq    = 0.0;
     mosaic[j].dMsys = image[i].flags;
     mosaic[j].flags = image[i].flags;
@@ -566 +592 @@
       Mcal  += image[m].Mcal;
       dMcal += image[m].dMcal;
-      Xm    += image[m].Xm;
+      Xm    += image[m].Mchisq;
 
       // for ubercal images, we (elsewhere) keep Mcal frozen
@@ -573 +599 @@
       image[m].Mcal  = 0.0;
       image[m].dMcal = NAN;
-      image[m].Xm    = NAN_S_SHORT;
+      image[m].Mchisq    = NAN_S_SHORT;
     }
 
@@ -617 +643 @@
     mosaic[i].coords.pc2_1  = 0.0;
 
-    mosaic[i].Mcal  = Mcal / MosaicN_Image[i];
-    mosaic[i].dMcal = dMcal / MosaicN_Image[i];
-    mosaic[i].Xm    = Xm / MosaicN_Image[i];
+    mosaic[i].Mcal   = Mcal / MosaicN_Image[i];
+    mosaic[i].dMcal  = dMcal / MosaicN_Image[i];
+    mosaic[i].Mchisq = Xm / MosaicN_Image[i];
   }
   return;
@@ -669 +695 @@
       Mcal  += image[m].Mcal;
       dMcal += image[m].dMcal;
-      Xm    += image[m].Xm;
+      Xm    += image[m].Mchisq;
 
       // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but 
@@ -679 +705 @@
 
       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-      image[m].Mcal  = 0.0;
-      image[m].dMcal = NAN;
-      image[m].Xm    = NAN_S_SHORT;
+      image[m].Mcal   = 0.0;
+      image[m].dMcal  = NAN;
+      image[m].Mchiqs = NAN;
     }
     dS /= MosaicN_Image[i];
@@ -697 +723 @@
     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].Mcal   = Mcal / MosaicN_Image[i];
+    mosaic[i].dMcal  = dMcal / MosaicN_Image[i];
+    mosaic[i].Mchisq = Xm / MosaicN_Image[i];
   }
   if (!USE_GRID) return;
@@ -729 +755 @@
       image[im].Mcal += mosaic[i].Mcal;
       image[im].dMcal = mosaic[i].dMcal;
-      image[im].Xm = mosaic[i].Xm;
+      image[im].Mchisq = mosaic[i].Mchisq;
       image[im].dMagSys = mosaic[i].dMsys;
       image[im].nFitPhotom = mosaic[i].nFitPhotom;
@@ -840 +866 @@
   mosID = ImageToMosaic[idx];
   if (mosID < 0) {
-    Image *image = getimage(idx);
-    fprintf (stderr, "unmatched image %s\n", image[0].name);
+    // Image *image = getimage(idx);
+    // fprintf (stderr, "unmatched image %s\n", image[0].name);
+    // skip measurements from simple chips (not matched to a mosaic by definition)
     return;
   }
@@ -882 +909 @@
   if (!MOSAIC_ZEROPT) return (0);
 
+  // unassigned measurements belong to simple chips
   i = MeasureToMosaic[cat][meas];
-  if (i == -1) return (NAN);
+  if (i == -1) return (0.0);
 
   if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
@@ -1219 +1247 @@
   myMosaic[0].dMcal = stats.error;
   myMosaic[0].nFitPhotom = N;
-  myMosaic[0].Xm    = 100.0*log10(stats.chisq);
+  myMosaic[0].Mchisq    = stats.chisq;
 
   if (testImage) {
@@ -1239 +1267 @@
 
   if (testImage) {
-    fprintf (stderr, "%f %f  :  %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm));
+    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Mchisq);
   }
 
@@ -1699 +1727 @@
     if (mosaic[i].flags & IMAGE_BAD) continue;
     if (mosaic[i].skipCal) continue;
-    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
+    list[n] = mosaic[i].Mchisq;
     dlist[n] = 1;
     n++;
@@ -1784 +1812 @@
   off_t i, j, m, c, N, Nimage;
   double *xlist, *ylist;
-  // double Xmin, Xmax, Ymin, Ymax;
   char string[64];
   Graphdata graphdata;
@@ -1803 +1830 @@
   for (i = 0; i < Nmosaic; i++) {
     N = 0;
-    // Xmin = Ymin = +360.0;
-    // Xmax = Ymax = -360.0;
     for (j = 0; j < N_onMosaic[i]; j++) {
       

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/StarOps.c

@@ -153 +153 @@
   int i;
 
-  fprintf (stderr, "this version of setMrel is invalid for now\n");
-  exit (1);
+  // fprintf (stderr, "this version of setMrel is invalid for now\n");
+  // exit (1);
         
   int Nsecfilt = GetPhotcodeNsecfilt ();
@@ -673 +673 @@
       liststats (list, dlist, wlist, N, &stats);
 
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M      = stats.mean;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM     = stats.error;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN;
 
       // when running -averages, we have no information about the images, so we cannot set this
@@ -953 +953 @@
 
   int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
-  float dM, Xm;
+  float dM;
   double Chisq, MaxScatter, MaxChisq;
   double *xlist, *slist, *dlist;
@@ -981 +981 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
-        Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        if (Xm == -1) continue;
-        Chisq = pow (10.0, 0.01*Xm);
-        xlist[Ntot] = Chisq;
+        Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        xlist[Ntot] = Mchisq;
         slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
         dlist[Ntot] = 1;
@@ -1002 +1000 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
-        Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        Chisq = pow (10.0, 0.01*Xm);
-        mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
+        Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        mark = (dM > MaxScatter) || (isnan(Mchisq) || (Mchisq > MaxChisq);
         if (mark) {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
           Ndel ++;
-          if (dM > MaxScatter) { Nscat ++; }
-          if (Xm == NAN_S_SHORT) { Nnan ++; }
-          if (Chisq > MaxChisq) { Nchi ++; }
+          if (dM > MaxScatter)   { Nscat ++; }
+          if (isnan(Mchisq))     { Nnan ++; }
+          if (Mchisq > MaxChisq) { Nchi ++; }
         } else {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
@@ -1236 +1233 @@
 
   off_t j, Ntot;
-  int i, n, Xm;
+  int i, n;
   double *list, *dlist;
   StatType stats;
@@ -1257 +1254 @@
       if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
         
-      Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-      if (Xm == NAN_S_SHORT) continue;
-      list[n] = pow (10.0, 0.01*Xm);
+      Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+      if (isnan(Mchisq)) continue;
+      list[n] = Mchisq;
       dlist[n] = 1;
       n++;
@@ -1380 +1377 @@
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
         xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
-        value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        if (value == NAN_S_SHORT) continue;
-        ylist[N] = 0.01*value;
+        value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        if (isnan(value)) continue;
+        ylist[N] = value;
         N++;
       }

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/args.c

@@ -227 +227 @@
   }
 
-  SHOW_PARAMS = FALSE;
+  SHOW_PARAMS = TRUE;
   if ((N = get_argument (argc, argv, "-params"))) {
     remove_argument (N, &argc, argv);
-    SHOW_PARAMS = TRUE;
+    SHOW_PARAMS = FALSE;
   }
 

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/assign_images.c

@@ -27 +27 @@
     makeMosaics (image, Nimage);
 
-    // center coords and Mcal, dMcal, Xm for the mosaics
+    // center coords and Mcal, dMcal, Mchisq for the mosaics
     setMosaicCenters (image, Nimage);
     MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/bcatalog.c

@@ -22 +22 @@
   ALLOCATE (subcatalog[0].secfilt,  SecFilt,     NAVERAGE*Nsecfilt);
   Nmeasure = Naverage = 0;
+
+  int NmMax = 0;
+  int NmMin = catalog[0].Nmeasure;
 
   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
@@ -82 +85 @@
         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
-        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm     = NAN;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN;
         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
         // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
@@ -187 +190 @@
       }
       myAssert (Nmeasure < NMEASURE, "realloc failure");
-    }
+    } // end of catalog.average.Nmeasure loop
 
     // skip object if it is likely to be a galaxy
@@ -212 +215 @@
       continue; 
     }
+    NmMin = MIN (Nm, NmMin);
+    NmMax = MAX (Nm, NmMax);
 
     // for w-band photometry (& other cases?) convert gri photometry
@@ -242 +247 @@
   subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
 
+  if (VERBOSE) {
+    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" ("OFF_T_FMT" measures) for catalog %s, %d < Nm < %d\n", 
+             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
+    fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+             Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
+  }
+
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
@@ -247 +259 @@
   }
 
-  if (VERBOSE) {
-    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", 
-             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
-    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
-             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
-  }
   return (TRUE);
 }

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/initialize.c

@@ -36 +36 @@
 
   if (SHOW_PARAMS) {
-    fprintf (stderr, "current parameter settings:\n");
+    int Ns;
+    fprintf (stderr, "subset selection criteria:\n");
+    fprintf (stderr, "  photcodes ");
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
+      if (Ns == Nphotcodes - 1) {
+        fprintf (stderr, "%s\n", photcodes[Ns][0].name);
+      } else {
+        fprintf (stderr, "%s, ", photcodes[Ns][0].name);
+      }
+    }
     if (TimeSelect) {
       fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", ohana_sec_to_date (TSTART), ohana_sec_to_date (TSTOP));
@@ -42 +51 @@
       fprintf (stderr, "TimeSelect: FALSE\n");
     }
+    if (DophotSelect) {
+      fprintf (stderr, "DophotSelect: TRUE (%d)\n", DophotValue);
+    } else {
+      fprintf (stderr, "DophotSelect: FALSE\n");
+    }
+    fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n");
+
+    // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n");
+
+    fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM);
+    fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax);
+
+    fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW);
+
     fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
     fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", 
@@ -49 +72 @@
              (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
 
-    fprintf (stderr, "MAG_LIM                %lf\n", MAG_LIM);
     fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);
     fprintf (stderr, "IMAGE_SCATTER          %lf\n", IMAGE_SCATTER);
@@ -56 +78 @@
     fprintf (stderr, "GSCFILE                %s\n",  GSCFILE);
     fprintf (stderr, "CATDIR                 %s\n",  CATDIR);
-    exit (0);
   }
 

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/setMrelCatalog.c

@@ -354 +354 @@
     liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
 
-    secfilt[Nsec].M  = psfstats->mean;
-    secfilt[Nsec].dM = psfstats->error;
-    secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT;
+    secfilt[Nsec].M      = psfstats->mean;
+    secfilt[Nsec].dM     = 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
@@ -613 +613 @@
     secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
 
-    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? zp -2.5*log10(measure[k].FluxPSF) : NAN;
-    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? zp -2.5*log10(measure[k].FluxKron) : NAN;
-    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? zp -2.5*log10(measure[k].FluxAp) : NAN;
+    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].stackDetectID = ID;
@@ -716 +716 @@
         if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
       }
-      float Fcal = MagToFlux(Mcal);
+      float Fcal = MagToFlux(-Mcal);
 
       // in the calculations below, 
@@ -767 +767 @@
     secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
     secfilt[Nsec].NusedWrp = psfstats->Nmeas;
-    secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : -2.5*log(secfilt[Nsec].FpsfWrp); // XXX need to add zero-point
+    secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
 
     // NOTE : use the modified weight for apmags as well as psf mags
@@ -775 +775 @@
     secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN; 
     secfilt[Nsec].NusedApWrp = Nap;
-    secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : -2.5*log(secfilt[Nsec].FapWrp); // XXX need to add zero-point
+    secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
 
     liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
@@ -782 +782 @@
     secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN; 
     secfilt[Nsec].NusedKronWrp = Nkron;
-    secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : -2.5*log(secfilt[Nsec].FkronWrp); // XXX need to add zero-point
+    secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
   }
   return (TRUE);

branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/share_mean_mags.c

@@ -123 +123 @@
   meanmags->M  = secfilt->M;
   meanmags->dM = secfilt->dM;
-  meanmags->Xm = secfilt->Xm;
+  meanmags->Mchisq = secfilt->Mchisq;
   meanmags->Nsec = Nsec; // key to secfilt entry