IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 4, 2011, 1:04:41 PM (15 years ago)
Author:
eugene
Message:

updates to pmPeak to better distinguish peak flux versions; updates to visualization; add bits for substantial suspect masking; consolidate assignment of source position and flux based on peak, moments, etc; improve footprint culling process; fix PSF_QF and PSF_QF_PERFECT calculations; fix source model chisq values

Location:
trunk/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:ignore

      •  

        old new  
        2828ChangeLog
        2929psmodules-*.tar.*
         30a.out.dSYM

  • trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r30621 r31153  
    5050
    5151#include "pmSourceIO.h"
     52#include "pmSourceOutputs.h"
    5253
    5354// panstarrs-style FITS table output (header + table in 1st extension)
    … …  
    6667    psArray *table;
    6768    psMetadata *row;
    68     psF32 *PAR, *dPAR;
    69     psEllipseAxes axes;
    70     psF32 xPos, yPos;
    71     psF32 xErr, yErr;
    72     psF32 errMag, chisq, apRadius;
    73     psS32 nPix, nDOF;
    7469
    7570    pmChip *chip = readout->parent->parent;
    76     pmFPA  *fpa  = chip->parent;
    77 
    78     bool status1 = false;
    79     bool status2 = false;
    80     float magOffset = NAN;
    81     float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
    82     float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    83     if (!isfinite(zeropt)) {
    84         zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    85     }
    86     if (status1 && status2 && (exptime > 0.0)) {
    87         magOffset = zeropt + 2.5*log10(exptime);
    88     }
    89     float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
    90 
    91     // if the sequence is defined, write these in seq order; otherwise
    92     // write them in S/N order:
     71
     72    // if the sequence is defined, write these in seq order; otherwise write them in S/N order.
     73    // Careful: if we are working with child sources, then we need to sort by the parent info,
     74    // not our info
    9375    if (sources->n > 0) {
    94         pmSource *source = (pmSource *) sources->data[0];
     76        pmSource *source = sources->data[0];
    9577        if (source->seq == -1) {
    96             // let's write these out in S/N order
    97             sources = psArraySort (sources, pmSourceSortBySN);
     78            sources = psArraySort (sources, pmSourceSortByFlux);
    9879        } else {
    9980            sources = psArraySort (sources, pmSourceSortBySeq);
    … …  
    10384    table = psArrayAllocEmpty (sources->n);
    10485
    105 # if (0)
    106     // we use this just to define the output vectors (which must be present for all objects)
    107     bool status = false;
    108     psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    109     psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    110     psAssert (radMax, "this must have been defined and tested earlier!");
    111     psAssert (radMax->n, "this must have been defined and tested earlier!");
    112     psAssert (radMin->n == radMax->n, "inconsistent annular bins");
    113 
    114     // write the radial profile apertures to header
    115     for (int i = 0; i < radMax->n; i++) {
    116       sprintf (keyword1, "RMIN_%02d", i);
    117       sprintf (keyword2, "RMAX_%02d", i);
    118       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
    119       psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
    120     }
    121 # endif
     86    short nImageOverlap;
     87    float magOffset;
     88    float zeroptErr;
     89    float fwhmMajor;
     90    float fwhmMinor;
     91    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    12292
    12393    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    12494    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    12595    for (int i = 0; i < sources->n; i++) {
    126         pmSource *source = (pmSource *) sources->data[i];
    127 
    128         // If source->seq is -1, source was generated in this analysis.  If source->seq is
    129         // not -1, source was read from elsewhere: in the latter case, preserve the source
    130         // ID.  source.seq is used instead of source.id since the latter is a const
    131         // generated on Alloc, and would thus be wrong for read in sources.
     96        // this is the source associated with this image
     97        pmSource *thisSource = sources->data[i];
     98
     99        // this is the "real" version of this source
     100        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     101
     102        // If source->seq is -1, source is unique and generated in this analysis.  If
     103        // source->seq is not -1, source was read from elsewhere or tied to other source (eg
     104        // from another image): in the latter case, preserve the source ID.  source.seq is used
     105        // instead of source.id since the latter is a const generated on Alloc, and would thus
     106        // be wrong for read in sources.
    132107        if (source->seq == -1) {
    133108            source->seq = i;
    134109        }
    135110
    136         // no difference between PSF and non-PSF model
    137         pmModel *model = source->modelPSF;
    138 
    139         if (model != NULL) {
    140             PAR = model->params->data.F32;
    141             dPAR = model->dparams->data.F32;
    142             xPos = PAR[PM_PAR_XPOS];
    143             yPos = PAR[PM_PAR_YPOS];
    144             if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
    145                 xErr = dPAR[PM_PAR_XPOS];
    146                 yErr = dPAR[PM_PAR_YPOS];
    147             } else {
    148                 // in linear-fit mode, there is no error on the centroid
    149                 xErr = source->peak->dx;
    150                 yErr = source->peak->dy;
    151             }
    152             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
    153                 axes = pmPSF_ModelToAxes (PAR, 20.0);
    154             } else {
    155                 axes.major = NAN;
    156                 axes.minor = NAN;
    157                 axes.theta = NAN;
    158             }
    159             chisq = model->chisq;
    160             nDOF = model->nDOF;
    161             nPix = model->nPix;
    162             apRadius = source->apRadius;
    163             errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    164         } else {
    165             xPos = source->peak->xf;
    166             yPos = source->peak->yf;
    167             xErr = source->peak->dx;
    168             yErr = source->peak->dy;
    169             axes.major = NAN;
    170             axes.minor = NAN;
    171             axes.theta = NAN;
    172             chisq = NAN;
    173             nDOF = 0;
    174             nPix = 0;
    175             apRadius = NAN;
    176             errMag = NAN;
    177         }
    178 
    179         float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
    180         float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
    181         psS16 nImageOverlap = 1;
    182 
    183         psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    184         float posAngle = 0.0;
    185         float pltScale = 0.0;
    186         pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     111        // set the 'best' values for various output fields:
     112        pmSourceOutputs outputs;
     113        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     114
     115        pmSourceOutputsMoments moments;
     116        pmSourceOutputsSetMoments (&moments, source);
    187117
    188118        row = psMetadataAlloc ();
    189119        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    190         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
    191         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
    192         psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
    193         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
    194         psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
    195         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
     120        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     121        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     122        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
     123        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     124        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     125        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    196126        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    197         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
     127        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
    198128        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
    199129        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
    200130        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    201131        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    202         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
    203         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
    204         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
     132        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     133        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     134        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    205135        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    206         psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
    207         psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
     136        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     137        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    208138        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    209139        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    210140
    211         psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
     141        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    212142        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    213143        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    214144
    215         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
    216         psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
    217         psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
     145        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     146        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     147        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    218148        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    219149        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    220         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
    221         psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
    222 
    223         // distinguish moments measure from window vs S/N > XX ??
    224         float Mxx = source->moments ? source->moments->Mxx : NAN;
    225         float Mxy = source->moments ? source->moments->Mxy : NAN;
    226         float Myy = source->moments ? source->moments->Myy : NAN;
    227 
    228         float Mrf  = source->moments ? source->moments->Mrf : NAN;
    229         float Mrh  = source->moments ? source->moments->Mrh : NAN;
    230         float Krf  = source->moments ? source->moments->KronFlux : NAN;
    231         float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    232 
    233         float Kinner = source->moments ? source->moments->KronFinner : NAN;
    234         float Kouter = source->moments ? source->moments->KronFouter : NAN;
    235 
    236         float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
    237         float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
    238         float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
    239         float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
    240 
    241         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
    242         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
    243         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
    244 
    245         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
    246         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
    247         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
    248         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
    249 
    250         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
    251         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
    252         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
    253         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
    254 
    255         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
    256         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
     150        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     151        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     152
     153        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     154        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     155        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     156
     157        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     158        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     159        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     160        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     161
     162        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     163        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     164        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     165        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     166        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
     167        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
    257168
    258169        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
    259170        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
    260 
    261 # if (0)
    262         // XXX if we have raw radial apertures, write them out here
    263         psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    264         psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
    265         psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    266         psVectorInit (radFlux,    NAN);
    267         psVectorInit (radFluxErr, NAN);
    268         psVectorInit (radFill,    NAN);
    269         if (!source->radial) goto empty_annuli;
    270         if (!source->radial->flux) goto empty_annuli;
    271         if (!source->radial->fill) goto empty_annuli;
    272         psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths");
    273         psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths");
    274 
    275         // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    276         for (int j = 0; j < source->radial->flux->n; j++) {
    277             radFlux->data.F32[j]    = source->radial->flux->data.F32[j];
    278             radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j];
    279             radFill->data.F32[j]    = source->radial->fill->data.F32[j];
    280         }
    281 
    282     empty_annuli:
    283         psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
    284         psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
    285         psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
    286         psFree (radFlux);
    287         psFree (radFluxErr);
    288         psFree (radFill);
    289 # endif
    290171
    291172        // XXX not sure how to get this : need to load Nimages with weight?
    … …  
    406287        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
    407288        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
    408         source->peak->flux = peakFlux;
     289        source->peak->rawFlux = peakFlux;
     290        source->peak->smoothFlux = peakFlux;
    409291        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
    410292        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
    411293        source->peak->dx   = dPAR[PM_PAR_XPOS];
    412294        source->peak->dy   = dPAR[PM_PAR_YPOS];
    413         if (isfinite (source->errMag) && (source->errMag > 0.0)) {
    414           source->peak->SN = 1.0 / source->errMag;
    415         } else {
    416           source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
    417         }
    418295
    419296        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    … …  
    429306
    430307        source->moments = pmMomentsAlloc ();
     308        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
     309        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
     310
    431311        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
    432312        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
    … …  
    500380
    501381    // let's write these out in S/N order
    502     sources = psArraySort (sources, pmSourceSortBySN);
     382    sources = psArraySort (sources, pmSourceSortByFlux);
    503383
    504384    table = psArrayAllocEmpty (sources->n);
    … …  
    522402    // we write out all sources, regardless of quality.  the source flags tell us the state
    523403    for (int i = 0; i < sources->n; i++) {
    524         pmSource *source = sources->data[i];
     404        // this is the source associated with this image
     405        pmSource *thisSource = sources->data[i];
     406
     407        // this is the "real" version of this source
     408        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    525409
    526410        // skip sources without measurements
    … …  
    556440        }
    557441        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
    558         psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
     442        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",            AxialTheta);
    559443
    560444        // Petrosian measurements
    … …  
    567451                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
    568452                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
    569                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
    570                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
    571                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
    572                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
     453                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude", mag);
     454                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", magErr);
     455                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
     456                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
    573457                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
    574458                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
    575459                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
    576460                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
     461                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill);
    577462            } else {
    578                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    579                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
    580                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
    581                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
     463                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude",       NAN);
     464                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", NAN);
     465                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius",          NAN);
     466                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error",    NAN);
    582467                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
    583468                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN);
    584469                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
    585470                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN);
     471                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", NAN);
    586472            }
    587473        }
    … …  
    672558
    673559    // let's write these out in S/N order
    674     sources = psArraySort (sources, pmSourceSortBySN);
     560    sources = psArraySort (sources, pmSourceSortByFlux);
    675561
    676562    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
    677563    int nParamMax = 0;
    678564    for (int i = 0; i < sources->n; i++) {
    679         pmSource *source = sources->data[i];
     565        // this is the source associated with this image
     566        pmSource *thisSource = sources->data[i];
     567
     568        // this is the "real" version of this source
     569        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     570
    680571        if (source->modelFits == NULL) continue;
    681572        for (int j = 0; j < source->modelFits->n; j++) {
    … …  
    691582    for (int i = 0; i < sources->n; i++) {
    692583
    693         pmSource *source = sources->data[i];
     584        pmSource *thisSource = sources->data[i];
     585
     586        // this is the "real" version of this source
     587        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    694588
    695589        // XXX if no model fits are saved, write out modelEXT?
    … …  
    747641
    748642                if (k < model->params->n) {
    749                     psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
     643                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
    750644                } else {
    751                     psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
     645                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN);
    752646                }
    753647            }
    754648
    755             // XXX other parameters which may be set.
    756             // XXX flag / value to define the model
    757             // XXX write out the model type, fit status flags
    758 
     649            // optionally, write out the covariance matrix values
     650            // XXX do I need to pad this to match the biggest covar matrix?
     651            if (model->covar) {
     652                for (int iy = 0; iy < model->covar->numCols; iy++) {
     653                    for (int ix = iy; ix < model->covar->numCols; ix++) {
     654                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
     655                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
     656
     657                    }
     658                }                   
     659            }
    759660            psArrayAdd (table, 100, row);
    760661            psFree (row);
    … …  
    790691bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
    791692{
    792 
     693    bool status = false;
    793694    psArray *table;
    794695    psMetadata *row;
    … …  
    796697    char keyword1[80], keyword2[80];
    797698
     699    // perform full non-linear fits / extended source analysis?
     700    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
     701        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
     702        return true;
     703    }
     704
    798705    // create a header to hold the output data
    799706    psMetadata *outhead = psMetadataAlloc ();
    … …  
    803710
    804711    // we use this just to define the output vectors (which must be present for all objects)
    805     bool status = false;
    806712    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    807713    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    … …  
    818724    }
    819725
     726    // the FWHM values are available if we measured a psf-matched convolved set
    820727    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    821     if (!fwhmValues) {
    822         psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
    823         return false;
    824     }
    825728
    826729    // let's write these out in S/N order
    827     sources = psArraySort (sources, pmSourceSortBySN);
     730    sources = psArraySort (sources, pmSourceSortByFlux);
    828731
    829732    table = psArrayAllocEmpty (sources->n);
    … …  
    832735    for (int i = 0; i < sources->n; i++) {
    833736
    834         pmSource *source = sources->data[i];
     737        // this is the source associated with this image
     738        pmSource *thisSource = sources->data[i];
     739
     740        // this is the "real" version of this source
     741        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    835742
    836743        // skip sources without radial aper measurements (or insufficient)
    837         if (source->radialAper == NULL) continue;
    838         psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
    839 
    840         for (int entry = 0; entry < fwhmValues->n; entry++) {
     744        if (source->radialAper == NULL) continue;
     745
     746        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
     747
     748        for (int entry = 0; entry < source->radialAper->n; entry++) {
    841749
    842750            // choose the convolved EXT model, if available, otherwise the simple one
    … …  
    844752            assert (radialAper);
    845753
    846             bool useMoments = true;
    847             useMoments = (useMoments && source->moments);          // can't if there are no moments
    848             useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
    849             useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
    850 
    851             if (useMoments) {
     754            if (pmSourcePositionUseMoments(source)) {
    852755                xPos = source->moments->Mx;
    853756                yPos = source->moments->My;
    … …  
    863766            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
    864767            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
    865             psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
     768            if (fwhmValues) {
     769                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
     770            } else {
     771                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
     772            }
    866773
    867774            // XXX if we have raw radial apertures, write them out here
    868             psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
    869             psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
    870             psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
     775            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     776            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
     777            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     778            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
    871779            psVectorInit (radFlux,    NAN);
    872780            psVectorInit (radFluxErr, NAN);
    … …  
    879787            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    880788            for (int j = 0; j < radialAper->flux->n; j++) {
    881                 radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
    882                 radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
    883                 radFill->data.F32[j]    = radialAper->fill->data.F32[j];
     789                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
     790                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
     791                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
     792                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
    884793            }
    885794
    886795        write_annuli:
    887             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     PS_DATA_VECTOR, "flux within annuli",    radFlux);
    888             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli",  radFluxErr);
    889             psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
     796            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
     797            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
     798            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
     799            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
    890800            psFree (radFlux);
    891801            psFree (radFluxErr);
     802            psFree (radFluxStdev);
    892803            psFree (radFill);
    893804
Note: See TracChangeset for help on using the changeset viewer.