Changeset 41359 for trunk/psphot/src/psphotSourceStats.c

Timestamp:

May 19, 2020, 2:53:41 PM (6 years ago)

Author:

eugene

Message:

better method to handle case where windowed second-momemts measurements yield no stars

File:

• trunk/psphot/src/psphotSourceStats.c (modified) (3 diffs)

trunk/psphot/src/psphotSourceStats.c

@@ -560 +560 @@
         }
 
-#if (0)
+#if (1)
         if (psTraceGetLevel("psphot.moments.save")) {
           char name[64];
@@ -604 +604 @@
 #endif
 
-        Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+        // a clump with no stars is not a valid clump
         Nout[i] = psfClump.nStars;
+        Sout[i] = (Nout[i] == 0) ? NAN : sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     }
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)
+    // EAM 2020.05.19 : this analysis fails if any of the measurements above have no valid stars
+    // A result with Nout[] = 0 should be skipped in the analysis
+    
     int Nstars = 0;
     float minKronRadius = NAN;
     float Sigma = NAN;
-    float minS = Sout[0];
-    float maxS = Sout[0];
+
+    float minS = NAN;
+    float maxS = NAN;
+    int maxN = -1;
+    int minN = -1;
     for (int i = 0; i < nsigma; i++) {
-        minS = PS_MIN(Sout[i], minS);
-        maxS = PS_MAX(Sout[i], maxS);
-    }
-    if (minS > 0.65) Sigma = sigma[nsigma-1];
-    if (maxS < 0.65) Sigma = sigma[0];
-
+      if (!isfinite(Sout[i])) continue;
+      if (!isfinite(minS)) {
+        minS = Sout[i];
+        maxS = Sout[i];
+        minN = i;
+        maxN = i;
+      }
+      if (Sout[i] < minS) {
+          minS = Sout[i];
+          minN = i;
+      }
+      if (Sout[i] > maxS) {
+          maxS = Sout[i];
+          maxN = i;
+      }
+    }
+    if (!isfinite(minS)) {
+      psLogMsg ("psphot", 3, "failed to find a valid window function\n");
+      return false;
+    }
+
+    // if all Sout values are above (or below) the threshold, saturate on the minimum (or maximum) valid sigma value
+    if (minS > 0.65) { Sigma = sigma[minN]; Nstars = Nout[minN]; minKronRadius = Rmin[minN]; }
+    if (maxS < 0.65) { Sigma = sigma[maxN]; Nstars = Nout[maxN]; minKronRadius = Rmin[maxN]; }
+
+    // this loop will be skipped if Sigma is set above
     for (int i = 0; i < nsigma - 1 && isnan(Sigma); i++) {
         if (!isfinite(Sout[i]) || !isfinite(Sout[i+1])) continue;
@@ -629 +656 @@
         minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]);
     }
-    psAssert (isfinite(Sigma), "did we miss a case?");
+    // EAM: this test will fail if only a single entry above is valid
+    if (!isfinite(Sigma)) {
+      psLogMsg ("psphot", 3, "failed to find a valid window function (case 2)\n");
+      return false;
+    }
     if (!isfinite(minKronRadius)) {
         minKronRadius = Sigma;