IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 19, 2014, 4:05:27 PM (12 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ps2-tc3-20130727
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ps2-tc3-20130727

  • branches/eam_branches/ps2-tc3-20130727/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1515pmSourceIO_CMF_PS1_SV1.c
        1616pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
         18pmSourceIO_CMF_PS1_SV4.c
        1719pmSourceIO_CMF_PS1_DV1.c
        1820pmSourceIO_CMF_PS1_DV2.c
        1921pmSourceIO_CMF_PS1_DV3.c
         22pmSourceIO_CMF_PS1_DV4.c
        2023

  • branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMF.c.in

    r36680 r37403  
    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"
    … …  
    142142        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    143143        @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    144         @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);
    145145        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
    146146        @>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);
    147148
    148149        @<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);
    … …  
    163164        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    164165
     166        // PSF shape parameters:
    165167        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
    166168        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
    167169        @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
    168177        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    169178        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    … …  
    197206        }
    198207
    199         if (source->lensingOBJ && source->lensingPSF->smear) {
     208        if (source->lensingPSF && source->lensingPSF->smear) {
    200209          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
    201210          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
    … …  
    205214        }
    206215
    207         if (source->lensingOBJ && source->lensingPSF->shear) {
     216        if (source->lensingPSF && source->lensingPSF->shear) {
    208217          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
    209218          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
    … …  
    212221          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2);
    213222        }
     223
     224        // if lensing params exist also include the backmapped chipID and chip coordinates
     225        if (source->lensingPSF && source->lensingPSF->shear) {
     226            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",   PS_DATA_S16, "id of warp input chip",     source->chipNum);
     227            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",    PS_DATA_S16, "x coord in warp input chip",     source->chipX);
     228            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y",    PS_DATA_S16, "y coord in warp input chip",     source->chipY);
     229            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "PADDING3",      PS_DATA_S16, "more padding", 0);
     230        }
    214231
    215232        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
    … …  
    302319
    303320    // define PSF model type
    304     // XXX need to carry the extra model parameters
    305321    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
    306322
    307     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     323    // if header does not define the model, default to a gaussian
     324    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
    308325    if (PSF_NAME != NULL) {
    309326        modelType = pmModelClassGetType (PSF_NAME);
    310327    }
    311328    assert (modelType > -1);
     329
     330    // do we expect to find lensing parameters?
     331    bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ");
     332    bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF");
    312333
    313334    // We get the size of the table, and allocate the array of sources first because the table
    … …  
    344365        @ALL@     axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");
    345366        @ALL@     axes.theta        = axes.theta * PS_RAD_DEG;
     367       
     368        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
     369        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
     370        @>PS1_V4,>PS1_SV2,>PS1_DV3@ }
    346371
    347372        @ALL@     PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
    … …  
    355380        @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
    356381        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
     382        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX");
     383        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
    357384
    358385        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
    … …  
    365392        @ALL@     dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
    366393
    367         pmPSF_AxesToModel (PAR, axes, modelType);
     394        pmPSF_AxesToModel (PAR, axes, model->class->useReff);
    368395
    369396        @ALL@     float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
    … …  
    383410        @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    384411        @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    385         @ALL@     source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
     412        @ALL@     source->apRadius  = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS");
     413        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
    386414
    387415        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    410438        @>PS1_V2,PS1_SV?@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
    411439        @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0;
     440
     441        // Lensing parameters (on read if PS1_V5+)
     442        if (haveLensOBJ) {
     443          source->lensingOBJ = pmSourceLensingAlloc ();
     444          source->lensingOBJ->smear = pmLensingParsAlloc();
     445          source->lensingOBJ->shear = pmLensingParsAlloc();
     446
     447          @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
     448          @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
     449          @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
     450          @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
     451          @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
     452          @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
     453          @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
     454          @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
     455          @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
     456          @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
     457        }
     458
     459        if (haveLensPSF) {
     460          source->lensingPSF = pmSourceLensingAlloc ();
     461          source->lensingPSF->smear = pmLensingParsAlloc();
     462          source->lensingPSF->shear = pmLensingParsAlloc();
     463
     464          @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
     465          @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
     466          @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
     467          @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
     468          @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
     469          @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
     470          @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
     471          @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
     472          @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
     473          @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
     474        }
    412475
    413476        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
    … …  
    482545    }
    483546
     547#ifdef SORT_EXTENSIONS_BY_FLUX
    484548    // let's write these out in S/N order
    485549    sources = psArraySort (sources, pmSourceSortByFlux);
     550#endif
    486551
    487552    table = psArrayAllocEmpty (sources->n);
    … …  
    563628                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point
    564629                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
     630                // XXX NOTE EAM 20140806 : PETRO_MAG_ERR is inverted!! oops!!
    565631                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
    566632                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
    … …  
    821887    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
    822888
     889#ifdef SORT_EXTENSIONS_BY_FLUX
    823890    // let's write these out in S/N order
    824891    sources = psArraySort (sources, pmSourceSortByFlux);
     892#endif
    825893
    826894    float magOffset;
    … …  
    847915    }
    848916
    849     @>PS1_DV2@ pmChip *chip = readout->parent->parent;
     917    @>PS1_DV2,>PS1_SV3@ pmChip *chip = readout->parent->parent;
    850918
    851919    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
    … …  
    906974            }
    907975
    908             @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    909             @>PS1_DV2@ float posAngle = 0.0;
    910             @>PS1_DV2@ float pltScale = 0.0;
    911             @>PS1_DV2@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
    912             @>PS1_DV2@ double raPos = ptSky.r*PS_DEG_RAD;
    913             @>PS1_DV2@ double decPos = ptSky.d*PS_DEG_RAD;
    914             @>PS1_DV2@ posAngle *= PS_DEG_RAD;
    915             @>PS1_DV2@ pltScale *= PS_DEG_RAD*3600.0;
     976            @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
     977            @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0;
     978            @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0;
     979            @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     980            @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD;
     981            @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD;
     982            @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD;
     983            @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0;
    916984
    917985            float kronFlux = source->moments ? source->moments->KronFlux : NAN;
    … …  
    927995            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
    928996            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
    929             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model ra coordinate",                    raPos);
    930             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model dec coordinate",                   decPos);
     997            @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model ra coordinate",                    raPos);
     998            @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model dec coordinate",                   decPos);
    931999            @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
    9321000            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
    … …  
    9421010
    9431011            // EAM : adding for PV2 outputs:
    944             @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0);
     1012            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
    9451013
    9461014            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
    … …  
    9781046                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
    9791047                } else {
    980                     psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     1048                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
    9811049                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
    9821050                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
    … …  
    11091177        model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ");
    11101178        model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF");
     1179        model->flags = psMetadataLookupS16(&status, row, "EXT_FLAGS");
    11111180
    11121181        // EXT_MODEL_TYPE gives the model chosen by psphot as the best.
    … …  
    11271196        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
    11281197        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
    1129         if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
     1198        if (!pmPSF_AxesToModel(PAR, axes, model->class->useReff)) {
    11301199            // Do we need to fail here or can this happen?
    11311200            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
    … …  
    12251294    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
    12261295
     1296#ifdef SORT_EXTENSIONS_BY_FLUX
    12271297    // let's write these out in S/N order
    12281298    sources = psArraySort (sources, pmSourceSortByFlux);
     1299#endif
    12291300
    12301301    table = psArrayAllocEmpty (sources->n);
    … …  
    14391510    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname);
    14401511
    1441     psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE");
    1442 
    1443     // let's write these out in S/N order
    1444     sources = psArraySort (sources, pmSourceSortByFlux);
     1512    psF32 Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q");
     1513    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1514
     1515    psF32 NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA");
     1516    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1517
     1518    psF32 clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN");
     1519    psMetadataAddF32 (outhead, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1520
     1521    // They are probably already in this order but ...
     1522    sources = psArraySort (sources, pmSourceSortBySeq);
    14451523
    14461524    psArray *table = psArrayAllocEmpty (sources->n);
    … …  
    14591537        if (source->galaxyFits == NULL) continue;
    14601538
    1461         pmModel *model = source->modelFits->data[0];
    1462         if (!model) return false;
    1463 
    1464         // X,Y coordinates are stored with the model parameters
    1465         psF32 *PAR = model->params->data.F32;
    1466 
    1467         psMetadata *row = psMetadataAlloc ();
    1468 
    1469         // we write out the x,y positions so people can link to the psf either way (position or ID)
    1470         psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
    1471         psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
    1472         psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
    1473         psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      source->galaxyFits->nPix);
    1474 
    1475         psVector *Flux = source->galaxyFits->Flux;
    1476         psVector *dFlux = source->galaxyFits->dFlux;
    1477         psVector *chisq = source->galaxyFits->chisq;
    1478 
    1479         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
    1480         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
    1481         psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
    1482 
    1483         psArrayAdd (table, 100, row);
    1484         psFree (row);
     1539        for (int iModel = 0; iModel < source->modelFits->n; iModel++) {
     1540            pmModel *model = source->modelFits->data[iModel];
     1541            if (!model) continue;
     1542
     1543            pmSourceGalaxyFits *galaxyFits = NULL;
     1544            for (int iFit = 0; iFit < source->galaxyFits->n; iFit++) {
     1545                galaxyFits = source->galaxyFits->data[iFit];
     1546                if (model->type == galaxyFits->modelType) break;
     1547                galaxyFits = NULL;
     1548            }
     1549
     1550            if (!galaxyFits) continue;
     1551
     1552            // X,Y coordinates are stored with the model parameters
     1553            psF32 *PAR = model->params->data.F32;
     1554
     1555            psMetadata *row = psMetadataAlloc ();
     1556
     1557            // we write out the x,y positions so people can link to the psf either way (position or ID)
     1558            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
     1559            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",                     galaxyFits->modelType);
     1560            psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
     1561            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
     1562            psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      galaxyFits->nPix);
     1563            // psMetadataAddS32 (row, PS_LIST_TAIL, "REDUCED_TRIALS",   0, "source has reduced number of trials",      galaxyFits->reducedTrials);
     1564
     1565            psVector *Flux = galaxyFits->Flux;
     1566            psVector *dFlux = galaxyFits->dFlux;
     1567            psVector *chisq = galaxyFits->chisq;
     1568
     1569            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
     1570            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
     1571            psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
     1572
     1573            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MIN",    0, "fractional major axis min",      galaxyFits->fRmajorMin);
     1574            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_MAX",    0, "fractional major axis max",      galaxyFits->fRmajorMax);
     1575            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MAJOR_DEL",    0, "fractional major axis max",      galaxyFits->fRmajorDel);
     1576            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MIN",    0, "fractional minor axis min",      galaxyFits->fRminorMin);
     1577            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_MAX",    0, "fractional minor axis max",      galaxyFits->fRminorMax);
     1578            psMetadataAddF32 (row, PS_LIST_TAIL, "FR_MINOR_DEL",    0, "fractional minor axis max",      galaxyFits->fRminorDel);
     1579
     1580            psArrayAdd (table, 100, row);
     1581            psFree (row);
     1582        }
    14851583    }
    14861584
    … …  
    15201618        return false;
    15211619    }
     1620    if (!readout->analysis) {
     1621        readout->analysis = psMetadataAlloc();
     1622    }
     1623    psF32 Q = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_Q");
     1624    // XXX: turn this into an assert once we're done
     1625    if (status) {
     1626        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_Q", PS_META_REPLACE, "", Q);
     1627
     1628        psF32 NSigma = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_NSIGMA");
     1629        psAssert(status, "missing GALAXY_SHAPES_NSIGMA");
     1630        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_NSIGMA", PS_META_REPLACE, "", NSigma);
     1631
     1632        psF32 clampSN = psMetadataLookupF32(&status, tableHeader, "GALAXY_SHAPES_CLAMP_SN");
     1633        psAssert(status, "missing GALAXY_SHAPES_CLAMP_SN");
     1634        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "GALAXY_SHAPES_CLAMP_SN", PS_META_REPLACE, "", clampSN);
     1635    }
    15221636
    15231637    for (long i = 0; i < numSources; i++) {
    … …  
    15411655        }
    15421656
     1657        int modelType = psMetadataLookupS32 (&status,    row, "MODEL_TYPE");
    15431658        psVector *Flux  = psMetadataLookupVector(&status, row, "GAL_FLUX");
    15441659        psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR");
    … …  
    15461661
    15471662        if (Flux && Flux->n > 0) {
    1548             psFree(source->galaxyFits);
    1549             source->galaxyFits = pmSourceGalaxyFitsAlloc();
    1550             source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
    1551 
    1552             psFree(source->galaxyFits->Flux);
    1553             source->galaxyFits->Flux  = psMemIncrRefCounter(Flux);
    1554             psFree(source->galaxyFits->dFlux);
    1555             source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
    1556             psFree(source->galaxyFits->chisq);
    1557             source->galaxyFits->chisq = psMemIncrRefCounter(chisq);
     1663            if (!source->galaxyFits) {
     1664                source->galaxyFits = psArrayAllocEmpty(1);
     1665            }
     1666
     1667            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
     1668
     1669            psArrayAdd(source->galaxyFits, 1, galaxyFits);
     1670
     1671            psFree(galaxyFits);
     1672            galaxyFits->modelType = modelType;
     1673            galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
     1674
     1675            galaxyFits->fRmajorMin = psMetadataLookupF32(&status, row, "FR_MAJOR_MIN");
     1676            galaxyFits->fRmajorMax = psMetadataLookupF32(&status, row, "FR_MAJOR_MAX");
     1677            galaxyFits->fRmajorDel = psMetadataLookupF32(&status, row, "FR_MAJOR_DEL");
     1678            galaxyFits->fRminorMin = psMetadataLookupF32(&status, row, "FR_MINOR_MIN");
     1679            galaxyFits->fRminorMax = psMetadataLookupF32(&status, row, "FR_MINOR_MAX");
     1680            galaxyFits->fRminorDel = psMetadataLookupF32(&status, row, "FR_MINOR_DEL");
     1681
     1682            psFree(galaxyFits->Flux);
     1683            galaxyFits->Flux  = psMemIncrRefCounter(Flux);
     1684            psFree(galaxyFits->dFlux);
     1685            galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
     1686            psFree(galaxyFits->chisq);
     1687            galaxyFits->chisq = psMemIncrRefCounter(chisq);
    15581688        }
    15591689
Note: See TracChangeset for help on using the changeset viewer.