Changeset 35521 for branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c

Timestamp:

May 6, 2013, 10:32:51 AM (13 years ago)

Author:

eugene

Message:

remove superfulous ppStats.h includes; generate output psf model

File:

• branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c (modified) (3 diffs)

branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c

@@ -11 +11 @@
 
 #include "pswarp.h"
-#include <ppStats.h>
-#include "pswarpFileNames.h"            // Lists of file rules used at different stages
 
 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
@@ -19 +17 @@
 // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
 // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
-bool pswarpMakePSF (pmConfig *config) {
+
+bool pswarpMakePSF (pmConfig *config, pmFPAfile *output, psMetadata *stats) {
+
+    bool status = false;
 
     // load the recipe
@@ -28 +29 @@
     }
 
-    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
-        pswarpFileActivation(config, photFiles, true);
-        if (!pswarpIOChecksBefore(config)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
+    if (!psMetadataLookupBool(&status, recipe, "PSF")) {
+        psLogMsg("pswarp", PS_LOG_INFO, "Skipping PSF measurement");
+        return true;
+    }
 
-        // supply the readout and fpa of interest to psphot
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outFPA);
+    // XXX What was being read here??  skip this for now
+    if (false && !pswarpIOChecksBefore(config)) {
+        psError(psErrorCodeLast(), false, "Unable to read files.");
+        return false;
+    }
 
-        pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell
-        view->chip = view->cell = view->readout = 0;
+    // supply the readout and fpa of interest to psphot
+    pmFPAfile *photFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.INPUT");
+    pmFPACopy(photFile->fpa, output->fpa);
 
-        // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
-        psArray *sources = psphotLoadPSFSources (config, view);
-        if (!sources) {
-            psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
-            goto DONE;
-        }
+    pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell
 
-        pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, photFile->fpa, 1)) != NULL) {
+        psTrace ("pswarpMakePSF", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
 
-        // measure the PSF using these sources
-        if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
-            // This is likely a data quality issue
-            // XXX Split into multiple cases using error codes?
-            psErrorStackPrint(stderr, "Unable to determine PSF");
-            psWarning("Unable to determine PSF --- suspect bad data quality.");
-            if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) {
-                psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                                 "Unable to determine PSF", psErrorCodeLast());
-            }
-            psErrorClear();
-            psphotFilesActivate(config, false);
-        }
+        pmCell *cell;
+        while ((cell = pmFPAviewNextCell (view, photFile->fpa, 1)) != NULL) {
+            psTrace ("pswarpMakePSF", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
 
-        // Ensure seeing is carried over
-        pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);                 // Chip with seeing
-        psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing
-        item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING");
+            // process each of the readouts
+            pmReadout *readout;
+            while ((readout = pmFPAviewNextReadout(view, photFile->fpa, 1)) != NULL) {
+                if (!readout->data_exists) {
+                    continue;
+                }
 
-// XXX EAM : put this in a visualization function
-#if (TESTING)
-        {
-            #define PSF_SIZE 20         ///< Half-size of PSF
-            #define PSF_FLUX 10000      ///< Central flux for PSF
-            pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
-            pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
-            psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
-            psImageInit(image, 0);
-            pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
-            pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
-            psFree(model);
-            psFits *fits = psFitsOpen("psf.fits", "w");
-            psFitsWriteImage(fits, NULL, image, 0, NULL);
-            psFitsClose(fits);
-            psFree(image);
-        }
-#endif
+                // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
+                psArray *sources = psphotLoadPSFSources (config, view);
+                if (!sources) {
+                    // psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
+                    psLogMsg ("psphot", 3, "no psf sources for this readout (%d %d %d)", view->chip, view->cell, view->readout);
+                    continue;
+                }
 
-        psFree(view);
+                // XXX move above the loop
+                pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
+                
+                // measure the PSF using these sources
+                if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
+                    // A failure is likely a data quality issue.  raise a quality flag, but do not skip other readouts
+                    // XXX Split into multiple cases using error codes?
+                    psErrorStackPrint(stderr, "Unable to determine PSF for readout (%d %d %d)", view->chip, view->cell, view->readout);
+                    psWarning("Unable to determine PSF --- suspect bad data quality : readout (%d %d %d)", view->chip, view->cell, view->readout);
+                    if (stats) {
+                        int currentQuality = psMetadataLookupS32(&status, stats, "QUALITY");
+                        if (currentQuality == 0) {
+                            psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "Unable to determine PSF", psErrorCodeLast());
+                        }
+                    }
+                    psErrorClear();
+                }
+            }
+        }
+
+        // Ensure seeing is carried over from the PSPHOT.INPUT chip to the output chip
+        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Chip with seeing
+        psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing
+        item->data.F32 = psMetadataLookupF32(&status, chip->concepts, "CHIP.SEEING");
     }
+
+    // XXX test function found below
+    // pswarpMakePSF_test ();
+
+    psFree(view);
+    return true;
+}
+
+# if (0)
+bool pswarpMakePSF_test () {
+
+#define PSF_SIZE 20         ///< Half-size of PSF
+#define PSF_FLUX 10000      ///< Central flux for PSF
+    pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
+    pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
+    psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
+    psImageInit(image, 0);
+    pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
+    pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
+    psFree(model);
+    psFits *fits = psFitsOpen("psf.fits", "w");
+    psFitsWriteImage(fits, NULL, image, 0, NULL);
+    psFitsClose(fits);
+    psFree(image);
     return true;
 }
 # endif
-