Index: branches/eam_branches/ipp-20110404/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphot.h (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphot.h (revision 31362) @@ -313,5 +313,5 @@ bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf); -bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf); +bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources); bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psImageMaskType maskVal, psImageMaskType markVal); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotApResid.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotApResid.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotApResid.c (revision 31362) @@ -206,5 +206,5 @@ // XXX make this user-configurable? - if (source->errMag > 0.01) continue; + if (source->psfMagErr > 0.01) continue; // aperture residual for this source @@ -222,5 +222,5 @@ source->peak->xf, source->peak->yf, source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], - source->psfMag, source->apMag, source->errMag, + source->psfMag, source->apMag, source->psfMagErr, source->modelPSF->params->data.F32[PM_PAR_I0], source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], @@ -228,5 +228,5 @@ # endif if (!isfinite(source->psfMag)) psAbort ("nan in psfMag"); - if (!isfinite(source->errMag)) psAbort ("nan in errMag"); + if (!isfinite(source->psfMagErr)) psAbort ("nan in psfMagErr"); if (!isfinite(source->apMag)) psAbort ("nan in apMag"); if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos"); @@ -234,5 +234,5 @@ psVectorAppend (mag, source->psfMag); - psVectorAppend (dMag,source->errMag); + psVectorAppend (dMag,source->psfMagErr); psVectorAppend (apResid, dap); psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotChoosePSF.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotChoosePSF.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotChoosePSF.c (revision 31362) @@ -350,5 +350,5 @@ // build curve-of-growth vector for this psf - if (!psphotMakeGrowthCurve (readout, recipe, try->psf)) { + if (!psphotMakeGrowthCurve (readout, recipe, try->psf, try->sources)) { psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF"); psFree (models); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotEfficiency.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotEfficiency.c (revision 31362) @@ -468,8 +468,8 @@ source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], - magRef, source->psfMag, source->errMag); + magRef, source->psfMag, source->psfMagErr); #endif magDiff->data.F32[j] = source->psfMag - magRef; - magErr->data.F32[j] = source->errMag; + magErr->data.F32[j] = source->psfMagErr; } magDiff->n = sources->n; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotExtendedSourceFits.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotExtendedSourceFits.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotExtendedSourceFits.c (revision 31362) @@ -333,5 +333,5 @@ // this uses the footprint to judge both radius and aperture? // XXX save the psf-based moments for output - if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) { + if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) { fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf); // subtract the best fit from the object, leave local sky Index: branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c (revision 31362) @@ -131,5 +131,5 @@ } } - psphotSaveImage (NULL, kronMask, "kronmask.fits"); + // psphotSaveImage (NULL, kronMask, "kronmask.fits"); psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c (revision 31362) @@ -67,5 +67,5 @@ source->psfMag = 0.0; - source->errMag = 0.0; + source->psfMagErr = 0.0; source->apMag = 0.0; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotMakeGrowthCurve.c (revision 31362) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf) { +bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) { bool status; @@ -22,10 +22,18 @@ // measure the aperture loss as a function of radius for PSF - bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS"); float PSF_FIT_PAD = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING"); + + float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); + float PSF_APERTURE = (int)(apScale*gaussSigma); psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS); +# if (0) + bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) { psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections"); @@ -33,4 +41,15 @@ return false; } +# else + bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); + if (!pmGrowthCurveGenerateFromSources (readout, psf, sources, INTERPOLATE_AP, maskVal, markVal)) { + psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections"); + psFree(psf->growth); psf->growth = NULL; + return false; + } +# endif + + float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE); + fprintf (stderr, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset); psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth")); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotOutput.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotOutput.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotOutput.c (revision 31362) @@ -75,5 +75,5 @@ source->peak->xf, source->peak->yf, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], - source->psfMag, source->apMag, source->errMag, + source->psfMag, source->apMag, source->psfMagErr, model->params->data.F32[PM_PAR_I0], model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY], Index: branches/eam_branches/ipp-20110404/psphot/src/psphotSetThreads.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotSetThreads.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotSetThreads.c (revision 31362) @@ -25,5 +25,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10); + task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 11); task->function = &psphotSourceStats_Threaded; psThreadTaskAdd(task); Index: branches/eam_branches/ipp-20110404/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotSourceFits.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotSourceFits.c (revision 31362) @@ -235,5 +235,5 @@ // XXX save the psf-based moments for output psfMoments = *source->moments; - if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) { + if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) { *source->moments = psfMoments; return false; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotSourceSize.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotSourceSize.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotSourceSize.c (revision 31362) @@ -213,5 +213,5 @@ psVectorAppend (ApOff, dMag); - psVectorAppend (ApErr, source->errMag); + psVectorAppend (ApErr, source->psfMagErr); } if (num == 0) { @@ -520,5 +520,5 @@ // set nSigmaMAG to include both systematic and poisson error terms. we include a hard // floor on the Ap Sys Err (to be a bit generous). XXX put the floor in the recipe... - float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02)); + float nSigmaMAG = (dMag - options->ApResid) / hypot(source->psfMagErr, hypot(options->ApSysErr, 0.02)); source->extNsigma = nSigmaMAG; Index: branches/eam_branches/ipp-20110404/psphot/src/psphotSourceStats.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotSourceStats.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotSourceStats.c (revision 31362) @@ -165,4 +165,8 @@ RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); } + float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON"); + if (!status) { + MIN_KRON_RADIUS = 0.75*SIGMA; + } // threaded measurement of the source magnitudes @@ -192,4 +196,5 @@ PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, SIGMA, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); @@ -221,9 +226,9 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[7]; + scalar = job->args->data[8]; Nmoments += scalar->data.S32; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Nfail += scalar->data.S32; - scalar = job->args->data[9]; + scalar = job->args->data[10]; Nfaint += scalar->data.S32; } @@ -382,7 +387,8 @@ float RADIUS = PS_SCALAR_VALUE(job->args->data[3],F32); float SIGMA = PS_SCALAR_VALUE(job->args->data[4],F32); - - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); // if no valid pixels, massive swing or very large Mrf, likely saturated source, @@ -440,5 +446,5 @@ if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) { // measure basic source moments (no S/N clipping on input pixels) - status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); + status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, MIN_KRON_RADIUS, maskVal); } else { // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but @@ -456,5 +462,5 @@ psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal); + status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, MIN_KRON_RADIUS, maskVal); source->mode |= PM_SOURCE_MODE_BIG_RADIUS; } @@ -507,4 +513,5 @@ float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0}; float Sout[NSIGMA]; + float Rmin[NSIGMA]; int Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass @@ -528,5 +535,5 @@ // measure basic source moments (no S/N clipping on input pixels) - status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, maskVal); + status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal); } @@ -539,4 +546,6 @@ pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); + + Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM); #if 0 @@ -564,4 +573,5 @@ // we are looking for sigma for which Sout = 0.65 (or some other value) int Nstars = 0; + float minKronRadius = NAN; float Sigma = NAN; float minS = Sout[0]; @@ -580,4 +590,5 @@ Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); Nstars = 0.5*(Nout[i] + Nout[i+1]); + minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]); } psAssert (isfinite(Sigma), "did we miss a case?"); @@ -590,6 +601,7 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars); - - psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma); + psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_MIN_KRON", PS_META_REPLACE, "minimum Kron Radius", minKronRadius); + + psLogMsg ("psphot", 3, "using window function with sigma = %f (Min Kron Radius = %f, Nstars = %d)\n", Sigma, minKronRadius, Nstars); return true; } Index: branches/eam_branches/ipp-20110404/psphot/src/psphotVisual.c =================================================================== --- branches/eam_branches/ipp-20110404/psphot/src/psphotVisual.c (revision 31337) +++ branches/eam_branches/ipp-20110404/psphot/src/psphotVisual.c (revision 31362) @@ -1867,5 +1867,5 @@ continue; } - if (source->errMag > 0.1) { + if (source->psfMagErr > 0.1) { xLOW->data.F32[nLOW] = source->moments->Mxx; yLOW->data.F32[nLOW] = source->moments->Myy; @@ -2591,5 +2591,5 @@ x->data.F32[n] = source->psfMag; y->data.F32[n] = dMag; - dy->data.F32[n] = source->errMag; + dy->data.F32[n] = source->psfMagErr; graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);