Merge pull request #86 from GalSim-developers/fb/psf_normalization

Update Euclid-like PSF normalization to remove obscuration

euclidlike/euclidlike_psf.py

@@ -5,7 +5,7 @@
 import astropy.io.fits as pyfits
 from importlib.resources import files
 from galsim.utilities import LRU_Cache
-from . import n_ccd, n_pix_col, n_pix_row, pixel_scale, diameter, obscuration, det2ccd
+from . import n_ccd, n_pix_col, n_pix_row, pixel_scale, diameter, obscuration, det2ccd, collecting_area
 from .bandpass import getBandpasses
 from galsim.config import LoggerWrapper
 import pathlib
@@ -31,9 +31,14 @@ def _make_psf_list(psf_file):
     scale = pixel_scale/3  # images are oversampled by a factor of 3
     im_list = []
     nsample = len(image_array )
+    # Provided PSF is normalized such that it includes the effects of obscuration; we want to remove it as obscuration is internally handled by Galsim by using the effective collecting area to compute the measured flux.
+    # Diameter is in meters, collecting area in cm^2, need to convert diameter to cm
+    m2cm_conv = 1e2
+    obscuration_norm = collecting_area / ((m2cm_conv*diameter/2)**2*np.pi)
     for i in range(nsample):
+        psf_arr = image_array[i]/obscuration_norm
         im_list.append(
-            galsim.Image(image_array[i], scale=scale)
+            galsim.Image(psf_arr, scale=scale)
         )
     return im_list
 
@@ -82,6 +87,11 @@ def getPSF(
     PSF images are stored using the focal plane position format. Therefore,
     we convert the CCD detector ID to the appropiate focal plane position
     internally.
+
+    In addition, the provided PSF images are normalized for obscuration, 
+    vignetting and baffle effects. However, GalSim internally handles the 
+    obscuration, so we remove this part of the normalization by dividing the
+    PSF images by collecting_area / ((diameter/2)**2*np.pi).  As a result, the sum of the pixel values in the renormalized PSF images is very close to 1.
     
 
     Args:

tests/test_psf.py

@@ -37,6 +37,12 @@ def test_get_psf_function():
     np.testing.assert_allclose(psfobj.image.array, trueobj.array, atol = 0,
         err_msg = 'getPSF() fails to initialize input images correctly')
 
+    # check sum of PSF image pixels after normalizing by obscuration is within 1% of 1.
+    # It won't add to exactly 1 since PSF includes vignetting and baffle effects.
+    for i in range(len(psfobjs)):
+          np.testing.assert_allclose(np.sum(psfobjs[i].array), 1.0, atol = 0.01,
+              err_msg = 'Sum of PSF pixels != 1 by more than 1%') 
+
     # check different ccd works
     psfobj = getPSF(ccd=test_adj_ccd, bandpass="VIS", psf_dir = psf_dir)