pyxsim.PhotonList.from_data_source() needs to be run twice in order to get results

[SamSweere]

Hi @jzuhone,

I think I found a bug with the PhotonList.from_data_source(). After hours of debugging and getting inconsistent results I found that I need to run twice in order to get any generated photons.

I am currently in the progress of converting Illustrius TNG data to x-ray simput files. The code to re-create the bug is shown below. In order to run it you need the Illustrius TNG cutout which you can download from: https://drive.google.com/file/d/17Easa8wQ8PTG2Sb_tOqeGxH8cG-SLgz5/view?usp=sharing

import yt
import pyxsim

cutout = 'cutout_0.hdf5'
ds = yt.load(cutout)
sp = ds.sphere((43718.8, 48813.6, 147595.), (1000.0, "code_length"))

thermal_model = pyxsim.ThermalSourceModel(
    "apec",
    emin=0.5,
    emax=7.0,
    nchan=10000,
    temperature_field=("PartType0", "Temperature"),
    kT_min=0.1,
    emission_measure_field=("PartType0", "ElectronAbundance"))

redshift = 0.05  # The redshift to the object.
area = (3000., "cm**2")  # A constant collecting effective area to generate the photons with.
exp_time = (100., "ks")  # The exposure time to generate the photons with.
center = sp.center  # A center in 3D for the photon positions. If not specified,
# the center of the `data_source` will be chosen.

# TODO: I believe here is a bug, I need to run it twice to get it to work
# TODO: if you remove the for loop it does not work
for run_num in range(2):
    photons = pyxsim.PhotonList.from_data_source(sp, redshift, area, exp_time,
                                                 thermal_model, center=center)

nH = 0.04  # [10^22 cm-2]
sky_center = (0, 0)  # sky center in RA,DEC
events = photons.project_photons("y", sky_center, absorb_model="wabs", nH=nH)
fov = (20.008, "arcmin")  # the field of view / width of the image
nx = 2440  # The resolution of the image on a side
simput_file = "pyxsim_demonstration_events"
events.write_fits_image(simput_file + '.fits', fov, nx, overwrite=True, emin=0.5, emax=7.0)

# To display the simput image
import aplpy
fig = aplpy.FITSFigure(simput_file + '.fits')
fig.show_colorscale(vmax=100.0, stretch='sqrt', cmap="afmhot")
fig.save("pyxsim_demonstration.png")

Thanks in advance! Sam

[jzuhone]

@SamSweere what versions of yt and pyxsim are you using?

[SamSweere]

yt: Version: 3.6.1 pyxsim: Version: 2.3.0

[jzuhone]

@SamSweere First off, this is not correct:

    emission_measure_field=("PartType0", "ElectronAbundance")

This argument needs the actual emission measure field and not the electron abundance field. If you have the ("PartType0", "ElectronAbundance") field then you should not need to set the emission_measure_field at all because it will be auto-detected.

[SamSweere]

That indeed fixes the problem. My bad. Thank you a lot for the help!