Changeset 6825

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

Location:

branches/rel10_ifa/psModules/src/objects

Files:

• pmModel.h (modified) (2 diffs)
• pmPSF.c (modified) (2 diffs)
• pmPSF.h (modified) (3 diffs)
• pmPSFtry.c (modified) (2 diffs)
• pmSourceFitModel.c (modified) (2 diffs)
• pmSourceFitSet.c (modified) (3 diffs)

branches/rel10_ifa/psModules/src/objects/pmModel.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-01 02:47:20 $
+ *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-10 20:21:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -42 +42 @@
     psVector *dparams;   ///< Parameter errors.
     float chisq;   ///< Fit chi-squared.
+    float chisqNorm;   ///< re-normalized fit chi-squared.
     int nDOF;    ///< number of degrees of freedom
     int nIter;    ///< number of iterations to reach min

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

@@ -6 +6 @@
  *  @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
@@ -97 +97 @@
     pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
 
+    if (PM_PSF_POISSON_WEIGHTS) {
+        psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+    } else {
+        psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    }
+
     // don't define a growth curve : user needs to choose radius bins
     psf->growth = NULL;

branches/rel10_ifa/psModules/src/objects/pmPSF.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.22.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-09 03:14:23 $
+ *  @version $Revision: 1.1.22.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-10 20:21:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -15 +15 @@
 # ifndef PM_PSF_H
 # define PM_PSF_H
+
+# define PM_PSF_POISSON_WEIGHTS 1
 
 /** pmPSF data structure
@@ -32 +34 @@
     pmModelType type;   ///< PSF Model in use
     psArray *params;   ///< Model parameters (psPolynomial2D)
+    psPolynomial1D *ChiTrend;  ///< Chisq vs flux fit (correction for systematic errors)
     psPolynomial4D *ApTrend;  ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)
     pmGrowthCurve *growth;  ///< apMag vs Radius
-    float ApResid;   ///< ???
-    float dApResid;   ///< ???
-    float skyBias;   ///< ???
-    float skySat;   ///< ???
-    float chisq;   ///< PSF goodness statistic
+    float ApResid;   ///< apMag - psfMag (for PSF stars)
+    float dApResid;   ///< scatter of ApResid
+    float skyBias;   ///< implied residual sky offset from ApResid fit
+    float skySat;   ///< roll-over of ApResid fit
+    float chisq;   ///< PSF goodness statistic (unused??)
     int nPSFstars;   ///< number of stars used to measure PSF
     int nApResid;   ///< number of stars used to measure ApResid

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

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

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-09 18:45:42 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-10 20:21:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -96 +96 @@
             y->data.F32[nPix] = source->pixels->data.F32[i][j];
             // psMinimizeLMChi2 takes wt = 1/dY^2
-            if (POISSON_WEIGHT) {
+            if (PM_PSF_POISSON_WEIGHTS) {
                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
             } else {

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

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-01 02:47:20 $
+ *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-10 20:21:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -182 +182 @@
     psS32 nPix = source->pixels->numRows * source->pixels->numCols;
 
+    // local sky variance
+    psF32 dSky = source->moments->dSky;
+
     // construct the coordinate and value entries
     psArray *x = psArrayAlloc(nPix);
@@ -206 +209 @@
             y->data.F32[nPix] = source->pixels->data.F32[i][j];
             // psMinimizeLMChi2 takes wt = 1/dY^2
-            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            if (PM_PSF_POISSON_WEIGHTS) {
+                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            } else {
+                yErr->data.F32[nPix] = 1.0 / dSky;
+            }
             nPix++;
         }