Changeset 6825 for branches/rel10_ifa/psModules/src/objects/pmPSFtry.c

Timestamp:

Apr 10, 2006, 10:21:36 AM (20 years ago)

Author:

magnier

Message:

adding chisq-flux trend corrections; adding option to use non-poisson errors for fits

File:

• branches/rel10_ifa/psModules/src/objects/pmPSFtry.c (modified) (2 diffs)

branches/rel10_ifa/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.4.4.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-09 03:14:23 $
+ *  @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-10 20:21:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -184 +184 @@
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
 
-    // make this optional
-    // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
+    // measure the chi-square trend as a function of flux (PAR[1])
+    // this should be linear for Poisson errors and quadratic for constant sky errors
+    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_MASK);
+
+    // write sources with models first
+    for (int i = 0; i < psfTry->sources->n; i++) {
+        pmModel *model = psfTry->modelPSF->data[i];
+        if (model == NULL) {
+            flux->data.F64[i] = 0.0;
+            chisq->data.F64[i] = 0.0;
+            mask->data.U8[i] = 1;
+        } else {
+            flux->data.F64[i] = model->params->data.F32[1];
+            chisq->data.F64[i] = model[0].chisq/model[0].nDOF;
+            mask->data.U8[i] = 0;
+        }
+    }
+    // use 3hi/3lo sigma clipping on the r2rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->clipSigma = 3.0;
+    stats->clipIter = 3;
+    // linear clipped fit of ApResid to r2rflux
+    psVectorClipFitPolynomial1D (psfTry->psf->ChiTrend, stats, mask, 1, chisq, NULL, flux);
+    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
+        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), psfTry->psf->ChiTrend->coeffErr[i]);
+    }
+    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    psFree (stats);
+    psFree (flux);
+    psFree (chisq);
 
     // XXX this function wants aperture radius for pmSourcePhotometry