Changeset 25536

Timestamp:

Sep 23, 2009, 7:19:31 PM (17 years ago)

Author:

eugene

Message:

working on the better systematic error analysis

File:

• branches/eam_branches/20090715/psphot/src/psphotApResid.c (modified) (4 diffs)

branches/eam_branches/20090715/psphot/src/psphotApResid.c

@@ -192 +192 @@
         mask->data.PS_TYPE_VECTOR_MASK_DATA[Npsf] = 0;
 
+        // XXX don't I have errMag at this point??
         dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
 
@@ -410 +411 @@
     int Nx, Ny;
 
+    // XXX use the psf trend scale?
     if (readout->image->numCols > readout->image->numRows) {
         Nx = scale;
@@ -422 +424 @@
     }
 
-    // require at least 10 stars per spatial bin
-    if (Npsf < 10*Nx*Ny) {
-        return false;
-    }
-
     // the mask marks the values not used to calculate the ApTrend
     psVectorInit (mask, 0);
 
     // XXX stats structure for use by ApTrend : make parameters user setable
+    // XXX another stat?
     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     stats->min = 2.0;
@@ -514 +512 @@
 }
 
-# if (0)
-bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
-
-    int Nskip = 0;
-    int Nfail = 0;
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
-
-        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
-        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
-        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
-        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
-        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
-
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
-        }
-
-        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
-            continue;
-        }
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nskip = Nskip;
-    *nfail = Nfail;
+bool psphotApResidTrend_old (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
+
+    int Nx, Ny;
+
+    if (readout->image->numCols > readout->image->numRows) {
+        Nx = scale;
+        float AR = readout->image->numRows / (float) readout->image->numCols;
+        Ny = (int) (Nx * AR + 0.5);
+        Ny = PS_MAX (1, Ny);
+    } else {
+        Ny = scale;
+        float AR = readout->image->numRows / (float) readout->image->numCols;
+        Nx = (int) (Ny * AR + 0.5);
+        Nx = PS_MAX (1, Nx);
+    }
+
+    // require at least 10 stars per spatial bin
+    if (Npsf < 10*Nx*Ny) {
+        return false;
+    }
+
+    // the mask marks the values not used to calculate the ApTrend
+    psVectorInit (mask, 0);
+
+    // XXX stats structure for use by ApTrend : make parameters user setable
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    stats->min = 2.0;
+    stats->max = 3.0;
+
+    // measure Trend2D for the current spatial scale
+    psFree (psf->ApTrend);
+    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
+
+    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
+    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
+    for (int i = 0; i < dMag->n; i++) {
+        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
+    }
+
+    // XXX test for errors here
+    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
+
+    // construct the fitted values and the residuals
+    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
+    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
+
+    // measure systematic errorFloor & systematic / photon scale factor
+    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
+    int nGroup = PS_MAX (3*Nx*Ny, 10);
+    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
+
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
+
+    psFree (stats);
+    psFree (dMagSoft);
+    psFree (apResidFit);
+    psFree (apResidRes);
 
     return true;
 }
-# endif
+