Changeset 28185

Timestamp:

Jun 1, 2010, 4:41:14 PM (16 years ago)

Author:

eugene

Message:

modify measurement rejection process

File:

• trunk/Ohana/src/relphot/src/StarOps.c (modified) (6 diffs)

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

@@ -391 +391 @@
 }
 
+# define NSIGMA_CLIP 3.0
+# define NSIGMA_REJECT 5.0
 void clean_measures (Catalog *catalog, int Ncatalog, int final) {
 
@@ -396 +398 @@
   int i, N, image_bad, TOOFEW;
   off_t *ilist;
-  double *tlist, *list, *dlist, Ns;
+  double *tlist, *list, *dlist;
   float Msys, Mcal, Mmos, Mgrid;
   StatType stats;
+  int Ncal, Nmos, Ngrid, Nfew;
 
   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
@@ -416 +419 @@
   TOOFEW = MAX (5, STAR_TOOFEW);
 
-  Ns = 3;
   Ndel = Nave = 0;
+  Ncal = Nmos = Ngrid = Nfew = 0;
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
@@ -433 +436 @@
         /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
         Mcal  = getMcal  (m, i);
-        if (isnan(Mcal)) continue;
+        if (isnan(Mcal)) { Ncal ++; continue; }
         Mmos  = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
+        if (isnan(Mmos)) { Nmos ++; continue; }
         Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
+        if (isnan(Mgrid)) { Ngrid ++; continue; }
 
         Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
@@ -444 +447 @@
         N++;
       }
-      if (N <= TOOFEW) continue;
+      if (N <= TOOFEW) { Nfew ++; continue; }
 
       /* 3-sigma clip based on stats of inner 50% */
+
+      // calculated mean of inner 50%
       initstats ("INNER_MEAN");
       liststats (list, dlist, N, &stats);
       stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
+
+      // ignore entries > 3sigma from inner mean
       for (k = m = 0; k < N; k++) {
-        if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
+        if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) {
           list[m] = list[k];
           m++;
         }
       }
+      // recalculate the mean & sigma of the accepted measurements
       initstats ("MEAN");
       liststats (list, dlist, m, &stats);
@@ -483 +491 @@
       if (N < TOOFEW) continue;
 
-      /* mark bad measures */
+      /* mark bad measures (> 3 sigma deviant) */
       for (k = 0; k < N; k++) {
-        if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
+        if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
           catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
           Ndel ++;