How to accurately compute magnitude and magnitude errors for detected sources with IterativePSFPhotometry

[ferromatteo]

Hello,

I'm trying to get flux estimate (and magnitudes) from IterativePSFPhotometry and related errors, but I have some doubts, this is part of my code:

_, _, std = sigma_clipped_stats(astroImage)
fit_shape = (5,5)
finder = DAOStarFinder(2.0 * std, max(densest_x, densest_y))
psf_model = GaussianPSF(flux=1, x_fwhm=densest_x, y_fwhm=densest_y)
sigma_clip = SigmaClip(sigma=3.0)
bkgstat = MMMBackground(sigma_clip=sigma_clip)
localbkg_estimator = LocalBackground(1.5 * densest_ave, 3 * densest_ave, bkgstat)
AR = max(densest_x, densest_y)
grouper = SourceGrouper(min_separation = 1.5 * max(densest_x, densest_y))

psfphot = IterativePSFPhotometry(psf_model, fit_shape, finder=finder,
                                 aperture_radius=AR, grouper=grouper, 
                                 localbkg_estimator=localbkg_estimator)

bkg_estimator = MedianBackground()
bkg = Background2D(astroImage, (50, 50), 
                   filter_size=(3, 3),
                   sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)

data = astroImage - bkg.background
phot = psfphot(data)

Where densest_x, densest_y, densest_ave are just my estimate of the FWHM(x,y,mean) in pixels (from my previous code), astroImage is my fits data matrix (not background subtracted), and the fit_shape is put by hand since an odd number is needed and I prefer to check first, usually I put something a little larger than my estimate of the PSF.

I would like to know if this kind of approach makes sense, since I read here ([https://photutils.readthedocs.io/en/latest/user_guide/psf.html]) that my data should be already background_subtracted prior PSF photometry, and this is what I basically do in the second part with a median Background2D estimation. My question is:

In this case it still makes sense to include a local background estimation in the IterativePSFPhotometry call, as I do, Can be beneficial in cases of crowded part of the image, or it is just handled togheter by the grouping part and consequent simultaneous fit?

I would also like to know if the general guesses for aperture, fit_shape and so on make sense with respect to the measured FWHM and if, once the photometry is done, the flux_err in the table of the results can be used to compute the magnitude error as:

def flux2mag(flux,flux_err):
    mag = - 2.5 * np.log10(flux)
    mag_err = (2.5 * flux_err) / (flux * np.log(10))
    return mag, mag_err

Or it needs some other considerations? Such as being dependent on the error that one can include in phot = psfphot(data), that if I understood correctly could as a first approximation be set to the standard deviation of the global background estimation.

Thank you very much in advance,

Matteo

[larrybradley]

In this case it still makes sense to include a local background estimation in the IterativePSFPhotometry call, as I do, Can be beneficial in cases of crowded part of the image, or it is just handled togheter by the grouping part and consequent simultaneous fit?

In very crowded fields, I would be a bit worried about using a local background estimation because it'll likely be biased by close neighboring sources. The local background is subtracted over the fit_shape region prior to performing the PSF fitting. As noted in the docs (see the Notes section), this can also lead to over-subtraction in regions where the fit_shape of sources overlap:

For sources where their fit_shape regions overlap, the local background will effectively be 
subtracted twice in the overlapping fit_shape regions, even if the source grouper is input. 
This is not an issue if the sources are well-separated. However, for crowded fields, please 
use the localbkg_estimator (or local_bkg column in init_params) with care.

I would also like to know if the general guesses for aperture, fit_shape and so on make sense with respect to the measured FWHM and if, once the photometry is done,

That sounds like a reasonable for the initial aperture size. The input is aperture_radius (not diameter) so you'd probably want to use FWHM/2. For the fit_shape, it is recommended to use a small region where the signal to noise is largest. For space-based images, this is typically 5x5 pixels.

if, once the photometry is done, the flux_err in the table of the results can be used to compute the magnitude error as:

Since magnitude uses a log transformation, symmetric flux errors become asymmetric magnitude errors. However, if the SNR is large, they can be approximated as a symmetric error bar with mag_err = 2.5 * np.log10(1.0 + flux_err / flux).

The input error array should represent all sources of error. This error array is fully propagated in the output errors if you are using astropy 5.3+.