Changeset 27532 for trunk/psphot/src/psphotEfficiency.c

Timestamp:

Mar 30, 2010, 1:33:47 PM (16 years ago)

Author:

eugene

Message:

some API re-work to support psphotForced and detection efficiency measurements; better visualization; CR masking (optional); better source size analysis; some motion towards stack photometry

File:

• trunk/psphot/src/psphotEfficiency.c (modified) (14 diffs)

trunk/psphot/src/psphotEfficiency.c

@@ -4 +4 @@
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
 
-//#define TESTING
-
-
-# if 0
+# define TESTING 0
 
 // Calculate the limiting magnitude for an image
@@ -19 +16 @@
                      float *covarFactor,// Covariance factor
                      const pmReadout *ro,     // Readout of interest
-                     const pmPSF *psf,        // Point-spread function
+                     pmPSF *psf,              // Point-spread function
                      float thresh,            // Threshold for source identification
                      float smoothSigma,       // Gaussian smoothing sigma
@@ -152 +149 @@
 }
 
-# endif
-
 bool psphotEfficiency (pmConfig *config, const pmFPAview *view)
 {
@@ -178 +173 @@
 bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
 {
-# if 0
     bool status = true;
 
@@ -206 +200 @@
 
     // Collect recipe information
-    float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x
-    float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y
+    float smoothNsigma = psMetadataLookupF32(&status, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
+    psAssert (status && isfinite(smoothNsigma), "Unable to find PEAKS_SMOOTH_NSIGMA in recipe (or invalid value)");
+
+    float thresh = psMetadataLookupF32(&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
+    psAssert (status && isfinite(thresh), "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe (or invalid value)");
+
+    psVector *magOffsets = psMetadataLookupVector(&status, recipe, "EFF.MAG"); // Magnitude offsets
+    psAssert (status, "Unable to find EFF.MAG F32 vector in recipe");
+    psAssert (magOffsets->type.type == PS_TYPE_F32, "Unable to find EFF.MAG F32 vector in recipe");
+
+    int numSources = psMetadataLookupS32(&status, recipe, "EFF.NUM"); // Number of sources for each bin
+    psAssert (status && (numSources > 0), "Unable to find EFF.NUM in recipe (or invalid value)");
+
+    float minGauss = psMetadataLookupF32(&status, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
+    psAssert (status && isfinite(minGauss), "PEAKS_MIN_GAUSS is not set in recipe (or invalid)");
+
+    // find the PSF size information (why is this not part of the psf structure?)
+    float fwhmMajor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MAJ"); // PSF size in x
+    float fwhmMinor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MIN"); // PSF size in y
     if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) {
-        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe");
-        return false;
-    }
-    float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit
-    if (!isfinite(smoothNsigma)) {
-        psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe");
-        return false;
-    }
-    float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2");
-    if (!isfinite(thresh)) {
-        psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe");
-        return false;
-    }
-    psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets
-    if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) {
-        psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe");
-        return NULL;
-    }
-    int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin
-    if (numSources == 0) {
-        psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe");
-        return NULL;
-    }
-    float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel
-    if (!isfinite(minGauss)) {
-        psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value");
-        minGauss = 0.5;
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in readout->analysis");
+        return false;
     }
 
@@ -246 +232 @@
     // remove all sources, adding noise for subtracted sources
     psphotRemoveAllSources(realSources, recipe);
-    // psphotAddNoise(readout, realSources, recipe);
-
-
-#if defined(TESTING) && 0
+
+#if TESTING
     {
         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
@@ -278 +262 @@
     }
 
-#ifdef TESTING
+#if TESTING
     psphotSaveImage(NULL, readout->image, "orig_image.fits");
     psphotSaveImage(NULL, readout->variance, "orig_variance.fits");
@@ -293 +277 @@
     }
 
-#ifdef TESTING
+#if TESTING
     psphotSaveImage(NULL, readout->image, "fake_image.fits");
     psphotSaveImage(NULL, readout->variance, "fake_variance.fits");
@@ -408 +392 @@
     psFree(significance);
 
-    if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) {
+    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
@@ -415 +399 @@
     }
 
-    // Disable aperture corrections; casting away const!
+    // Disable aperture corrections (save current value)
     pmTrend2D *apTrend = psf->ApTrend;  // Aperture trend
-    ((pmPSF*)psf)->ApTrend = NULL;
-
-    if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) {
+    psf->ApTrend = NULL;
+
+    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
         psFree(fakeSources);
         psFree(count);
-        ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const!
+        psf->ApTrend = apTrend; // Casting away const!
         return false;
     }
     psFree(fakeSourcesAll);
 
-    // Re-enable aperture corrections; casting away const!
-    ((pmPSF*)psf)->ApTrend = apTrend;
+    // Replace aperture corrections
+    psf->ApTrend = apTrend;
 
     psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin
@@ -443 +427 @@
         psVectorInit(magMask, 0);
 
-#ifdef TESTING
+#if TESTING
         psString name = NULL;
         psStringAppend(&name, "fake_%d.dat", i);
@@ -459 +443 @@
             }
 
-#ifdef TESTING
+#if TESTING
             fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf,
                     source->modelPSF->params->data.F32[PM_PAR_XPOS],
@@ -495 +479 @@
         }
 
-#ifdef TESTING
+#if TESTING
         fclose(file);
 #endif
@@ -520 +504 @@
     psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
 
-# endif
-
     return true;
 }