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