Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 25256) +++ trunk/psphot/src/Makefile.am (revision 25383) @@ -128,5 +128,6 @@ psphotCheckStarDistribution.c \ psphotThreadTools.c \ - psphotAddNoise.c + psphotAddNoise.c \ + psphotEfficiency.c # dropped? psphotGrowthCurve.c Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 25256) +++ trunk/psphot/src/psphot.h (revision 25383) @@ -49,9 +49,9 @@ bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final); bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe); -bool psphotRemoveAllSources (psArray *sources, psMetadata *recipe); +bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe); bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf); @@ -64,5 +64,5 @@ bool psphotGuessModel_Threaded (psThreadJob *job); -bool psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf); +bool psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf); bool psphotMagnitudes_Threaded (psThreadJob *job); @@ -76,4 +76,5 @@ bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); // thread-related: @@ -147,7 +148,7 @@ pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe); bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources); -bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add); +bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe); +bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe); +bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add); bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source); bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 25256) +++ trunk/psphot/src/psphotAddNoise.c (revision 25383) @@ -1,13 +1,13 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { return psphotAddOrSubNoise (readout, sources, recipe, true); } -bool psphotSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { return psphotAddOrSubNoise (readout, sources, recipe, false); } -bool psphotAddOrSubNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) { +bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) { bool status = false; @@ -42,5 +42,5 @@ PS_ASSERT (status, false); if (isfinite(GAIN)) { - FACTOR /= GAIN; + FACTOR /= GAIN; } @@ -70,7 +70,7 @@ oldshape.sxy = PAR[PM_PAR_SXY]; - // XXX can this be done more intelligently? - if (oldI0 == 0.0) continue; - if (!isfinite(oldI0)) continue; + // XXX can this be done more intelligently? + if (oldI0 == 0.0) continue; + if (!isfinite(oldI0)) continue; // increase size and height of source Index: trunk/psphot/src/psphotEfficiency.c =================================================================== --- trunk/psphot/src/psphotEfficiency.c (revision 25383) +++ trunk/psphot/src/psphotEfficiency.c (revision 25383) @@ -0,0 +1,476 @@ +# include "psphotInternal.h" + +#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ + PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models + +//#define TESTING + + +// Calculate the limiting magnitude for an image +// +// We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo +// non-Gaussian PSF) to the limiting flux in the un-smoothed original image. +static bool effLimit(float *magLim, // Limiting magntiude, to return + int *radius, // Radius for fake sources, to return + float *minFlux, // Minimum flux for fake sources, to return + float *norm, // Normalisation of PSF (conversion: peak --> integrated flux) + float *covarFactor,// Covariance factor + const pmReadout *ro, // Readout of interest + const pmPSF *psf, // Point-spread function + float thresh, // Threshold for source identification + float smoothSigma, // Gaussian smoothing sigma + float smoothNsigma, // Smoothing limit + psImageMaskType maskVal // Value to mask + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_VARIANCE(ro, false); + assert(magLim); + assert(radius); + assert(minFlux); + assert(norm); + assert(covarFactor); + + psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing + psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix + psFree(kernel); + *covarFactor = psImageCovarianceFactor(newCovar); // Variance factor + psFree(newCovar); + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, ro->variance, ro->mask, maskVal, rng) || + !isfinite(stats->robustMedian)) { + psError(PSPHOT_ERR_DATA, false, "Unable to determine mean variance"); + psFree(stats); + psFree(rng); + return false; + } + psFree(rng); + float meanVar = stats->robustMedian; // Mean variance + psFree(stats); + + // Need to normalise out difference between Gaussian and real PSF + pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols, + 0.5 * ro->variance->numRows, 1.0); // Model for normalisation + if (!normModel || (normModel->flags & MODEL_MASK)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model."); + psFree(normModel); + return false; + } + *norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 + psFree(normModel); + + // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor) + // since the variance factor tells us the variance in the smoothed image. Now, the trick is working + // out what the intensity in the smoothed image is, and how it is related to the flux. We are + // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star, + // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the + // unsmoothed image, and w is the width of the Gaussian. Now, a normalised 2D Gaussian convolved + // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2. Therefore: + // I_smooth = Flux / 2*norm. + float peakLim = thresh * sqrtf(meanVar * *covarFactor); // Limiting peak value in smoothed image + float fluxLim = 2.0 * *norm * peakLim; // Limiting flux in original + *magLim = -2.5 * log10f(fluxLim); + psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim); + psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim); + psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim); + + *radius = smoothSigma * smoothNsigma; + + *minFlux = 0.1 * sqrtf(meanVar); + + return true; +} + +/// Generate a fake image and add it in to the existing readout +static bool effGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources + const pmReadout *ro, // Readout of interest + const pmPSF *psf, // Point-spread function + const psVector *magOffsets, // Magnitude offsets for fake sources + int numSources, // Number of fake sources for each bin + float refMag, // Reference magnitude + int radius, // Radius for fake sources + float minFlux // Minimum flux for fake sources + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PS_ASSERT_VECTOR_NON_NULL(magOffsets, false); + PS_ASSERT_VECTOR_TYPE(magOffsets, PS_TYPE_F32, false); + assert(xSrc); + assert(ySrc); + + int numBins = magOffsets->n; // Number of bins + int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image + + *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32); + *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32); + + psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + psVector *magAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32); + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + for (int i = 0, index = 0; i < numBins; i++) { + float mag = refMag + magOffsets->data.F32[i]; // Instrumental magnitude of sources + + for (int j = 0; j < numSources; j++, index++) { + xAll->data.F32[index] = psRandomUniform(rng) * numCols; + yAll->data.F32[index] = psRandomUniform(rng) * numRows; + magAll->data.F32[index] = mag; + } + memcpy((*xSrc)->data.F32[i], &xAll->data.F32[i * numSources], + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy((*ySrc)->data.F32[i], &yAll->data.F32[i * numSources], + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psFree(rng); + + pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout + if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll, + NULL, NULL, psf, minFlux, radius, false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image"); + psFree(fakeRO); + psFree(magAll); + return false; + } + psFree(magAll); + + psBinaryOp(ro->image, ro->image, "+", fakeRO->image); + psFree(fakeRO); + + return true; +} + + +// Determine detection efficiency +bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, + psMetadata *recipe, const psArray *realSources) +{ + PM_ASSERT_READOUT_NON_NULL(readout, false); + PM_ASSERT_READOUT_IMAGE(readout, false); + PS_ASSERT_PTR_NON_NULL(psf, false); + PS_ASSERT_METADATA_NON_NULL(recipe, false); + PS_ASSERT_ARRAY_NON_NULL(realSources, false); + + psTimerStart("psphot.fake"); + + // Collect recipe information + float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x + float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y + if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe"); + return false; + } + float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit + if (!isfinite(smoothNsigma)) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe"); + return false; + } + float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2"); + if (!isfinite(thresh)) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe"); + return false; + } + psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets + if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe"); + return NULL; + } + int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin + if (numSources == 0) { + psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe"); + return NULL; + } + float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel + if (!isfinite(minGauss)) { + psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value"); + minGauss = 0.5; + } + + float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image + int numBins = magOffsets->n; // Number of bins + + psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask + + // remove all sources, adding noise for subtracted sources + psphotRemoveAllSources(realSources, recipe); + // psphotAddNoise(readout, realSources, recipe); + + +#if defined(TESTING) && 0 + { + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + readout->image->data.F32[y][x] = psRandomGaussian(rng); + readout->variance->data.F32[y][x] = 1.0; + readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + } + } + psFree(readout->covariance); + readout->covariance = NULL; + psFree(rng); + } +#endif + + + + float magLim; // Guess at limiting magnitude + int radius; // Radius for fake sources + float minFlux; // Minimum flux for fake sources + float norm; // Normalisation of PSF + float covarFactor; // Covariance factor + if (!effLimit(&magLim, &radius, &minFlux, &norm, &covarFactor, readout, + psf, thresh, smoothSigma, smoothNsigma, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image"); + return false; + } + +#ifdef TESTING + psphotSaveImage(NULL, readout->image, "orig_image.fits"); + psphotSaveImage(NULL, readout->variance, "orig_variance.fits"); + psphotSaveImage(NULL, readout->mask, "orig_mask.fits"); +#endif + + psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row + if (!effGenerate(&xFake, &yFake, readout, psf, magOffsets, + numSources, magLim, radius, minFlux)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources"); + psFree(xFake); + psFree(yFake); + return false; + } + +#ifdef TESTING + psphotSaveImage(NULL, readout->image, "fake_image.fits"); + psphotSaveImage(NULL, readout->variance, "fake_variance.fits"); + psphotSaveImage(NULL, readout->mask, "fake_mask.fits"); +#endif + + // XXX Could speed this up significantly by only convolving the central pixels of each fake source + psVector *significance = NULL; // Significance image + { + int num = numSources * numBins; // Total number of sources + // Pixel coordinates of sources + psVector *xPix = psVectorAlloc(num, PS_TYPE_S32); + psVector *yPix = psVectorAlloc(num, PS_TYPE_S32); + for (int i = 0, index = 0; i < numBins; i++) { + for (int j = 0; j < numSources; j++, index++) { + float x = xFake->data.F32[i][j]; + float y = yFake->data.F32[i][j]; + xPix->data.S32[index] = PS_MIN(x + 0.5, numCols - 1); + yPix->data.S32[index] = PS_MIN(y + 0.5, numRows - 1); + } + } + + psVector *convImage = psImageSmoothMaskPixels(readout->image, readout->mask, maskVal, xPix, yPix, + smoothSigma, smoothNsigma, minGauss); // Convolved image + if (!convImage) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth image pixels"); + psFree(xPix); + psFree(yPix); + psFree(xFake); + psFree(yFake); + return false; + } + psVector *convVar = psImageSmoothMaskPixels(readout->variance, readout->mask, maskVal, xPix, yPix, + smoothSigma, smoothNsigma, minGauss); // Convolved varianc + if (!convVar) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to smooth variance pixels"); + psFree(xPix); + psFree(yPix); + psFree(xFake); + psFree(yFake); + return false; + } + + float factor = 1.0 / covarFactor; // Correction for covariance + + const psImage *mask = readout->mask; // Mask for readout + for (int i = 0; i < num; i++) { + float imageVal = convImage->data.F32[i]; // Image value + float varVal = convVar->data.F32[i]; // Variance value + int x = xPix->data.S32[i], y = yPix->data.S32[i]; // Coordinates + if (imageVal < 0 || varVal <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { + convImage->data.F32[i] = 0.0; + } else { + convImage->data.F32[i] = factor * PS_SQR(imageVal) / varVal; + } + } + + psFree(xPix); + psFree(yPix); + + significance = convImage; + psFree(convVar); + } + + thresh *= thresh; // "Significance" is actually significance-squared + + psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin + psArray *fakeSources = psArrayAlloc(numSources); // Fake sources in each bin + psArray *fakeSourcesAll = psArrayAllocEmpty(numSources * numBins); // All fake sources + for (int i = 0, index = 0; i < numBins; i++) { + int numFound = 0; // Number found + + // Determine extraction size + float mag = magLim + magOffsets->data.F32[i]; // Magnitude for bin + float peak = powf(10.0, -0.4 * mag) / norm; // Peak flux + pmModel *model = pmModelFromPSFforXY(psf, numCols / 2.0, numRows / 2.0, peak); // Model for source + if (!model || (model->flags & MODEL_MASK)) { + psError(PS_ERR_UNKNOWN, true, "Unable to generate model for bin %d", i); + psFree(model); + return false; + } + float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source + psFree(model); + + psArray *sources = psArrayAllocEmpty(numSources); // Sources in this bin + for (int j = 0; j < numSources; j++, index++) { + // Coordinates of interest + float sig = significance->data.F32[index]; // Significance of pixel + if (sig > thresh) { + pmSource *source = pmSourceAlloc(); // Fake source + float x = xFake->data.F32[i][j]; + float y = yFake->data.F32[i][j]; + source->peak = pmPeakAlloc(x, y, sig, PM_PEAK_LONE); + if (!pmSourceDefinePixels(source, readout, x, y, sourceRadius)) { + psErrorClear(); + continue; + } + source->peak->xf = x; + source->peak->yf = y; + source->modelPSF = pmModelFromPSFforXY(psf, x, y, 2.0 * sqrtf(sig)); + source->type = PM_SOURCE_TYPE_STAR; + + numFound++; + psArrayAdd(sources, sources->n, source); + psArrayAdd(fakeSourcesAll, fakeSourcesAll->n, source); + psFree(source); + } + } + fakeSources->data[i] = sources; + count->data.S32[i] = numFound; + } + psFree(xFake); + psFree(yFake); + psFree(significance); + + if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); + psFree(fakeSources); + psFree(count); + return false; + } + + // Disable aperture corrections; casting away const! + pmTrend2D *apTrend = psf->ApTrend; // Aperture trend + ((pmPSF*)psf)->ApTrend = NULL; + + if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); + psFree(fakeSources); + psFree(count); + ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const! + return false; + } + psFree(fakeSourcesAll); + + // Re-enable aperture corrections; casting away const! + ((pmPSF*)psf)->ApTrend = apTrend; + + psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin + psVector *magDiffStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of diff in magnitude for each bin + psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin + psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences + psVector *magErr = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude errors + psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors + psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | + PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + for (int i = 0; i < numBins; i++) { + psStatsInit(magStats); + psVectorInit(magMask, 0); + +#ifdef TESTING + psString name = NULL; + psStringAppend(&name, "fake_%d.dat", i); + FILE *file = fopen(name, "w"); + psFree(name); +#endif + + float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude + psArray *sources = fakeSources->data[i]; // Sources in bin + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source || !isfinite(source->psfMag)) { + magMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; + continue; + } + +#ifdef TESTING + fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], + source->modelPSF->params->data.F32[PM_PAR_YPOS], + magRef, source->psfMag, source->errMag); +#endif + magDiff->data.F32[j] = source->psfMag - magRef; + magErr->data.F32[j] = source->errMag; + } + magDiff->n = sources->n; + magMask->n = sources->n; + magErr->n = sources->n; + if (!psVectorStats(magStats, magDiff, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + + if (isfinite(magStats->robustMedian) && isfinite(magStats->robustStdev)) { + magDiffMean->data.F32[i] = magStats->robustMedian; + magDiffStdev->data.F32[i] = magStats->robustStdev; + } else { + magDiffMean->data.F32[i] = magStats->sampleMean; + magDiffStdev->data.F32[i] = magStats->sampleStdev; + } + + if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) { + // Probably because we don't have enough sources + psErrorClear(); + } + + if (isfinite(magStats->robustMedian)) { + magErrMean->data.F32[i] = magStats->robustMedian; + } else { + magErrMean->data.F32[i] = magStats->sampleMean; + } + +#ifdef TESTING + fclose(file); +#endif + } + + psFree(magStats); + psFree(magDiff); + psFree(magMask); + psFree(magErr); + + psFree(fakeSources); + + pmDetEff *de = pmDetEffAlloc(magLim, numSources, numBins); // Detection efficiency + de->magOffsets = psVectorCopy(NULL, magOffsets, PS_TYPE_F32); + de->counts = count; + de->magDiffMean = magDiffMean; + de->magDiffStdev = magDiffStdev; + de->magErrMean = magErrMean; + + psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PM_DETEFF_ANALYSIS, PS_META_REPLACE | PS_DATA_UNKNOWN, + "Detection efficiency", de); + psFree(de); + + psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake")); + + + return true; +} Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 25256) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 25383) @@ -12,5 +12,5 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) { bool status; Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 25256) +++ trunk/psphot/src/psphotMagnitudes.c (revision 25383) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) { +bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) { bool status = false; @@ -64,5 +64,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, (pmPSF*)psf); // Casting away const psArrayAdd(job->args, 1, binning); psArrayAdd(job->args, 1, backModel); Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 25256) +++ trunk/psphot/src/psphotReadout.c (revision 25383) @@ -61,5 +61,5 @@ // display the backsub and backgnd images psphotVisualShowBackground (config, view, readout); - + // run a single-model test if desired (exits from here if test is run) psphotModelTest (config, view, recipe); @@ -224,4 +224,9 @@ psphotMagnitudes(config, readout, view, sources, psf); + if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); + psErrorClear(); + } + // replace failed sources? // psphotReplaceUnfitSources (sources); Index: trunk/psphot/src/psphotReadoutMinimal.c =================================================================== --- trunk/psphot/src/psphotReadoutMinimal.c (revision 25256) +++ trunk/psphot/src/psphotReadoutMinimal.c (revision 25383) @@ -97,4 +97,9 @@ psphotMagnitudes(config, readout, view, sources, psf); + if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); + psErrorClear(); + } + // drop the references to the image pixels held by each source psphotSourceFreePixels (sources); Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 25256) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 25383) @@ -47,5 +47,5 @@ } -bool psphotRemoveAllSources (psArray *sources, psMetadata *recipe) { +bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) { bool status; Index: trunk/psphot/src/psphotSignificanceImage.c =================================================================== --- trunk/psphot/src/psphotSignificanceImage.c (revision 25256) +++ trunk/psphot/src/psphotSignificanceImage.c (revision 25383) @@ -22,20 +22,23 @@ } - bool status_x, status_y; - float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); - float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); - if (status_x && status_y) { - // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?) - SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); - guess = false; + bool statusMajor, statusMinor; + float fwhmMajor = psMetadataLookupF32(&statusMajor, recipe, "FWHM_MAJ"); + float fwhmMinor = psMetadataLookupF32(&statusMinor, recipe, "FWHM_MIN"); + if (statusMajor && statusMinor) { + // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?) + if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) { + psWarning("fwhmMajor (%f) or fwhmMinor (%f) is bad!", fwhmMajor, fwhmMinor); + } + SIGMA_SMTH = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0))); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + guess = false; } else { - // if we do not know the FWHM, use the guess smoothing kernel supplied. - // it is a configuration error if these are not supplied - SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); - PS_ASSERT (status, NULL); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); - PS_ASSERT (status, NULL); - guess = true; + // if we do not know the FWHM, use the guess smoothing kernel supplied. + // it is a configuration error if these are not supplied + SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); + PS_ASSERT (status, NULL); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + PS_ASSERT (status, NULL); + guess = true; } // record the actual smoothing sigma