Changeset 33394

Timestamp:

Mar 2, 2012, 12:20:43 PM (14 years ago)

Author:

eugene

Message:

fix ubercal distance; make mosaic work for ubercal; add M_20, M_80

Location:

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

Files:

• include/relphot.h (modified) (2 diffs)
• src/ImageOps.c (modified) (8 diffs)
• src/MosaicOps.c (modified) (8 diffs)
• src/StarOps.c (modified) (12 diffs)
• src/liststats.c (modified) (1 diff)
• src/relphot.c (modified) (1 diff)

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

@@ -25 +25 @@
   short Xm;
   float secz;
-  char flags;
+  float ubercalDist;
+  unsigned int flags;
   Coords coords;
 } Mosaic; 
@@ -37 +38 @@
   double min;
   double max;
+  double Upper80;
+  double Lower20;
   double total;
   int    Nmeas;

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

@@ -370 +370 @@
 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, Mflat, McalBright, McalBright2;
-  double *list, *dlist;
+  off_t i, j, m, c, n, N, Nmax;
+  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys, Nbright;
+  float Msys, Mrel, Mmos, Mgrid, Mflat;
+  double *list, *dlist, *Mlist, *dMlist;
   StatType stats;
 
@@ -383 +383 @@
 
   if (PoorImages) {
-    IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
+    // XXX use bad stars and measurements for PoorImages? or not?
+    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
+    IMAGE_BAD = 0;
   }
 
@@ -392 +394 @@
   ALLOCATE (list, double, Nmax);
   ALLOCATE (dlist, double, Nmax);
+  ALLOCATE (Mlist, double, Nmax);
+  ALLOCATE (dMlist, double, Nmax);
 
   Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
 
-  // counters to measure the bright-end scatter
-  McalBright = McalBright2 = 0.0;
-  Nbright = 0;
-
   for (i = 0; i < Nimage; i++) {
     
+    // XXXX TEST
+    Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
+
     /* on PoorImages run, skip good images */
     if (PoorImages) {
@@ -415 +418 @@
 
     int minUbercalDist = 1000;
+
+    // number of stars to measure the bright-end scatter
+    Nbright = 0;
 
     N = 0;
@@ -469 +475 @@
       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-        McalBright += list[N];
-        McalBright2 += SQ(list[N]);
+        Mlist[Nbright] = list[N];
+        dMlist[Nbright] = dlist[N];
         Nbright ++;
       }
@@ -476 +482 @@
     }
     /* N_onImage[i] is all measurements, N is good measurements */
+
+    // XXX TEST
+    fprintf (stderr, "meas skipped: (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ngrid, Nrel, Nsys);
 
     /* too few good measurements or too many bad measurements */
@@ -487 +496 @@
       }      
     }
-    
+
     liststats (list, dlist, N, &stats);
-
     image[i].Mcal  = stats.mean;
     image[i].dMcal = stats.error;
-    image[i].dMagSys = stats.sigma;
     image[i].nFitPhotom = N;
     image[i].Xm    = 100.0*log10(stats.chisq);
@@ -498 +505 @@
     plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
 
+    // bright end scatter
+    liststats (Mlist, dMlist, Nbright, &stats);
+    image[i].dMagSys = stats.sigma;
+
     if (image[i].Mcal < -CLOUD_TOLERANCE) {
       image[i].Mcal = 0.0;
     }
 
-    image[i].dummy3 = minUbercalDist;
-    fprintf (stderr, "%s : %d\n", image[i].name, image[i].dummy3);
+    // minUbercalDist calculated here is the min value for any star owned by this image
+    // since this particular image is tied to that star, bump its distance by 1
+    image[i].dummy3 = minUbercalDist + 1;
   }
   free (list);
   free (dlist);
+  free (Mlist);
+  free (dMlist);
 
   fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys);

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

@@ -1 +1 @@
 # include "relphot.h"
+void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 
 // array of mosaic definition structures
@@ -391 +392 @@
 int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
-  off_t i, j, m, c, N, Nmax;
-  int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt;
+  off_t i, j, m, c, n, N, Nmax;
+  int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright;
   float Msys, Mrel, Mcal, Mgrid, Mflat;
   double *list, *dlist, *Mlist, *dMlist;
@@ -403 +404 @@
   image = getimages (&N, NULL);
 
-  Nsecfilt = GetPhotcodeNsecfilt ();
+  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   if (PoorImages) {
+    // XXX use bad stars and measurements for PoorImages? or not?
+    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     IMAGE_BAD = 0;
   }
@@ -428 +433 @@
     }      
 
-    /* we optionally do not recalibrate images with UBERCAL zero points */
-    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+    // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
+    // we optionally do not recalibrate images with UBERCAL zero points 
+    if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
+      mosaic[i].ubercalDist = 0;
+      // propagate ubercalDist to the images
+      for (j = 0; j < MosaicN_Image[i]; j++) {
+        off_t im = MosaicToImage[i][j];
+        image[im].dummy3 = mosaic[i].ubercalDist;
+        fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].dummy3);
+      }
+      if (KEEP_UBERCAL) continue;
+    }
+
+    int minUbercalDist = 1000;
+
+    // number of stars to measure the bright-end scatter
+    Nbright = 0;
 
     N = 0;
@@ -470 +490 @@
         continue;
       }
+
+      PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
+      if (!code) goto skip;
+      if (code->equiv < 1) goto skip;
+      int Nsec = GetPhotcodeNsec (code->equiv);
+      if (Nsec == -1) goto skip;
+      minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].M_80, minUbercalDist); // M_80 == ubercalDist
+
+    skip:
       list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      Mlist[N] = Msys;
-      dMlist[N] = list[N];
+      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
+        Mlist[Nbright] = list[N];
+        dMlist[Nbright] = dlist[N];
+        Nbright ++;
+      }
       N++;
     }
@@ -479 +511 @@
 
     /* too few good measurements or too many bad measurements (skip in PoorImages run) */
-
     if (!PoorImages) {
       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[i]);
@@ -493 +524 @@
     if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
 
-    float CLOUD_TOLERANCE = 0.01;
     mosaic[i].Mcal  = stats.mean;
     mosaic[i].dMcal = stats.error;
-    mosaic[i].dMsys = stats.sigma;
     mosaic[i].nFitPhotom = N;
     mosaic[i].Xm    = 100.0*log10(stats.chisq);
 
+    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
+
+    // bright end scatter
+    liststats (Mlist, dMlist, Nbright, &stats);
+    mosaic[i].dMsys = stats.sigma;
+
     if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
       mosaic[i].Mcal = 0.0;
+    }
+
+    // minUbercalDist calculated here is the min value for any star owned by this image
+    // since this particular image is tied to that star, bump its distance by 1
+    mosaic[i].ubercalDist = minUbercalDist + 1;
+
+    // propagate ubercalDist to the images
+    for (j = 0; j < MosaicN_Image[i]; j++) {
+      off_t im = MosaicToImage[i][j];
+      image[im].dummy3 = mosaic[i].ubercalDist;
+      fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].dummy3);
     }
   }
@@ -875 +921 @@
   Nmark = Nscatter = Noffset = 0;
   for (i = 0; i < Nmosaic; i++) {
+    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+
     mark = FALSE;
     if (mosaic[i].dMcal > MaxScatter) {

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

@@ -61 +61 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-
       // XXX accumulate all secfilt values in a single pass?
 
@@ -100 +99 @@
 
           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) {
-            Nsys++;
-            continue;
-          }
+          if (isnan(Msys)) { Nsys++; continue; }
           list[N] = Msys - Mcal - Mmos - Mgrid;
 
@@ -194 +190 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
+      int minUbercalDist = 1000;
+
       int Ns;
       for (Ns = 0; Ns < Nphotcodes; Ns++) {
@@ -229 +227 @@
 
           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
+          if (isnan(Msys)) { continue; }
           list[N] = Msys - Mcal - Mmos - Mgrid;
+
+          int myUbercalDist = getUbercalDist(m,i);
+          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
 
           // dlist gives the error, which is used as the weight in WT_MEAN.
@@ -258 +260 @@
         if (N < 1) continue;
 
-        // force WT_MEAN or MEAN here?
-
+        // XXX force WT_MEAN or MEAN here?
         liststats (list, dlist, N, &stats);
         catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
@@ -286 +287 @@
         catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
         catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
+
+        catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = stats.Upper80; // M_80 == ubercalDist
+        catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = stats.Lower20; // M_80 == ubercalDist
+        catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = minUbercalDist; // M_80 == ubercalDist
       }
     }
@@ -465 +470 @@
 void clean_stars (Catalog *catalog, int Ncatalog) {
 
-  int i, j, Ndel, Nave, Ntot, mark, Ns;
+  int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
   float dM, Xm;
   double Chisq, MaxScatter, MaxChisq;
@@ -482 +487 @@
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // XX int oldPLOTSTUFF = PLOTSTUFF;
+  // XX PLOTSTUFF = TRUE;
+  // XX plot_chisq (catalog, Ncatalog);
+  // XX PLOTSTUFF = oldPLOTSTUFF;
 
   // eliminate bad stars using the stats for a single secfilt at a time
@@ -510 +520 @@
     fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
 
-    Ndel = Nave = 0;
+    Ndel = Nave = Nscat = Nnan = Nchi = 0;
     for (i = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
@@ -520 +530 @@
           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 ++; }
         } else {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
@@ -526 +539 @@
       }
     }
-    fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
+    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
     initstats (STATMODE);
   }
@@ -564 +577 @@
   TOOFEW = MAX (5, STAR_TOOFEW);
 
-  FILE *f = fopen ("poor.dat", "w");
+  FILE *f = final ? fopen ("poor.dat", "a") : fopen ("poor.dat", "w");
 
   Ndel = Nave = 0;

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

@@ -41 +41 @@
   stats[0].min    = value[0];
   stats[0].max    = value[N-1];
+  int N80 = MIN (N-1, 0.8*N);
+  int N20 = MAX (0, 0.2*N);
+  stats[0].Upper80 = value[N80];
+  stats[0].Lower20 = value[N20];
 
   switch (statmode) {

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

@@ -237 +237 @@
   freeGridBins (Ncatalog);
 
+  // XXX
+  FILE *f = fopen ("poor.dat", "w");
+  fclose (f);
+
   /* load catalog data from region files, update Mrel include all data */
   reload_catalogs (skylist, flatcorr, 0, NULL);