Index: trunk/psModules/src/objects/pmSourceIO_CMF.c.in =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 35768) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c.in (revision 36375) @@ -789,4 +789,6 @@ char name[64]; + pmModelType modelTypeTrail = pmModelClassGetType("PS_MODEL_TRAIL"); + // create a header to hold the output data psMetadata *outhead = psMetadataAlloc (); @@ -798,9 +800,9 @@ sources = psArraySort (sources, pmSourceSortByFlux); - @>PS1_DV2@ float magOffset; - @>PS1_DV2@ float zeroptErr; - @>PS1_DV2@ float fwhmMajor; - @>PS1_DV2@ float fwhmMinor; - @>PS1_DV2@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); + float magOffset; + float zeroptErr; + float fwhmMajor; + float fwhmMinor; + pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader); // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams @@ -823,4 +825,12 @@ @>PS1_DV2@ pmChip *chip = readout->parent->parent; + pmModelStatus badModel = PM_MODEL_STATUS_NONE; + badModel |= PM_MODEL_STATUS_BADARGS; + badModel |= PM_MODEL_STATUS_OFFIMAGE; + badModel |= PM_MODEL_STATUS_NAN_CHISQ; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS; + badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID; + badModel |= PM_MODEL_PCM_FAIL_GUESS; + table = psArrayAllocEmpty (sources->n); @@ -847,5 +857,6 @@ // skip models which were not actually fitted - if (model->flags & PM_MODEL_STATUS_BADARGS) continue; + // XXX + if (model->flags & badModel) continue; PAR = model->params->data.F32; @@ -853,6 +864,21 @@ xPos = PAR[PM_PAR_XPOS]; yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; + + // for the extended source models, we do not always fit the centroid in the non-linear fitting process + // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM, + // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM) + // TRAIL : X,Y are fitted + // + + // XXX this should be based on what happened, not on the model type + if (model->type == modelTypeTrail) { + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + } else { + // this is definitely an underestimate since it does not + // account for the extent of the source + xErr = fwhmMajor * model->magErr / 2.35; + yErr = fwhmMinor * model->magErr / 2.35; + } @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0}; @@ -870,5 +896,4 @@ row = psMetadataAlloc (); - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment" @@ -888,6 +913,10 @@ @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN; @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration", calMag); - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Magnitude using supplied calibration", model->chisq); - @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Magnitude using supplied calibration", model->nDOF); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ", PS_DATA_F32, "EXT Model Chisq", model->chisq); + @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF", PS_DATA_S32, "EXT Model num degrees of freedom", model->nDOF); + @>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); + + // EAM : adding for PV2 outputs: + @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0); @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE", 0, "position angle at source (degrees)", posAngle); @@ -946,5 +975,5 @@ snprintf (name, 64, "EXT_PAR_%02d", k); - + if (k < model->params->n) { psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]); @@ -956,5 +985,5 @@ // optionally, write out the covariance matrix values // XXX do I need to pad this to match the biggest covar matrix? - if (model->covar) { + if (false && model->covar) { for (int iy = 0; iy < model->covar->numCols; iy++) { for (int ix = iy; ix < model->covar->numCols; ix++) { @@ -1058,4 +1087,17 @@ model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ"); + model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF"); + + // EXT_MODEL_TYPE gives the model chosen by psphot as the best. + // Putting this into the XFIT table makes 3 copies of it (one for each model) + // but since we have many fewer XFIT rows than psf rows that is cheaper than putting it + // in the psf table. + psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE"); + if (!status) { + // older cmfs don't have this column + extModelType = -1; + } + psEllipseAxes axes; axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); @@ -1072,20 +1114,40 @@ if (model->params->n > 7) { PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); - } - // read the covariance matrix - int nparams = model->params->n; - psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); - for (int y = 0; y < nparams; y++) { - for (int x = 0; x < nparams; x++) { - char name[64]; - snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); - covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + // XXX add an error: + // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_"); + } + + // XXX : make this depend on what is in the cmf + if (0) { + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + } + + if (modelType == extModelType) { + // The software that created this source picked this model as the best of the fits. + // Set the extModel to point to it. + // This is important for programs like psastro (skycal) so that its output cmfs + // will have valid EXT_MODEL_TYPE + psFree(source->modelEXT); + source->modelEXT = psMemIncrRefCounter(model); + source->type = PM_SOURCE_TYPE_EXTENDED; + if (0) { + // since FLAGS were read we don't need to do this + source->mode |= PM_SOURCE_MODE_EXTMODEL; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; } } - model->covar = covar; psArrayAdd(source->modelFits, 1, model); psFree(model); - psFree(row); } @@ -1209,8 +1271,8 @@ write_annuli: - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_DATA_VECTOR, "flux error in annuli", radFluxErr); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation", radFluxStdev); - psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX", PS_META_REPLACE, "flux within annuli", radFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR", PS_META_REPLACE, "flux error in annuli", radFluxErr); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation", radFluxStdev); + psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL", PS_META_REPLACE, "fill factor of annuli", radFill); psFree (radFlux); psFree (radFluxErr); @@ -1343,2 +1405,88 @@ return true; } + +// XXX where should I record the number of columns?? +bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) +{ + bool status = false; + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) { + psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n"); + return true; + } + + // create a header to hold the output data + psMetadata *outhead = psMetadataAlloc (); + + // write the links to the image header + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname); + + // let's write these out in S/N order + sources = psArraySort (sources, pmSourceSortByFlux); + + psArray *table = psArrayAllocEmpty (sources->n); + + for (int i = 0; i < sources->n; i++) { + + pmSource *thisSource = sources->data[i]; + + // this is the "real" version of this source + pmSource *source = thisSource->parent ? thisSource->parent : thisSource; + + // if we did not fit the galaxy model, modelFits will be NULL + if (source->modelFits == NULL) continue; + + // if we did not fit the galaxy model, galaxyFits will also be NULL + if (source->galaxyFits == NULL) continue; + + pmModel *model = source->modelFits->data[0]; + if (!model) return false; + + // X,Y coordinates are stored with the model parameters + psF32 *PAR = model->params->data.F32; + + psMetadata *row = psMetadataAlloc (); + + // we write out the x,y positions so people can link to the psf either way (position or ID) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT", 0, "model x coordinate", PAR[PM_PAR_XPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT", 0, "model y coordinate", PAR[PM_PAR_YPOS]); + psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX", 0, "number of pixels for fits", source->galaxyFits->nPix); + + psVector *Flux = source->galaxyFits->Flux; + psVector *dFlux = source->galaxyFits->dFlux; + psVector *chisq = source->galaxyFits->chisq; + + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX", PS_META_REPLACE, "normalization for galaxy flux", Flux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux); + psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ", PS_META_REPLACE, "galaxy fit chisq", chisq); + + psArrayAdd (table, 100, row); + psFree (row); + } + + if (table->n == 0) { + if (!psFitsWriteBlank (fits, outhead, extname)) { + psError(psErrorCodeLast(), false, "Unable to write empty sources file."); + psFree(outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; + } + + psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname); + if (!psFitsWriteTable (fits, outhead, table, extname)) { + psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); + psFree (table); + return true; +} +