Changeset 33157

Timestamp:

Jan 26, 2012, 5:21:00 AM (14 years ago)

Author:

eugene

Message:

relphot runs with ubercal & flatcorr, signs seem right

Location:

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

Files:

• include/relphot.h (modified) (3 diffs)
• src/ImageOps.c (modified) (11 diffs)
• src/MosaicOps.c (modified) (8 diffs)
• src/StarOps.c (modified) (5 diffs)
• src/args.c (modified) (1 diff)
• src/help.c (modified) (1 diff)
• src/plot_scatter.c (modified) (2 diffs)
• src/plotstuff.c (modified) (1 diff)
• src/relphot.c (modified) (2 diffs)

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

@@ -79 +79 @@
 int    FREEZE_MOSAICS;
 int    USE_GRID;
+int    KEEP_UBERCAL;
 char  *OUTROOT;
 int    PLOTDELAY;
@@ -157 +158 @@
 off_t         getImageEntry       PROTO((off_t meas, int cat));
 float         getMcal             PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
+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));
@@ -196 +198 @@
 int           reload_images       PROTO((FITS_DB *db));
 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
-void          setMcal             PROTO((Catalog *catalog, int Poor));
+void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
 void          setMcalFinal        PROTO((void));
 int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));

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

@@ -1 +1 @@
 # include "relphot.h"
+void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 
 // we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage)
@@ -253 +254 @@
 }
 
+// returns image.Mcal - ff(x,y)
 float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
 
@@ -268 +270 @@
   int flat_id = image[i].photom_map_id;
   if (flat_id) {
-    offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measure[meas].Xccd, catalog[cat].measure[meas].Yccd);
-  }
-  value += offset;
+    offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
+  }
+  value -= offset;
 
   return (value);
+}
+
+float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
+
+  off_t i = MeasureToImage[cat][meas];
+  if (i == -1) return (NAN);
+
+  float offset = 0.0;
+
+  // to do this, I need to pass in the catalog and flatcorr pointers
+  int flat_id = image[i].photom_map_id;
+  if (flat_id) {
+    offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
+  }
+
+  return (offset);
 }
 
@@ -285 +303 @@
 
 /* determine Mcal values for all images */
-void setMcal (Catalog *catalog, int PoorImages) {
+void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
   off_t i, j, m, c, n, N, Nmax, mark, bad;
   int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright;
-  float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2;
+  float Msys, Mrel, Mmos, Mgrid, Mflat, McalBright, McalBright2;
   double *list, *dlist;
   StatType stats;
@@ -321 +339 @@
       if (!bad) continue;
     }      
+
+    /* we optionally do not recalibrate images with UBERCAL zero points */
+    if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
 
     N = 0;
@@ -348 +369 @@
       }
       
+      // image.Mcal is not supposed to include the flat-field correction, so we need to
+      // apply that offset as well here for this image (in other words, each detection is
+      // being compared to the model, excluding the zero point, Mcal.  The model includes
+      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
+
+      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]);
@@ -354 +382 @@
         continue;
       }
-      list[N] = Msys - Mrel - Mmos - Mgrid;
+
+      fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %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);
+
+      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
@@ -385 +416 @@
     image[i].Xm    = 100.0*log10(stats.chisq);
 
+    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+
     if (image[i].Mcal < -CLOUD_TOLERANCE) {
       image[i].Mcal = 0.0;
@@ -454 +487 @@
 }
 
+void plot_setMcal (double *list, int Npts, StatType *stats, float clouds) {
+
+  off_t i;
+  double *xlist;
+  Graphdata graphdata;
+
+  if (!PLOTSTUFF) return;
+
+  ALLOCATE (xlist, double, Npts);
+
+  /**** Mcal vs seq ****/
+  float minMcal    = +100.0;
+  float maxMcal    = -100.0;
+  for (i = 0; i < Npts; i++) {
+    xlist[i] = i;
+    minMcal = MIN (list[i], minMcal);
+    maxMcal = MAX (list[i], maxMcal);
+  }
+
+  float McalRange  = MAX(1.2*(maxMcal - minMcal), 0.21);
+  float McalCenter = 0.5*(maxMcal + minMcal);
+
+  plot_defaults (&graphdata);
+  graphdata.xmin = -1;
+  graphdata.xmax = Npts + 1;
+  graphdata.ymin = McalCenter - 0.5*McalRange;
+  graphdata.ymax = McalCenter + 0.5*McalRange;
+  plot_list (&graphdata, xlist, list, Npts, "sequence vs Mcal", "%s.seq.png", OUTROOT);
+
+  free (xlist);
+}
+
 void plot_images () {
 
@@ -467 +532 @@
 
   /**** dMcal vs airmass ****/
+  float minAirmass = 1000.0;
+  float maxAirmass =    0.0;
+  float minMcal    = +100.0;
+  float maxMcal    = -100.0;
+  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;
-  }
+    minAirmass = MIN (image[i].secz, minAirmass);
+    maxAirmass = MAX (image[i].secz, maxAirmass);
+    minMcal = MIN (image[i].Mcal, minMcal);
+    maxMcal = MAX (image[i].Mcal, maxMcal);
+    mindMcal = MIN (image[i].dMcal, mindMcal);
+    maxdMcal = MAX (image[i].dMcal, maxdMcal);
+  }
+
+  float AirmassRange  = MAX(1.2*(maxAirmass - minAirmass), 0.25);
+  float AirmassCenter = 0.5*(maxAirmass + minAirmass);
+
+  float McalRange  = MAX(1.2*(maxMcal - minMcal), 0.21);
+  float McalCenter = 0.5*(maxMcal + minMcal);
+
+  float dMcalRange  = MAX(1.2*(maxdMcal - mindMcal), 0.21);
+  float dMcalCenter = 0.5*(maxdMcal + mindMcal);
 
   plot_defaults (&graphdata);
-  graphdata.ymin = PlotdMmin; 
-  graphdata.ymax = PlotdMmax;
+  graphdata.xmin = AirmassCenter - 0.5*AirmassRange;
+  graphdata.xmax = AirmassCenter + 0.5*AirmassRange;
+  graphdata.ymin = McalCenter - 0.5*McalRange;
+  graphdata.ymax = McalCenter + 0.5*McalRange;
   plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
+
   plot_defaults (&graphdata);
+  graphdata.xmin = McalCenter - 0.5*McalRange;
+  graphdata.xmax = McalCenter + 0.5*McalRange;
+  graphdata.ymin = dMcalCenter - 0.5*dMcalRange;
+  graphdata.ymax = dMcalCenter + 0.5*dMcalRange;
   plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", "%s.Mcal.dMcal.png", OUTROOT);
 
@@ -495 +588 @@
   plot_defaults (&graphdata);
   graphdata.style = 1;
+  graphdata.xmin = dMcalCenter - 0.5*dMcalRange;
+  graphdata.xmax = dMcalCenter + 0.5*dMcalRange;
   plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT);
 

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

@@ -102 +102 @@
     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?
 
     /* add image to mosaic image list */
@@ -175 +177 @@
       dMcal += image[m].dMcal;
       Xm    += image[m].Xm;
+
+      // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but 
+      // solve for a single offset for each exposure (Mosaic.Mcal).  
+
+      // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
+
+      // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
+
       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
       image[m].Mcal  = 0.0;
@@ -215 +225 @@
   image = getimages (&Nimage, NULL);
 
+  // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL
+  // images? (no need to update them)
+
   // copy the mosaic results to the images.  set the mosaic Mcal to 0.0 since we have moved its
   // impact to the images
@@ -380 +393 @@
   off_t i, j, m, c, N, Nmax;
   int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt;
-  float Msys, Mrel, Mcal, Mgrid;
+  float Msys, Mrel, Mcal, Mgrid, Mflat;
   double *list, *dlist, *Mlist, *dMlist;
   StatType stats;
@@ -414 +427 @@
       if (!bad) continue;
     }      
+
+    /* we optionally do not recalibrate images with UBERCAL zero points */
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
 
     N = 0;
@@ -441 +457 @@
       }
       
+      // image.Mcal is not supposed to include the flat-field correction, so we need to
+      // apply that offset as well here for this image (in other words, each detection is
+      // being compared to the model, excluding the zero point, Mcal.  The model includes
+      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
+
+      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]);
@@ -447 +470 @@
         continue;
       }
-      list[N]  = Msys - Mrel - Mcal - Mgrid;
+      list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       Mlist[N] = Msys;
@@ -499 +522 @@
 // to 0.0.  At the same time, we make a guess to the effective impact on all other images,
 // driven by the coupling of common stars. 
+// XXX this function is no longer used because we force significantly negative clouds to 0.0
 int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
 

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

@@ -60 +60 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
+      // XXX accumulate all secfilt values in a single pass?
+
       int Ns;
       for (Ns = 0; Ns < Nphotcodes; Ns++) {
@@ -82 +84 @@
             continue;
           }
+
+          // XXX include the ubercal measurements?  yes, of course!  these are the well-determined values, supposedly!
+
           // XXX allow REF stars (no Image Entry) to be included in the calculation this
           // should be optionally set, and should allow for REF stars to be downweighted by
@@ -144 +149 @@
 }
 
+// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
+// * setMrel uses the internal Tiny structures only
+// * setMrelOutput skips stars for which there are too few good measurements 
+// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
 int setMrelOutput (Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr) {
 
@@ -167 +176 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
+
       /* skip stars already calibrated */
       if (catalog[i].found[j]) continue;  
@@ -172 +182 @@
       int Ns;
       for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
         int thisCode = photcodes[Ns][0].code;
         Nsec = GetPhotcodeNsec(thisCode);

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

@@ -189 +189 @@
   }
 
+  KEEP_UBERCAL = FALSE;
+  if ((N = get_argument (argc, argv, "-keep-ubercal"))) {
+    remove_argument (N, &argc, argv);
+    KEEP_UBERCAL = TRUE;
+  }
+
   MIN_ERROR = 0.001;
   if ((N = get_argument (argc, argv, "-minerror"))) {

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/help.c

@@ -37 +37 @@
   fprintf (stderr, "  -imfreeze\n");
   fprintf (stderr, "  -grid\n");
+  fprintf (stderr, "  -keep-ubercal : do not change zero point for ubercal-ed images\n");
   fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
   fprintf (stderr, "  -instmag min max\n");

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plot_scatter.c

@@ -11 +11 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Ntot += catalog[i].average[j].Nmeasure;
+      Ntot += catalog[i].averageT[j].Nmeasure;
     }
   }
@@ -32 +32 @@
             /* calculate the average value for a single star */
             if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-            m = catalog[i].average[j].measureOffset;
+            m = catalog[i].averageT[j].measureOffset;
 
-            for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+            for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
                 // skip measurements that do not match the current photcode
                 int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plotstuff.c

@@ -165 +165 @@
   graphdata[0].color = black;
   graphdata[0].lweight = 0;
-  graphdata[0].size = 0.5;
+  graphdata[0].size = 1.0;
 
   graphdata[0].xmin = dUNDEF;

branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c

@@ -120 +120 @@
   /* determine fit values */
   for (i = 0; i < NLOOP; i++) {
+    if (PLOTSTUFF) {
+      plot_scatter (catalog, Ncatalog, flatcorr); 
+    }
     setMrel  (catalog, Ncatalog, flatcorr);
-    setMcal  (catalog, FALSE);
+    if (PLOTSTUFF) {
+      plot_scatter (catalog, Ncatalog, flatcorr); 
+    }
+    setMcal  (catalog, FALSE, flatcorr);
     setMmos  (catalog, FALSE, flatcorr);
     setMgrid (catalog, flatcorr);
@@ -160 +166 @@
 
   /* set Mcal & Mmos for bad images */
-  setMcal  (catalog, TRUE);
+  setMcal  (catalog, TRUE, flatcorr);
   setMmos  (catalog, TRUE, flatcorr);
   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);