Changeset 30621 for trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

Timestamp:

Feb 13, 2011, 11:59:42 AM (15 years ago)

Author:

eugene

Message:

do not cull peaks with only one entry; skip the section that culls duplicate peaks -- this was dropping many valid peaks and seems unneeded; save radial apertures in their own extension; fixed logic on culling of faint peaks

Location:

trunk/psModules/src/objects

Files:

• . (modified) (1 prop)
• pmSourceIO_CMF_PS1_SV1.c (modified) (7 diffs)

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -20 +20 @@
 
 #include "pmConfig.h"
+#include "pmErrorCodes.h"
 #include "pmDetrendDB.h"
 
@@ -71 +72 @@
     psF32 errMag, chisq, apRadius;
     psS32 nPix, nDOF;
-    char keyword1[80], keyword2[80];
 
     pmChip *chip = readout->parent->parent;
@@ -103 +103 @@
     table = psArrayAllocEmpty (sources->n);
 
+# if (0)
     // we use this just to define the output vectors (which must be present for all objects)
     bool status = false;
@@ -118 +119 @@
       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     }
-
+# endif
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -258 +259 @@
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
 
+# if (0)
+        // XXX if we have raw radial apertures, write them out here
         psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
         psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
@@ -284 +287 @@
         psFree (radFluxErr);
         psFree (radFill);
+# endif
 
         // XXX not sure how to get this : need to load Nimages with weight?
@@ -781 +785 @@
     return true;
 }
+
+// **** write out the radial flux values for the sources for a given matched-PSF image
+// **** how do we distinguish the matched-PSF images from the non-matched version
+bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+{
+
+    psArray *table;
+    psMetadata *row;
+    psF32 xPos, yPos;
+    char keyword1[80], keyword2[80];
+
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname);
+
+    // we use this just to define the output vectors (which must be present for all objects)
+    bool status = false;
+    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (radMax, "this must have been defined and tested earlier!");
+    psAssert (radMax->n, "this must have been defined and tested earlier!");
+    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
+
+    // write the radial profile apertures to header
+    for (int i = 0; i < radMax->n; i++) {
+      sprintf (keyword1, "RMIN_%02d", i);
+      sprintf (keyword2, "RMAX_%02d", i);
+      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
+      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
+    }
+
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (!fwhmValues) {
+        psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
+        return false;
+    }
+
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    table = psArrayAllocEmpty (sources->n);
+
+    // we write out all sources, regardless of quality.  the source flags tell us the state
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        // skip sources without radial aper measurements (or insufficient)
+        if (source->radialAper == NULL) continue;
+        psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
+
+        for (int entry = 0; entry < fwhmValues->n; entry++) {
+
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
+            assert (radialAper);
+
+            bool useMoments = true;
+            useMoments = (useMoments && source->moments);          // can't if there are no moments
+            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
+            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
+
+            if (useMoments) {
+                xPos = source->moments->Mx;
+                yPos = source->moments->My;
+            } else {
+                xPos = source->peak->xf;
+                yPos = source->peak->yf;
+            }
+
+            row = psMetadataAlloc ();
+
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
+
+            // XXX if we have raw radial apertures, write them out here
+            psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVectorInit (radFlux,    NAN);
+            psVectorInit (radFluxErr, NAN);
+            psVectorInit (radFill,    NAN);
+            if (!radialAper->flux) goto write_annuli;
+            if (!radialAper->fill) goto write_annuli;
+            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
+
+            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
+            for (int j = 0; j < radialAper->flux->n; j++) {
+                radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
+                radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
+                radFill->data.F32[j]    = radialAper->fill->data.F32[j];
+            }
+
+        write_annuli:
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
+            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
+            psFree (radFlux);
+            psFree (radFluxErr);
+            psFree (radFill);
+
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
+    }
+
+    if (table->n == 0) {
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+    return true;
+}