Changeset 25246

Timestamp:

Sep 2, 2009, 9:58:54 AM (17 years ago)

Author:

Paul Price

Message:

Getting ready to run it.

Location:

branches/pap/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotFake.c (modified) (10 diffs)
• psphotOutput.c (modified) (1 diff)
• psphotReadout.c (modified) (2 diffs)

branches/pap/psphot/src/Makefile.am

@@ -128 +128 @@
         psphotCheckStarDistribution.c  \
         psphotThreadTools.c            \
-        psphotAddNoise.c
+        psphotAddNoise.c               \
+        psphotFake.c
 
 # dropped? psphotGrowthCurve.c

branches/pap/psphot/src/psphot.h

@@ -76 +76 @@
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
+bool            psphotFake(pmReadout *readout, const pmPSF *psf, const psMetadata *recipe, const psArray *realSources);
 
 // thread-related:

branches/pap/psphot/src/psphotFake.c

@@ -20 +20 @@
                       float smoothNsigma,       // Smoothing limit
                       psImageMaskType maskVal   // Value to mask
-                      );
+                      )
 {
     PM_ASSERT_READOUT_NON_NULL(ro, false);
     PM_ASSERT_READOUT_VARIANCE(ro, false);
-    PS_ASSERT_METADATA_NON_NULL(recipe, false);
-
-    float smoothSigma  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrtf(2.0*log(2.0)));
+
+    float smoothSigma  = 0.5*(xFWHM + yFWHM) / (2.0*sqrtf(2.0*log(2.0)));
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
-    psKernel *newCovar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
+    psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix
     psFree(kernel);
     float factor = psImageCovarianceFactor(newCovar); // Variance factor
@@ -50 +49 @@
     }
     if (minFlux) {
-        int fudge = psMetadataLookupF32(NULL, recipe, "FAKE.MINFLUX"); // Fudge factor for minimum flux
-        *minFlux = sqrtf(meanVar) * fudge;
+        *minFlux = sqrtf(meanVar);
     }
 
@@ -60 +58 @@
 static bool fakeGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
                          const pmReadout *ro,            // Readout of interest
+                         const pmPSF *psf,               // Point-spread function
                          const psVector *magOffsets,     // Magnitude offsets for fake sources
                          int numSources,                 // Number of fake sources for each bin
@@ -68 +67 @@
     PM_ASSERT_READOUT_NON_NULL(ro, false);
     PM_ASSERT_READOUT_IMAGE(ro, false);
-    PS_ASSERT_METADATA_NON_NULL(recipe, false);
     PS_ASSERT_VECTOR_NON_NULL(magOffsets, false);
     PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false);
 
-    int numBins = magOffset->n;                                     // Number of bins
+    int numBins = magOffsets->n;                                    // Number of bins
     int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
 
@@ -84 +82 @@
 
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
-    for (int i = 0; i <= numBins; i++) {
-        psArray *data = psArrayAlloc(3); // Sources for bin
-
-        psVector *x = psVectorAlloc(numSources, PS_TYPE_F32);
-        psVector *y = psVectorAlloc(numSources, PS_TYPE_F32);
-
-        float mag = refMag + magOffset->data.F32[i]; // Instrumental magnitude of sources
-
-        for (int j = 0; j <= numSources; j++) {
-            x->data.F32[j] = psRandomUniform(rng) * numCols;
-            y->data.F32[j] = psRandomUniform(rng) * numRows;
-            magAll->data.F32[i * numSources + j] = mag;
+    for (int i = 0, index = 0; i <= numBins; i++) {
+        float mag = refMag + magOffsets->data.F32[i]; // Instrumental magnitude of sources
+
+        for (int j = 0; j <= numSources; j++, index++) {
+            xAll->data.F32[index] = psRandomUniform(rng) * numCols;
+            yAll->data.F32[index] = psRandomUniform(rng) * numRows;
+            magAll->data.F32[index] = mag;
         }
-        memcpy(xSrc->data.F32[i], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        memcpy(ySrc->data.F32[i], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        memcpy(&xAll->data.F32[i * numSources], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
-        memcpy(&yAll->data.F32[i * numSources], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy((*xSrc)->data.F32[i], &xAll->data.F32[i * numSources],
+               numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy((*ySrc)->data.F32[i], &yAll->data.F32[i * numSources],
+               numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
     }
     psFree(rng);
 
-    pmReadout *fakeRO = pmReadoutAlloc(); // Fake readout
-    if (!pmReadoutFakeFromVectors(fakeRO, xAll, yAll, mag, NULL, NULL, psf, minFlux, NAN, false, false)) {
+    pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
+    if (!pmReadoutFakeFromVectors(fakeRO, xAll, yAll, magAll, NULL, NULL, psf, minFlux, NAN, false, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image");
         psFree(fakeRO);
@@ -125 +118 @@
 
 // *** in this section, perform the photometry for fake sources ***
-bool psphotFake(pmConfig *config, pmReadout *readout, const pmPSF *psf,
-                const psMetadata *recipe, const psArray *realSources)
+bool psphotFake(pmReadout *readout, const pmPSF *psf, const psMetadata *recipe, const psArray *realSources)
 {
-    PS_ASSERT_PTR_NON_NULL(config, false);
     PM_ASSERT_READOUT_NON_NULL(readout, false);
     PM_ASSERT_READOUT_IMAGE(readout, false);
@@ -138 +129 @@
 
     // Collect recipe information
+    bool mdok;                                                 // Status of MD lookup
     float xFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_X"); // PSF size in x
     float yFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_Y"); // PSF size in y
@@ -155 +147 @@
     }
     psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "FAKE.MAG"); // Magnitude offsets
-    if (!magOffset || magOffset->type.type != PS_TYPE_F32) {
+    if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) {
         psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE.MAG F32 vector in recipe");
         return NULL;
@@ -164 +156 @@
         return NULL;
     }
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Value to mask
-
+    float minFluxFudge = psMetadataLookupF32(&mdok, recipe, "FAKE.MINFLUX"); // Fudge factor for minimum flux
+    if (!mdok) {
+        minFluxFudge = 1.0;
+    }
+
+    psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
 
     // remove all sources, adding noise for subtracted sources
@@ -177 +173 @@
         return false;
     }
+    minFlux *= minFluxFudge;
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
-    if (!fakeGenerate(&xFake, &yFake, readout, magOffsets, numSources, magLim, minFlux)) {
+    if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, minFlux)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
         psFree(xFake);

branches/pap/psphot/src/psphotOutput.c

@@ -228 +228 @@
     psMetadataItemSupplement (header, recipe, "MSKY_NY");
 
+    // fake source recovery statistics
+    psMetadataItemSupplement(header, recipe, "FAKE.EFF");
+    psMetadataItemSupplement(header, recipe, "FAKE.MAG");
+    psMetadataItemSupplement(header, recipe, "FAKE.REF");
+
     psMetadataAddF32 (header, PS_LIST_TAIL, "DT_PHOT", PS_META_REPLACE, "elapsed psphot time", psTimerMark ("psphotReadout"));
 

branches/pap/psphot/src/psphotReadout.c

@@ -61 +61 @@
     // display the backsub and backgnd images
     psphotVisualShowBackground (config, view, readout);
-    
+
     // run a single-model test if desired (exits from here if test is run)
     psphotModelTest (config, view, recipe);
@@ -224 +224 @@
     psphotMagnitudes(config, readout, view, sources, psf);
 
+    if (!psphotFake(readout, psf, recipe, sources)) {
+        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
+        psErrorClear();
+    }
+
     // replace failed sources?
     // psphotReplaceUnfitSources (sources);