Poisson noise added in incorrect units in accelerated PSF injection (v0.12.0)

[AstroJacobLi]

Thanks for releasing the new version of romanisim (v0.12.0). I was excited to test the newly introduced accelerated PSF injection feature. While validating it against the previous (slow) implementation, I identified a bug related to how Poisson noise is applied to point sources.

In romanisim.image.add_objects_to_image, the input fluxes for point sources are expected to be in maggies, and the conversion to detector counts is handled via flux_to_counts_factor.

In the original (slow) PSF injection path, the workflow is:

1. Draw the PSF into a stamp in electrons.
2. Add Poisson noise to the stamp
   

   romanisim/romanisim/image.py

   Line 332 in 25fe17e

   stamp.addNoise(galsim.PoissonNoise(rng))

3. Convert the stamp to the output unit (typically MJy/sr)
   

   romanisim/romanisim/image.py

   Line 337 in 25fe17e

   stamp /= outputunit_to_electrons[i]

This procedure is physically correct, since Poisson noise should be applied in count space.

However, in the newly added accelerated code, the flux is converted directly to MJy/sr before noise is added:

romanisim/romanisim/image.py

Lines 286 to 287 in 25fe17e

	if outputunit_to_electrons is not None:
	flux2counts /= np.array(outputunit_to_electrons)

As a result, both the PSF stamp and image_pointsources are already in MJy/sr when Poisson noise is applied

romanisim/romanisim/image.py

Lines 368 to 369 in 25fe17e

	if np.sum(pointsources) > 0 and add_noise:
	image_pointsources.addNoise(galsim.PoissonNoise(rng))

This leads to Poisson noise being added in surface-brightness units rather than electrons, which is inconsistent with the detector noise model and differs from the behavior of the original implementation.

Proposed fix

A minimal fix would be to delay the unit conversion until after the Poisson noise is applied. Concretely, this could be achieved by removing lines 286–287 and adding something like:

if np.sum(pointsources) > 0 and add_noise:
        image_pointsources.addNoise(galsim.PoissonNoise(rng)) # noise should be added when the unit is electrons
if outputunit_to_electrons is not None:
    image_pointsources /= np.array(outputunit_to_electrons).mean()

This assumes that outputunit_to_electrons is a constant for all sources, which should be true in most cases.

Again, thanks for creating this new PSF injection feature! Hope the above makes sense : )

[schlafly]

Thanks for the detailed issue report. I agree with your assessment. We will need to handle the outputunit_to_electrons case to handle source injection into mosaics where different parts of the mosaic have different number of input units, but I guess that can be another special case so that we don't slow down the ordinary case.

I'll take a look at fixing this later this week---good catch!

[schlafly]

Would you mind seeing if #313 fixes your case? Thanks!

[AstroJacobLi]

@schlafly I replied here #313 (comment). Do you mind checking it? Thanks!