IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36757 for trunk/psModules


Ignore:
Timestamp:
May 21, 2014, 9:37:09 AM (12 years ago)
Author:
bills
Message:

Changes to full force galaxy shapes analysis

  • Put parameters for all extended models into cff file
  • in psphotFullForce select models to force and those to skip
  • skip sources that are very likely to be stars
  • clean up a few "hacky" pieces of the implementation
Location:
trunk/psModules/src/objects
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r36623 r36757  
    12021202        break;
    12031203
    1204       case PM_FPA_FILE_CFF:
     1204      case PM_FPA_FILE_CFF: {
     1205        // determine the output table format
     1206        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
     1207        if (!status) {
     1208            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
     1209            return false;
     1210        }
    12051211        // read in header, if not yet loaded
    12061212        hdu = pmFPAviewThisHDU (view, file->fpa);
    … …  
    12441250        }
    12451251
    1246         sources = pmSourcesRead_CFF(file->fits, hdu->header);
     1252        sources = pmSourcesRead_CFF(file->fits, hdu->header, recipe);
    12471253
    12481254        psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname);
    … …  
    12501256        psFree (dataname);
    12511257        psFree (tableHeader);
     1258        }
    12521259        break;
    12531260

  • trunk/psModules/src/objects/pmSourceIO.h

    r36623 r36757  
    5555
    5656psArray *pmSourcesReadCMP (char *filename, psMetadata *header);
    57 psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header);
     57psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe);
    5858bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe);
    5959

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

    r36747 r36757  
    5454
    5555// read in sources readout from a cff fits file
    56 psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header)
     56psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe)
    5757{
    5858    PS_ASSERT_PTR_NON_NULL(fits, false);
    … …  
    7474    }
    7575    PS_ASSERT_INT_NONNEGATIVE(modelType, NULL);
     76
     77    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
     78    pmModelType devModelType    = pmModelClassGetType("PS_MODEL_DEV");
     79
     80    psString modelForce = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE");
     81    psF32 forceDevSersicMin = NAN;
     82    bool forceAll = false;
     83    int modelTypeForce = 0;
     84    if (!strcmp(modelForce, "ALL")) {
     85        forceAll = true;
     86    } else if (!strcmp(modelForce, "PS_MODEL_SERSIC")) {
     87        modelTypeForce = sersicModelType;
     88        forceDevSersicMin = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN");
     89    } else {
     90        modelTypeForce = pmModelClassGetType(modelForce);
     91        PS_ASSERT_INT_NONNEGATIVE(modelTypeForce, NULL);
     92    }
     93
     94    // skip extended model types for likely stars
     95    // max value of KronMag - psfMag to keep ...
     96    psF32 starCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_MAGDIFF_MAX");
     97    if (!status) {
     98        starCut = 0;
     99    }
     100    // ... unless SN is less than this value
     101    psF32 SNMinForCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_CUT_SN_MIN");
     102    if (!status) {
     103        SNMinForCut = 10;
     104    }
    76105
    77106    // We get the size of the table, and allocate the array of sources first because the table
    … …  
    116145        float Rminor     = psMetadataLookupF32 (&status, row, "R_MINOR");
    117146        float theta      = psMetadataLookupF32 (&status, row, "THETA");
     147        float chisq      = psMetadataLookupF32 (&status, row, "CHISQ");
     148        float nDOF       = psMetadataLookupF32 (&status, row, "NDOF");
     149        float magDiff    = psMetadataLookupF32 (&status, row, "MAG_DIFF");
     150        psS16 modelFlags = psMetadataLookupS32 (&status, row, "MODEL_FLAGS");
    118151
    119152        int   galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
    120         if (status) {
     153        if (status && galaxyModelType >= 0) {
    121154            galaxyModelType = pmModelClassGetLocalType(galaxyModelType);
    122155        } else {
    … …  
    201234
    202235        }
     236        bool saveExtModelParams = false;
    203237        if (fitGalaxy && galaxyModelType >= 0) {
     238            // skip likely stars
     239            if (magDiff < starCut || SN < SNMinForCut) {
     240                if (forceAll) {
     241                    saveExtModelParams = true;
     242                } else if (galaxyModelType == modelTypeForce) {
     243                    // This is the model type that we are looking for
     244                    // proceed
     245                    saveExtModelParams = true;
     246                } else if (modelTypeForce == sersicModelType && galaxyModelType == devModelType) {
     247                    // We're doing sersic models, if sersic index is greater than the recipe's minimum index
     248                    // do dev model as well
     249                    if (isfinite(forceDevSersicMin) || Sindex >= forceDevSersicMin) {
     250                        saveExtModelParams = true;
     251                    }
     252                } else {
     253                    // not interested in this model
     254                }
     255            }
     256        }
     257
     258        if (saveExtModelParams) {
    204259            if (!source->modelFits) {
    205260                source->modelFits = psArrayAllocEmpty (1);
    … …  
    223278            }
    224279
     280            model->chisq = chisq;
     281            model->nDOF = nDOF;
     282            model->flags = modelFlags;
     283
    225284            psArrayAdd (source->modelFits, 1, model);
    226285
    … …  
    288347        psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC];
    289348        psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC];
     349        psF32 chisq[MAX_ROWS_PER_SRC], nDOF[MAX_ROWS_PER_SRC];
     350        psS32 modelFlags[MAX_ROWS_PER_SRC];
    290351        psS32 modelType[MAX_ROWS_PER_SRC];
    291         psF32 sersicIndex[MAX_ROWS_PER_SRC];
     352        psF32 sersicIndex = NAN;
    292353        bool fitGalaxy = false;
    293354        bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false;
    294355        int n_rows = 0;
    295356
    296         psF32 SN = source->moments->KronFlux/source->moments->KronFluxErr;
     357        psF32 kronFlux = source->moments->KronFlux;
     358        psF32 SN = NAN;
     359        psF32 magDiff = NAN;
     360        if (isfinite(kronFlux) && isfinite(source->moments->KronFluxErr) && isfinite(source->psfMag)) {
     361            SN = kronFlux/source->moments->KronFluxErr;
     362            // kronMag - psfMag for use as star/glaxy separator
     363            magDiff = -2.5 * log10(kronFlux) - source->psfMag ;
     364        }
    297365
    298366        // start with psf model
    … …  
    312380        theta[0] = 0;
    313381        modelType[0] = -1;
    314         sersicIndex[0] = NAN;
     382        sersicIndex = NAN;
     383        chisq[0] = NAN;
     384        nDOF[0] = NAN;
     385        modelFlags[0] = 0;
    315386
    316387        if (source->modelFits != NULL) {
    … …  
    338409                        flux[n_rows] = pow(10.0, -0.4*model->mag);
    339410                        modelType[n_rows] = model->type;
     411                        chisq[n_rows] = model->chisq;
     412                        nDOF[n_rows] = model->nDOF;
     413                        modelFlags[n_rows] = model->flags;
    340414                        fitGalaxy = true;
    341415                        if (model->params->n == 8) {
    342416                            // this will save the sersic index in all models but that might
    343417                            // be useful
    344                             sersicIndex[n_rows] = 0.5 / PAR[PM_PAR_7];
     418                            sersicIndex = 0.5 / PAR[PM_PAR_7];
    345419                        } else {
    346                             sersicIndex[n_rows] = NAN;
     420                            sersicIndex = NAN;
    347421                        }
    348422
    … …  
    393467                        flux[0] = pow(10.0, -0.4*model->mag);
    394468                        modelType[0] = model->type;
     469                        chisq[0] = model->chisq;
     470                        nDOF[0] = model->nDOF;
     471                        modelFlags[0] = model->flags;
    395472                        fitGalaxy = true;
    396473                        if (model->type == sersicModelType) {
    397474                            PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false);
    398                             sersicIndex[0] = 0.5 / PAR[PM_PAR_7];
     475                            sersicIndex = 0.5 / PAR[PM_PAR_7];
    399476                        }
    400477
    … …  
    404481                        // sersic index above the recipe limit.
    405482                        if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) &&
    406                             sersicIndex[0] > sersicMinDev) {
     483                            sersicIndex > sersicMinDev) {
    407484
    408485                            model = source->modelFits->data[jModelDev];
    … …  
    419496                                flux[1] = pow(10.0, -0.4*model->mag);
    420497                                modelType[1] = model->type;
    421                                 sersicIndex[1] = NAN;
     498                                sersicIndex = NAN;
     499                                chisq[1] = model->chisq;
     500                                nDOF[1] = model->nDOF;
     501                                modelFlags[1] = model->flags;
    422502                                n_rows = 2;
    423503                            }
    … …  
    439519            psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX",             0, "flux per second",       flux[j]/exptime);
    440520            psMetadataAddF32 (row, PS_LIST_TAIL, "SN",               0, "kron flux signal to noise", SN);
     521            psMetadataAddF32 (row, PS_LIST_TAIL, "MAG_DIFF",         0, "psf mag - kron mag",    magDiff);
    441522            psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS",        0, "aperture radius",       source->apRadius);
    442523            psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS",      0, "Kron radius",           source->moments->Mrf * 2.5);
    … …  
    448529            psMetadataAddF32 (row, PS_LIST_TAIL, "THETA",            0, "theta",                 theta[j]);
    449530            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "model type",            modelType[j]);
    450             psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex[j]);
     531            psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX",            0, "sersic index",          sersicIndex);
     532            psMetadataAddF32 (row, PS_LIST_TAIL, "CHISQ",            0, "chisq",                 chisq[j]);
     533            psMetadataAddF32 (row, PS_LIST_TAIL, "NDOF",             0, "n degrees of freedom",  nDOF[j]);
     534            psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_FLAGS",      0, "model flags",           modelFlags[j]);
    451535
    452536            psArrayAdd(table, 100, row);
Note: See TracChangeset for help on using the changeset viewer.