IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects/pmSourceIO_CMF.c.in

    r35768 r37066  
    3737#include "pmMoments.h"
    3838#include "pmModelFuncs.h"
     39#include "pmModelClass.h"
    3940#include "pmModel.h"
    4041#include "pmModelUtils.h"
    41 #include "pmModelClass.h"
    4242#include "pmSourceMasks.h"
    4343#include "pmSourceExtendedPars.h"
    4444#include "pmSourceDiffStats.h"
    4545#include "pmSourceSatstar.h"
     46#include "pmSourceLensing.h"
    4647#include "pmSource.h"
    4748#include "pmSourceFitModel.h"
    … …  
    141142        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    142143        @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    143         @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     144        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              source->apRadius);
    144145        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
    145146        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
     147        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX",          PS_DATA_S32, "aperture unmasked pixels",                   source->apNpixels);
    146148
    147149        @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    … …  
    162164        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    163165
     166        // PSF shape parameters:
    164167        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
    165168        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
    166169        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
     170        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE",         PS_DATA_F32, "k term if defined",                          outputs.psfCore);
     171
     172        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
     173        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ",        PS_DATA_F32, "PSF FWHM (major axis)",                   outputs.psfMajorFWHM);
     174        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN",        PS_DATA_F32, "PSF FWHM (minor axis)",                   outputs.psfMinorFWHM);
     175
     176        // psf data quality
    167177        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    168178        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    … …  
    178188        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
    179189        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     190
     191        // Lensing parameters:
     192        if (source->lensingOBJ && source->lensingOBJ->smear) {
     193          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
     194          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
     195          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
     196          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
     197          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
     198        }
     199
     200        if (source->lensingOBJ && source->lensingOBJ->shear) {
     201          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
     202          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
     203          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
     204          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1);
     205          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2);
     206        }
     207
     208        if (source->lensingPSF && source->lensingPSF->smear) {
     209          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
     210          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
     211          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
     212          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1);
     213          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2);
     214        }
     215
     216        if (source->lensingPSF && source->lensingPSF->shear) {
     217          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
     218          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
     219          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
     220          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1);
     221          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2);
     222        }
    180223
    181224        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
    … …  
    268311
    269312    // define PSF model type
    270     // XXX need to carry the extra model parameters
    271313    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
    272314
    273     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     315    // if header does not define the model, default to a gaussian
     316    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
    274317    if (PSF_NAME != NULL) {
    275318        modelType = pmModelClassGetType (PSF_NAME);
    276319    }
    277320    assert (modelType > -1);
     321
     322    // do we expect to find lensing parameters?
     323    bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ");
     324    bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF");
    278325
    279326    // We get the size of the table, and allocate the array of sources first because the table
    … …  
    310357        @ALL@     axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");
    311358        @ALL@     axes.theta        = axes.theta * PS_RAD_DEG;
     359       
     360        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
     361        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
     362        @>PS1_V4,>PS1_SV2,>PS1_DV3@ }
    312363
    313364        @ALL@     PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
    … …  
    321372        @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
    322373        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
     374        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX");
     375        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
    323376
    324377        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
    … …  
    331384        @ALL@     dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
    332385
    333         pmPSF_AxesToModel (PAR, axes, modelType);
     386        pmPSF_AxesToModel (PAR, axes, model->class->useReff);
    334387
    335388        @ALL@     float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
    … …  
    349402        @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    350403        @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    351         @ALL@     source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
     404        @ALL@     source->apRadius  = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS");
     405        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
    352406
    353407        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    376430        @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
    377431        @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0;
     432
     433        // Lensing parameters (on read if PS1_V5+)
     434        if (haveLensOBJ) {
     435          source->lensingOBJ = pmSourceLensingAlloc ();
     436          source->lensingOBJ->smear = pmLensingParsAlloc();
     437          source->lensingOBJ->shear = pmLensingParsAlloc();
     438
     439          @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
     440          @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
     441          @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
     442          @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
     443          @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
     444          @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
     445          @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
     446          @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
     447          @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
     448          @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
     449        }
     450
     451        if (haveLensPSF) {
     452          source->lensingPSF = pmSourceLensingAlloc ();
     453          source->lensingPSF->smear = pmLensingParsAlloc();
     454          source->lensingPSF->shear = pmLensingParsAlloc();
     455
     456          @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
     457          @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
     458          @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
     459          @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
     460          @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
     461          @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
     462          @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
     463          @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
     464          @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
     465          @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
     466        }
    378467
    379468        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
    … …  
    448537    }
    449538
     539#ifdef SORT_EXTENSIONS_BY_FLUX
    450540    // let's write these out in S/N order
    451541    sources = psArraySort (sources, pmSourceSortByFlux);
     542#endif
    452543
    453544    table = psArrayAllocEmpty (sources->n);
    … …  
    708799            return false;
    709800        }
    710         // Find the source with this sequence number.
    711         // XXX: I am assuming that sources is sorted in order of seq
     801        // Find the source with this sequence number using the sourceIndex.
    712802        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
    713         pmSource *source = NULL;
    714 #ifndef ASSUME_SORTED
    715         long j = seq < sources->n ? seq : sources->n - 1;
    716         for (; j >= 0; j--) {
    717             source = sources->data[j];
    718             if (source->seq == seq) {
    719                 break;
    720             }
    721         }
    722 #else
    723803        long j = sourceIndex[seq];
    724804        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
    725         source = sources->data[j];
    726 #endif
     805        pmSource *source = sources->data[j];
    727806        if (!source) {
    728807            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
    … …  
    789868    char name[64];
    790869
     870    pmModelType modelTypeTrail = pmModelClassGetType("PS_MODEL_TRAIL");
     871
    791872    // create a header to hold the output data
    792873    psMetadata *outhead = psMetadataAlloc ();
    793874
     875    pmModelClassWriteHeader(outhead);
     876
    794877    // write the links to the image header
    795878    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
    796879
     880#ifdef SORT_EXTENSIONS_BY_FLUX
    797881    // let's write these out in S/N order
    798882    sources = psArraySort (sources, pmSourceSortByFlux);
    799 
    800     @>PS1_DV2@ float magOffset;
    801     @>PS1_DV2@ float zeroptErr;
    802     @>PS1_DV2@ float fwhmMajor;
    803     @>PS1_DV2@ float fwhmMinor;
    804     @>PS1_DV2@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     883#endif
     884
     885    float magOffset;
     886    float zeroptErr;
     887    float fwhmMajor;
     888    float fwhmMinor;
     889    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    805890
    806891    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
    … …  
    823908    @>PS1_DV2@ pmChip *chip = readout->parent->parent;
    824909
     910    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
     911    badModel |= PM_MODEL_STATUS_BADARGS;
     912    badModel |= PM_MODEL_STATUS_OFFIMAGE;
     913    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
     914    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
     915    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
     916    badModel |= PM_MODEL_PCM_FAIL_GUESS;
     917
    825918    table = psArrayAllocEmpty (sources->n);
    826919
    … …  
    847940
    848941            // skip models which were not actually fitted
    849             if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
     942            // XXX
     943            if (model->flags & badModel) continue;
    850944
    851945            PAR = model->params->data.F32;
    … …  
    853947            xPos = PAR[PM_PAR_XPOS];
    854948            yPos = PAR[PM_PAR_YPOS];
    855             xErr = dPAR[PM_PAR_XPOS];
    856             yErr = dPAR[PM_PAR_YPOS];
     949
     950            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
     951            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
     952            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
     953            // TRAIL : X,Y are fitted
     954            //
     955           
     956            // XXX this should be based on what happened, not on the model type
     957            if (model->type == modelTypeTrail) {
     958                xErr = dPAR[PM_PAR_XPOS];
     959                yErr = dPAR[PM_PAR_YPOS];
     960            } else {
     961                // this is definitely an underestimate since it does not
     962                // account for the extent of the source
     963                xErr = fwhmMajor * model->magErr / 2.35;
     964                yErr = fwhmMinor * model->magErr / 2.35;
     965            }
    857966
    858967            @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    … …  
    870979            row = psMetadataAlloc ();
    871980
    872             // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    873981            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
    874982            // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
    … …  
    888996            @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
    889997            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
    890             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Magnitude using supplied calibration",   model->chisq);
    891             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Magnitude using supplied calibration",   model->nDOF);
     998            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
     999            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
     1000            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
     1001
     1002            // EAM : adding for PV2 outputs:
     1003            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
    8921004
    8931005            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
    … …  
    9251037                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
    9261038                } else {
    927                     psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     1039                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
    9281040                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
    9291041                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
    … …  
    9461058
    9471059                snprintf (name, 64, "EXT_PAR_%02d", k);
    948 
     1060               
    9491061                if (k < model->params->n) {
    9501062                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
    … …  
    9561068            // optionally, write out the covariance matrix values
    9571069            // XXX do I need to pad this to match the biggest covar matrix?
    958             if (model->covar) {
     1070            if (false && model->covar) {
    9591071                for (int iy = 0; iy < model->covar->numCols; iy++) {
    9601072                    for (int ix = iy; ix < model->covar->numCols; ix++) {
    … …  
    10051117        return false;
    10061118    }
     1119    // set up the lookup table to translate between input model types and output model types
     1120    // if not defined it is assumed that the tables are the same
     1121    pmModelClassReadHeader(tableHeader);
    10071122
    10081123    for (long i = 0; i < numSources; i++) {
    … …  
    10131128            return false;
    10141129        }
    1015         // Find the source with this sequence number.
    1016         // XXX: I am assuming that sources is sorted in order of seq.
    10171130        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
    1018         long j = seq < sources->n ? seq : sources->n - 1;
    1019         pmSource *source = NULL;
    1020         for (; j >= 0; j--) {
    1021             source = sources->data[j];
    1022             if (source->seq == seq) {
    1023                 break;
    1024             }
    1025         }
     1131        long j = sourceIndex[seq];
     1132        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
     1133        pmSource *source = sources->data[j];
    10261134        if (!source) {
    10271135            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
    … …  
    10581166        model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
    10591167
     1168        model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ");
     1169        model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF");
     1170        model->flags = psMetadataLookupS16(&status, row, "EXT_FLAGS");
     1171
     1172        // EXT_MODEL_TYPE gives the model chosen by psphot as the best.
     1173        // Putting this into the XFIT table makes 3 copies of it (one for each model)
     1174        // but since we have many fewer XFIT rows than psf rows that is cheaper than putting it
     1175        // in the psf table.
     1176        psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE");
     1177        if (status) {
     1178            // translate between the type value in xfit and values used by this program
     1179            extModelType = pmModelClassGetLocalType(extModelType);
     1180        } else {
     1181            // older cmfs don't have this column
     1182            extModelType = -1;
     1183        }
     1184
    10601185        psEllipseAxes axes;
    10611186        axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
    10621187        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
    10631188        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
    1064         if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
     1189        if (!pmPSF_AxesToModel(PAR, axes, model->class->useReff)) {
    10651190            // Do we need to fail here or can this happen?
    10661191            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
    … …  
    10721197        if (model->params->n > 7) {
    10731198            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
    1074         }
    1075         // read the covariance matrix
    1076         int nparams = model->params->n;
    1077         psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
    1078         for (int y = 0; y < nparams; y++) {
    1079             for (int x = 0; x < nparams; x++) {
    1080                 char name[64];
    1081                 snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
    1082                 covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
     1199            // XXX add an error:
     1200            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
     1201        }
     1202
     1203        // XXX : make this depend on what is in the cmf
     1204        if (0) {
     1205            // read the covariance matrix
     1206            int nparams = model->params->n;
     1207            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
     1208            for (int y = 0; y < nparams; y++) {
     1209                for (int x = 0; x < nparams; x++) {
     1210                    char name[64];
     1211                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
     1212                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
     1213                }
     1214            }
     1215            model->covar = covar;
     1216        }
     1217
     1218        if (modelType == extModelType) {
     1219            // The software that created this source picked this model as the best of the fits.
     1220            // Set the extModel to point to it.
     1221            // This is important for programs like psastro (skycal) so that its output cmfs
     1222            // will have valid EXT_MODEL_TYPE
     1223            psFree(source->modelEXT);
     1224            source->modelEXT = psMemIncrRefCounter(model);
     1225            source->type = PM_SOURCE_TYPE_EXTENDED;
     1226            if (0) {
     1227                // since FLAGS were read we don't need to do this
     1228                source->mode |= PM_SOURCE_MODE_EXTMODEL;
     1229                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    10831230            }
    10841231        }
    1085         model->covar = covar;
    10861232
    10871233        psArrayAdd(source->modelFits, 1, model);
    10881234        psFree(model);
    1089 
    10901235        psFree(row);
    10911236    }
    … …  
    11401285    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    11411286
     1287#ifdef SORT_EXTENSIONS_BY_FLUX
    11421288    // let's write these out in S/N order
    11431289    sources = psArraySort (sources, pmSourceSortByFlux);
     1290#endif
    11441291
    11451292    table = psArrayAllocEmpty (sources->n);
    … …  
    12091356
    12101357        write_annuli:
    1211             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
    1212             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
    1213             psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
    1214             psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
     1358            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
     1359            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
     1360            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
     1361            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
    12151362            psFree (radFlux);
    12161363            psFree (radFluxErr);
    … …  
    12771424            return false;
    12781425        }
    1279         // Find the source with this sequence number.
    1280         // XXX: I am assuming that sources is sorted in order of seq.
    12811426        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
    1282         long j = seq < sources->n ? seq : sources->n - 1;
    1283         pmSource *source = NULL;
    1284         for (; j >= 0; j--) {
    1285             source = sources->data[j];
    1286             if (source->seq == seq) {
    1287                 break;
    1288             }
    1289         }
     1427        long j = sourceIndex[seq];
     1428        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
     1429        pmSource *source = sources->data[j];
    12901430        if (!source) {
    12911431            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
    … …  
    13431483    return true;
    13441484}
     1485
     1486// XXX where should I record the number of columns??
     1487bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
     1488{
     1489    bool status = false;
     1490
     1491    // perform full non-linear fits / extended source analysis?
     1492    if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
     1493        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
     1494        return true;
     1495    }
     1496
     1497    // create a header to hold the output data
     1498    psMetadata *outhead = psMetadataAlloc ();
     1499
     1500    // write the links to the image header
     1501    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname);
     1502
     1503    psF32 Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q");
     1504    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1505
     1506    psF32 NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA");
     1507    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1508
     1509    psF32 clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN");
     1510    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1511
     1512    // They are probably already in this order but ...
     1513    sources = psArraySort (sources, pmSourceSortBySeq);
     1514
     1515    psArray *table = psArrayAllocEmpty (sources->n);
     1516
     1517    for (int i = 0; i < sources->n; i++) {
     1518
     1519        pmSource *thisSource = sources->data[i];
     1520
     1521        // this is the "real" version of this source
     1522        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     1523
     1524        // if we did not fit the galaxy model, modelFits will be NULL
     1525        if (source->modelFits == NULL) continue;
     1526
     1527        // if we did not fit the galaxy model, galaxyFits will also be NULL
     1528        if (source->galaxyFits == NULL) continue;
     1529
     1530        for (int iModel = 0; iModel < source->modelFits->n; iModel++) {
     1531            pmModel *model = source->modelFits->data[iModel];
     1532            if (!model) continue;
     1533
     1534            pmSourceGalaxyFits *galaxyFits = NULL;
     1535            for (int iFit = 0; iFit < source->galaxyFits->n; iFit++) {
     1536                galaxyFits = source->galaxyFits->data[iFit];
     1537                if (model->type == galaxyFits->modelType) break;
     1538                galaxyFits = NULL;
     1539            }
     1540
     1541            if (!galaxyFits) continue;
     1542
     1543            // X,Y coordinates are stored with the model parameters
     1544            psF32 *PAR = model->params->data.F32;
     1545
     1546            psMetadata *row = psMetadataAlloc ();
     1547
     1548            // we write out the x,y positions so people can link to the psf either way (position or ID)
     1549            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
     1550            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",                     galaxyFits->modelType);
     1551            psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
     1552            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
     1553            psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      galaxyFits->nPix);
     1554            // psMetadataAddS32 (row, PS_LIST_TAIL, "REDUCED_TRIALS",   0, "source has reduced number of trials",      galaxyFits->reducedTrials);
     1555
     1556            psVector *Flux = galaxyFits->Flux;
     1557            psVector *dFlux = galaxyFits->dFlux;
     1558            psVector *chisq = galaxyFits->chisq;
     1559
     1560            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
     1561            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
     1562            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
     1563
     1564            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MIN",    0, "fractional major axis min",      galaxyFits->fRmajorMin);
     1565            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MAX",    0, "fractional major axis max",      galaxyFits->fRmajorMax);
     1566            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_DEL",    0, "fractional major axis max",      galaxyFits->fRmajorDel);
     1567            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MIN",    0, "fractional minor axis min",      galaxyFits->fRminorMin);
     1568            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MAX",    0, "fractional minor axis max",      galaxyFits->fRminorMax);
     1569            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_DEL",    0, "fractional minor axis max",      galaxyFits->fRminorDel);
     1570
     1571            psArrayAdd (table, 100, row);
     1572            psFree (row);
     1573        }
     1574    }
     1575
     1576    if (table->n == 0) {
     1577        if (!psFitsWriteBlank (fits, outhead, extname)) {
     1578            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     1579            psFree(outhead);
     1580            psFree(table);
     1581            return false;
     1582        }
     1583        psFree (outhead);
     1584        psFree (table);
     1585        return true;
     1586    }
     1587
     1588    psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname);
     1589    if (!psFitsWriteTable (fits, outhead, table, extname)) {
     1590        psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname);
     1591        psFree (outhead);
     1592        psFree(table);
     1593        return false;
     1594    }
     1595    psFree (outhead);
     1596    psFree (table);
     1597    return true;
     1598}
     1599
     1600bool pmSourcesRead_CMF_@CMFMODE@_XGAL(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, psArray *sources, long *sourceIndex)
     1601{
     1602    PS_ASSERT_PTR_NON_NULL(fits, false);
     1603    PS_ASSERT_PTR_NON_NULL(sources, false);
     1604
     1605    bool status;
     1606    long numSources = psFitsTableSize(fits); // Number of sources in table
     1607    if (numSources == 0) {
     1608        psError(psErrorCodeLast(), false, "XGAL Table contains no entries\n");
     1609        return false;
     1610    }
     1611    if (!readout->analysis) {
     1612        readout->analysis = psMetadataAlloc();
     1613    }
     1614    psF32 Q = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_Q");
     1615    // XXX: turn this into an assert once we're done
     1616    if (status) {
     1617        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1618
     1619        psF32 NSigma = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_NSIGMA");
     1620        psAssert(status, "missing GALAXY_SHAPES_NSIGMA");
     1621        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1622
     1623        psF32 clampSN = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_CLAMP_SN");
     1624        psAssert(status, "missing GALAXY_SHAPES_CLAMP_SN");
     1625        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1626    }
     1627
     1628    for (long i = 0; i < numSources; i++) {
     1629        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     1630        if (!row) {
     1631            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     1632            psFree(row);
     1633            return false;
     1634        }
     1635        // Find the source with this sequence number.
     1636        // XXX: I am assuming that sources is sorted in order of seq
     1637        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     1638        long j = sourceIndex[seq];
     1639        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
     1640
     1641        pmSource *source = sources->data[j];
     1642        if (!source) {
     1643            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     1644            psFree(row);
     1645            return false;
     1646        }
     1647
     1648        int modelType = psMetadataLookupS32 (&status,    row, "MODEL_TYPE");
     1649        psVector *Flux  = psMetadataLookupVector(&status, row, "GAL_FLUX");
     1650        psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR");
     1651        psVector *chisq = psMetadataLookupVector(&status, row, "GAL_CHISQ");
     1652
     1653        if (Flux && Flux->n > 0) {
     1654            if (!source->galaxyFits) {
     1655                source->galaxyFits = psArrayAllocEmpty(1);
     1656            }
     1657
     1658            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     1659
     1660            psArrayAdd(source->galaxyFits, 1, galaxyFits);
     1661
     1662            psFree(galaxyFits);
     1663            galaxyFits->modelType = modelType;
     1664            galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
     1665
     1666            galaxyFits->fRmajorMin = psMetadataLookupF32(&status, row, "FR_MAJOR_MIN");
     1667            galaxyFits->fRmajorMax = psMetadataLookupF32(&status, row, "FR_MAJOR_MAX");
     1668            galaxyFits->fRmajorDel = psMetadataLookupF32(&status, row, "FR_MAJOR_DEL");
     1669            galaxyFits->fRminorMin = psMetadataLookupF32(&status, row, "FR_MINOR_MIN");
     1670            galaxyFits->fRminorMax = psMetadataLookupF32(&status, row, "FR_MINOR_MAX");
     1671            galaxyFits->fRminorDel = psMetadataLookupF32(&status, row, "FR_MINOR_DEL");
     1672
     1673            psFree(galaxyFits->Flux);
     1674            galaxyFits->Flux  = psMemIncrRefCounter(Flux);
     1675            psFree(galaxyFits->dFlux);
     1676            galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
     1677            psFree(galaxyFits->chisq);
     1678            galaxyFits->chisq = psMemIncrRefCounter(chisq);
     1679        }
     1680
     1681        psFree(row);
     1682    }
     1683
     1684    return true;
     1685}
Note: See TracChangeset for help on using the changeset viewer.