Error in the flux when summing achromatic opticalPSF using photon shooting

[aguinot]

Hello,
I found that when using photon shooting with a Sum of opticalPSF the output flux is biased. The ChromaticSum works fine.
Bellow is a short example to reproduce the effect.

import galsim
import galsim.roman as roman

bp = roman.getBandpasses()["J129"]
vega_sed = galsim.SED('vega.txt', 'nm', 'flambda')

seed = 42
flux = 1000.

# Achromatic PSF
psf_1 = galsim.OpticalPSF(
    diam=roman.diameter,
    lam=bp.effective_wavelength,
)
psf_tot = galsim.Add([psf_1, psf_1])/2
# ChromaticPSF
psf_1_chrom = galsim.ChromaticOpticalPSF(
    diam=roman.diameter,
    lam=bp.effective_wavelength,
)
psf_tot_chrom = galsim.Add([psf_1_chrom, psf_1_chrom])/2

# Make example object
star = galsim.DeltaFunction(flux=flux)
star_chrom = (star * vega_sed).withFlux(flux, bp)

# Chromatic drawing
gal_rng = galsim.BaseDeviate(seed)
img_chrom = star_chrom.drawImage(
    bandpass=bp,
    nx=101,
    ny=101,
    method="phot",
    n_photons=flux,
    photon_ops=[psf_tot_chrom],
    rng=gal_rng,
    scale=0.1,
)

# Achromatic drawing
gal_rng = galsim.BaseDeviate(seed)
img_achrom = star.drawImage(
    nx=101,
    ny=101,
    method="phot",
    n_photons=flux,
    photon_ops=[psf_tot],
    rng=gal_rng,
    scale=0.1,
)

# Check fluxes
print("Chromatic PSF flux:", np.sum(img_chrom.array))
print("Achromatic PSF flux:", np.sum(img_achrom.array))
print("True flux:", flux)

Chromatic PSF flux: 1063.5648
Achromatic PSF flux: 1577.9474
True flux: 1000.0