Fix an error in how ChromaticConvolution handled poisson_flux=True

(#1229)

CHANGELOG.rst

@@ -44,3 +44,5 @@ Bug Fixes
 
 - Fixed a bug that prevented Eval types from generating lists in config files in some contexts.
 - Changed the SED class to correctly broadcast over waves when the SED is constant. (#1228)
+- Fixed some errors when drawing ChromaticTransformation objects with photon shooting. (#1229)
+- Fixed the flux drawn by ChromaticConvolution with photon shooting when poisson_flux=True. (#1229)

galsim/chromatic.py

@@ -2674,7 +2674,8 @@ def drawImage(self, bandpass, image=None, integrator='quadratic', iimult=None, *
             poisson_flux = kwargs.pop('poisson_flux', n_photons == 0.)
             max_extra_noise = kwargs.pop('max_extra_noise', 0.)
             rng = BaseDeviate(kwargs.get('rng', None))
-            n_photons, _ = prof1._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
+            n_photons, g = prof1._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
+            gal *= g
             return gal.drawImage(bandpass, image=image, integrator=integrator,
                                  n_photons=n_photons, **kwargs)
 

tests/test_chromatic.py

@@ -2967,6 +2967,12 @@ def test_shoot_transformation():
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux)
 
+    # We used to do the wrong thing with poisson_flux=True.  Check that the flux isn't
+    # just 1000 in that case.
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
+
     # Rotate
     psf = galsim.ChromaticObject(
         galsim.Gaussian(fwhm=1)
@@ -2976,6 +2982,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Expand
     psf = galsim.ChromaticObject(
@@ -2989,16 +2998,22 @@ def test_shoot_transformation():
     # With only 1000 photons, it only matches to better than 2.e-3.
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Shift
     psf = galsim.ChromaticObject(
         galsim.Gaussian(fwhm=1)
-    ).shift(lambda w: ((w/500)*0.3, (w/500)*0.9))
+    ).shift(lambda w: ((w/500)*0.3, (w/500)*-0.4))
     obj = galsim.Convolve(psf, star).withFlux(flux, bandpass)
     img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng,
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Shear
     psf = galsim.ChromaticObject(
@@ -3009,6 +3024,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Magnify
     psf = galsim.ChromaticObject(
@@ -3019,6 +3037,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Lens
     psf = galsim.ChromaticObject(
@@ -3029,6 +3050,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Transform
     psf = galsim.ChromaticObject(
@@ -3040,6 +3064,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux, rtol=2.e-3)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
     # Flux_scale
     psf = galsim.ChromaticObject(
@@ -3050,6 +3077,9 @@ def test_shoot_transformation():
                         poisson_flux=False)
     print(img.added_flux)
     np.testing.assert_allclose(img.added_flux, flux)
+    img = obj.drawImage(bandpass, nx=25, ny=25, scale=0.2, method='phot', rng=rng)
+    print(img.added_flux)
+    assert abs(img.added_flux - flux) > 0.1
 
 
 if __name__ == "__main__":