Index: /trunk/ippconfig/recipes/psphot.config =================================================================== --- /trunk/ippconfig/recipes/psphot.config (revision 36756) +++ /trunk/ippconfig/recipes/psphot.config (revision 36757) @@ -198,9 +198,18 @@ # Extended Model to use for forced photometry in full force # In addtion to PS_MODEL_SERSIC, PS_MODEL_DEV, and PS_MODEL_EXP the values BEST and ALL may be selected -EXT_MODEL_TYPE_FOR_CFF STR PS_MODEL_SERSIC +# This specifies what gets saved in cff file output by psastro +EXT_MODEL_TYPE_FOR_CFF STR ALL +# this model is forced +EXT_MODEL_TYPE_FORCE STR PS_MODEL_SERSIC # if the sersic model is selected, the dev model will also be forced for sources with sersic index above this limit # (Note: this is not used if BEST was selected) EXT_MODEL_FORCE_DEV_SERSIC_MIN F32 3.0 + +# to avoid forcing extended models for likely stars +# only force sources with kron mag - psf mag less than this value ... +EXT_MODEL_FORCE_MAGDIFF_MAX F32 0.0 +# ... unless SN is less than this value +EXT_MODEL_FORCE_CUT_SN_MIN F32 10.0 # Limits on extended source fits by galactic coordinates Index: /trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.c (revision 36756) +++ /trunk/psModules/src/objects/pmSourceIO.c (revision 36757) @@ -1202,5 +1202,11 @@ break; - case PM_FPA_FILE_CFF: + case PM_FPA_FILE_CFF: { + // determine the output table format + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } // read in header, if not yet loaded hdu = pmFPAviewThisHDU (view, file->fpa); @@ -1244,5 +1250,5 @@ } - sources = pmSourcesRead_CFF(file->fits, hdu->header); + sources = pmSourcesRead_CFF(file->fits, hdu->header, recipe); psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname); @@ -1250,4 +1256,5 @@ psFree (dataname); psFree (tableHeader); + } break; Index: /trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.h (revision 36756) +++ /trunk/psModules/src/objects/pmSourceIO.h (revision 36757) @@ -55,5 +55,5 @@ psArray *pmSourcesReadCMP (char *filename, psMetadata *header); -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header); +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe); bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe); Index: /trunk/psModules/src/objects/pmSourceIO_CFF.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36756) +++ /trunk/psModules/src/objects/pmSourceIO_CFF.c (revision 36757) @@ -54,5 +54,5 @@ // read in sources readout from a cff fits file -psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header) +psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header, psMetadata *recipe) { PS_ASSERT_PTR_NON_NULL(fits, false); @@ -74,4 +74,33 @@ } PS_ASSERT_INT_NONNEGATIVE(modelType, NULL); + + pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModelType devModelType = pmModelClassGetType("PS_MODEL_DEV"); + + psString modelForce = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + psF32 forceDevSersicMin = NAN; + bool forceAll = false; + int modelTypeForce = 0; + if (!strcmp(modelForce, "ALL")) { + forceAll = true; + } else if (!strcmp(modelForce, "PS_MODEL_SERSIC")) { + modelTypeForce = sersicModelType; + forceDevSersicMin = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_DEV_SERSIC_MIN"); + } else { + modelTypeForce = pmModelClassGetType(modelForce); + PS_ASSERT_INT_NONNEGATIVE(modelTypeForce, NULL); + } + + // skip extended model types for likely stars + // max value of KronMag - psfMag to keep ... + psF32 starCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_MAGDIFF_MAX"); + if (!status) { + starCut = 0; + } + // ... unless SN is less than this value + psF32 SNMinForCut = psMetadataLookupF32(&status, recipe, "EXT_MODEL_FORCE_CUT_SN_MIN"); + if (!status) { + SNMinForCut = 10; + } // We get the size of the table, and allocate the array of sources first because the table @@ -116,7 +145,11 @@ float Rminor = psMetadataLookupF32 (&status, row, "R_MINOR"); float theta = psMetadataLookupF32 (&status, row, "THETA"); + float chisq = psMetadataLookupF32 (&status, row, "CHISQ"); + float nDOF = psMetadataLookupF32 (&status, row, "NDOF"); + float magDiff = psMetadataLookupF32 (&status, row, "MAG_DIFF"); + psS16 modelFlags = psMetadataLookupS32 (&status, row, "MODEL_FLAGS"); int galaxyModelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); - if (status) { + if (status && galaxyModelType >= 0) { galaxyModelType = pmModelClassGetLocalType(galaxyModelType); } else { @@ -201,5 +234,27 @@ } + bool saveExtModelParams = false; if (fitGalaxy && galaxyModelType >= 0) { + // skip likely stars + if (magDiff < starCut || SN < SNMinForCut) { + if (forceAll) { + saveExtModelParams = true; + } else if (galaxyModelType == modelTypeForce) { + // This is the model type that we are looking for + // proceed + saveExtModelParams = true; + } else if (modelTypeForce == sersicModelType && galaxyModelType == devModelType) { + // We're doing sersic models, if sersic index is greater than the recipe's minimum index + // do dev model as well + if (isfinite(forceDevSersicMin) || Sindex >= forceDevSersicMin) { + saveExtModelParams = true; + } + } else { + // not interested in this model + } + } + } + + if (saveExtModelParams) { if (!source->modelFits) { source->modelFits = psArrayAllocEmpty (1); @@ -223,4 +278,8 @@ } + model->chisq = chisq; + model->nDOF = nDOF; + model->flags = modelFlags; + psArrayAdd (source->modelFits, 1, model); @@ -288,11 +347,20 @@ psF32 xPos[MAX_ROWS_PER_SRC], yPos[MAX_ROWS_PER_SRC], flux[MAX_ROWS_PER_SRC]; psF32 rMajor[MAX_ROWS_PER_SRC], rMinor[MAX_ROWS_PER_SRC], theta[MAX_ROWS_PER_SRC]; + psF32 chisq[MAX_ROWS_PER_SRC], nDOF[MAX_ROWS_PER_SRC]; + psS32 modelFlags[MAX_ROWS_PER_SRC]; psS32 modelType[MAX_ROWS_PER_SRC]; - psF32 sersicIndex[MAX_ROWS_PER_SRC]; + psF32 sersicIndex = NAN; bool fitGalaxy = false; bool psfStar = (source->mode & PM_SOURCE_MODE_PSFSTAR) ? true : false; int n_rows = 0; - psF32 SN = source->moments->KronFlux/source->moments->KronFluxErr; + psF32 kronFlux = source->moments->KronFlux; + psF32 SN = NAN; + psF32 magDiff = NAN; + if (isfinite(kronFlux) && isfinite(source->moments->KronFluxErr) && isfinite(source->psfMag)) { + SN = kronFlux/source->moments->KronFluxErr; + // kronMag - psfMag for use as star/glaxy separator + magDiff = -2.5 * log10(kronFlux) - source->psfMag ; + } // start with psf model @@ -312,5 +380,8 @@ theta[0] = 0; modelType[0] = -1; - sersicIndex[0] = NAN; + sersicIndex = NAN; + chisq[0] = NAN; + nDOF[0] = NAN; + modelFlags[0] = 0; if (source->modelFits != NULL) { @@ -338,11 +409,14 @@ flux[n_rows] = pow(10.0, -0.4*model->mag); modelType[n_rows] = model->type; + chisq[n_rows] = model->chisq; + nDOF[n_rows] = model->nDOF; + modelFlags[n_rows] = model->flags; fitGalaxy = true; if (model->params->n == 8) { // this will save the sersic index in all models but that might // be useful - sersicIndex[n_rows] = 0.5 / PAR[PM_PAR_7]; + sersicIndex = 0.5 / PAR[PM_PAR_7]; } else { - sersicIndex[n_rows] = NAN; + sersicIndex = NAN; } @@ -393,8 +467,11 @@ flux[0] = pow(10.0, -0.4*model->mag); modelType[0] = model->type; + chisq[0] = model->chisq; + nDOF[0] = model->nDOF; + modelFlags[0] = model->flags; fitGalaxy = true; if (model->type == sersicModelType) { PS_ASSERT_FLOAT_LARGER_THAN(PAR[PM_PAR_7], 0.0, false); - sersicIndex[0] = 0.5 / PAR[PM_PAR_7]; + sersicIndex = 0.5 / PAR[PM_PAR_7]; } @@ -404,5 +481,5 @@ // sersic index above the recipe limit. if (jModelSelected == -1 && jModelDev >= 0 && isfinite(sersicMinDev) && - sersicIndex[0] > sersicMinDev) { + sersicIndex > sersicMinDev) { model = source->modelFits->data[jModelDev]; @@ -419,5 +496,8 @@ flux[1] = pow(10.0, -0.4*model->mag); modelType[1] = model->type; - sersicIndex[1] = NAN; + sersicIndex = NAN; + chisq[1] = model->chisq; + nDOF[1] = model->nDOF; + modelFlags[1] = model->flags; n_rows = 2; } @@ -439,4 +519,5 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "FLUX", 0, "flux per second", flux[j]/exptime); psMetadataAddF32 (row, PS_LIST_TAIL, "SN", 0, "kron flux signal to noise", SN); + psMetadataAddF32 (row, PS_LIST_TAIL, "MAG_DIFF", 0, "psf mag - kron mag", magDiff); psMetadataAddF32 (row, PS_LIST_TAIL, "AP_RADIUS", 0, "aperture radius", source->apRadius); psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_RADIUS", 0, "Kron radius", source->moments->Mrf * 2.5); @@ -448,5 +529,8 @@ psMetadataAddF32 (row, PS_LIST_TAIL, "THETA", 0, "theta", theta[j]); psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "model type", modelType[j]); - psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "INDEX", 0, "sersic index", sersicIndex); + psMetadataAddF32 (row, PS_LIST_TAIL, "CHISQ", 0, "chisq", chisq[j]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NDOF", 0, "n degrees of freedom", nDOF[j]); + psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_FLAGS", 0, "model flags", modelFlags[j]); psArrayAdd(table, 100, row); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 36756) +++ /trunk/psphot/src/psphot.h (revision 36757) @@ -435,15 +435,5 @@ float NSigma; float clampSN; -#ifdef notdef - float fRmajorMin; - float fRmajorMax; - float fRmajorDel; - float fRminorMin; - float fRminorMax; - float fRminorDel; - float snMinAllTrials; - float reduceTrialsFactor; - bool reducedTrials; -#endif + int extModelType; } psphotGalaxyShapeOptions; Index: /trunk/psphot/src/psphotGalaxyShape.c =================================================================== --- /trunk/psphot/src/psphotGalaxyShape.c (revision 36756) +++ /trunk/psphot/src/psphotGalaxyShape.c (revision 36757) @@ -89,4 +89,10 @@ opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); psAssert (status, "missing GALAXY_SHAPES_NSIGMA"); + psString modelTypeToSave = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + if (modelTypeToSave && strcmp(modelTypeToSave, "ALL") && strcmp(modelTypeToSave, "BEST")) { + opt->extModelType = pmModelClassGetType(modelTypeToSave); + } else { + opt->extModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + } #ifdef notdef @@ -251,4 +257,5 @@ pmModelType modelType = model->type; + model->flags = PM_MODEL_STATUS_NONE; // we are using fitOptions->mode : be sure this makes sense @@ -322,5 +329,5 @@ // Replace modelEXT with the best model from the first of the model fits if one of them is good - if (isfinite(PAR[PM_PAR_I0]) && iModel == 0) { + if (isfinite(PAR[PM_PAR_I0]) && modelType == opt->extModelType) { psFree (source->modelEXT);