ChromaticSum.drawImage sometimes shoots too many photons

[jmeyers314]

If n_photons is explicitly set in ChromaticSum.drawImage when drawing with method='phot', then that many photons will be used for each summand, instead of that many photons being used amongst all summands. This can also happen if ChromaticConvolution includes a ChromaticSum as a convolutant even when n_photons is not explicitly specified. I think we may need some logic similar to galsim.Sum.shoot inside ChromaticSum to get this right.

Also, once photons have been drawn from each summand profile, I think it's always okay to pass them through subsequent photons ops together. I.e., there's no need to repeatedly call the same photon op for each summand, and in fact avoiding this may lead to some small speedup.

Code to reproduce:

import numpy as np
import galsim

sed1 = galsim.SED("CWW_E_ext.sed", wave_type='A', flux_type='flambda')
sed2 = galsim.SED("CWW_Sbc_ext.sed", wave_type='A', flux_type='flambda')
bandpass = galsim.Bandpass("LSST_r.dat", wave_type="nm")

obj1 = (galsim.Gaussian(fwhm=0.6) * sed1).withFlux(990.0, bandpass)
obj2 = (galsim.Gaussian(fwhm=0.3) * sed2).withFlux(10.0, bandpass)
obj = obj1 + obj2

class Counter(galsim.PhotonOp):
    def __init__(self):
        self.nphot = []
        self.meanflux = []

    def applyTo(self, photon_array, local_wcs=None, rng=None):
        self.nphot.append(len(photon_array))
        self.meanflux.append(np.mean(photon_array.flux))

# Looks okay when n_photons unspecified
counter = Counter()
img = obj.drawImage(
    bandpass, nx=24, ny=24, scale=0.2, method='phot',
    photon_ops=[counter]
)
print(f"{counter.nphot = }")
print(f"{counter.meanflux = }")
# RNG-dependent, but I got:
# counter.nphot = [962, 7]
# counter.meanflux = [1.0, 0.9999999999999999]

# When n_photons is explicit, shoots too many photons
counter = Counter()
img = obj.drawImage(
    bandpass, nx=24, ny=24, scale=0.2, method='phot',
    photon_ops=[counter], n_photons=101, poisson_flux=False
)
print(f"{counter.nphot = }")
print(f"{counter.meanflux = }")
# counter.nphot = [101, 101]  # Shot 202 photons instead of 101!
# counter.meanflux = [9.801980198019804, 0.09900990099009901]  weights produce right fluxes, but will have wrong Poisson statistics.

[rmjarvis]

@jmeyers314 I'm pretty sure I know what to do to fix this. But let me know if you're already working on it.

[jmeyers314]

Nope, haven't started. Please go ahead.