Changeset 5084 for trunk/psphot/src/pmPSFtry.c

Timestamp:

Sep 20, 2005, 6:48:16 PM (21 years ago)

Author:

eugene

Message:

alternative clipping for pmPSFtryMetric

File:

• trunk/psphot/src/pmPSFtry.c (modified) (4 diffs)

trunk/psphot/src/pmPSFtry.c

@@ -161 +161 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (try, RADIUS);
+    pmPSFtryMetric_Alt (try, RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
               modelName, try->psf->ApResid, try->psf->dApResid, try->psf->skyBias);
@@ -167 +167 @@
     return (try);
 }
-
 
 bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
@@ -195 +194 @@
   }
   fclose (f);
+
+  // XXX EAM : try 3hi/1lo sigma clipping on the rflux v dap fit
 
   // find min and max of (1/flux):
@@ -297 +298 @@
   return true;
 }
+
+bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) {
+
+  // the measured (aperture - fit) magnitudes (dA == try->metric)
+  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+  //     dA = dAo + dsky/flux
+  //   where flux is the flux of the star
+  // we fit this trend to find the infinite flux aperture correction (dAo),
+  //   the nominal sky bias (dsky), and the error on dAo
+  // the values of dA are contaminated by stars with close neighbors in the aperture
+  //   we use an outlier rejection to avoid this bias
+
+  // rflux = ten(0.4*fitMag);
+  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < try->sources->n; i++) {
+    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
+  }
+
+  // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+  psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+  // XXX EAM 
+  stats->min = 1.0;
+  stats->max = 3.0;
+  stats->clipIter = 3;
+
+  // linear fit to rfBin, daBin
+  psPolynomial1D *poly = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1);
+
+  // XXX EAM : this is the intended API (cycle 7? cycle 8?) 
+  poly = VectorClipFitPolynomial1D_EAM (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
+
+  fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+  // these stats should come back from the fit above...
+  // psVector *daFit   = Polynomial1DEvalVector_EAM (poly, rflux);
+  // psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) try->metric, "-", (void *) daFit);
+
+  // stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+  // stats->clipIter = 3;
+  // stats->clipSigma = 3;
+
+  // stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+  try->psf->ApResid = poly->coeff[0];
+  try->psf->dApResid = stats->sampleStdev;
+  try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+  FILE *f;
+  f = fopen ("apresid.dat", "w");
+  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
+
+  for (int i = 0; i < try->sources->n; i++) {
+    fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
+  }
+  fclose (f);
+
+  psFree (rflux);
+  psFree (poly);
+  psFree (stats);
+
+  // psFree (daFit);
+  // psFree (daResid);
+
+  return true;
+}