A bug in soxs.simulate_spectrum

[liuguanfu1120]

Hi, I think I found a bug in soxs.simulate_spectrum

Bug description

InLine 934 of the file instrument.py of version 4.8.3 (the latest version), it writes

bkgnd_area = np.sqrt(parse_value(bkgnd_area, "arcmin**2"))

Such a line converts the background area (bkgnd_area) in arcmin**2 to an angle in arcmin.

However, in lines 956-958,

        out_spec += generate_channel_spectrum(
            frgnd_spec, exp_time, bkgnd_area, prng=prng
        )

the angle (in arcmin!) is passed to the function generate_channel_spectrum.
The first several lines of generate_channel_spectrum

def generate_channel_spectrum(count_rate, t_exp, solid_angle, prng=None):
    """
    Generate photon energy channels from this diffuse
    background spectrum given an exposure time,
    effective area, and solid angle.

    Parameters
    ----------
    t_exp : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The exposure time in seconds.
    solid_angle : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The solid angle in arcmin**2.
    prng : :class:`~numpy.random.RandomState` object, integer, or None
        A pseudo-random number generator. Typically will only
        be specified if you have a reason to generate the same
        set of random numbers, such as for a test. Default is None,
        which sets the seed based on the system time.
    """

I think we cannot pass the angle to generate_channel_spectrum and we need to pass a background area in arcmin**2.

How to fix it

I think we could just delete the line 934.

What is the problem?

I would like to post a demo here to show the problem if the line 934 still exists, which is actually how I found such a bug.

import matplotlib
matplotlib.rc("font", size=18)
import matplotlib.pyplot as plt
import soxs
import os
os.environ["XDG_CONFIG_HOME"] = "/Users/liuguanfu/.config"
from astropy.io import fits
instrument = "lem_0.9eV"
exp_time = (218.7, "ks")  # in seconds
fov = soxs.instrument_registry[instrument]['fov']  # in arcmin
sky_center = [22.0, -27.0]  # in degrees, [RA, DEC]
output_file="Test.simput"
name="test"
soxs.make_point_sources_file(
    output_file,
    name,
    exp_time,
    fov,
    sky_center,
    prng=24,  # Random seed
    output_sources="test.dat",
    overwrite=True,
    area=(4000.0, "cm**2"),
)
input_events=output_file
foreground=True
out_file="Test-Galactic-events.fits"
soxs.instrument_simulator(
    input_events,
    out_file,
    exp_time,
    instrument,
    sky_center,
    overwrite=True,
    instr_bkgnd=False,
    foreground=foreground,
    ptsrc_bkgnd=False,  # The point source is input, this shoule be False
    prng=24
)
CXB_events=out_file
soxs.write_spectrum(CXB_events, 
                    "Test-Galactic-events.pi", 
                    overwrite=True)
bkg1 = fits.open("Test-Galactic-events.pi")
x, y = bkg1['SPECTRUM'].data['CHANNEL'], bkg1['SPECTRUM'].data['COUNTS']
fig, ax = plt.subplots()
ax.plot(x, y)
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_title('Background Spectrum from instrument_simulator')
plt.show()
out_file = "Test.pha"
soxs.simulate_spectrum(None, instrument, exp_time, out_file,
                  ptsrc_bkgnd=True, foreground=foreground,
                  instr_bkgnd=False, overwrite=True,
                  prng=24,  # Random seed
                  bkgnd_area=(32*32, "arcmin**2")) # 1024=32*32
bkg2 = fits.open(out_file)
x, y = bkg2['SPECTRUM'].data['CHANNEL'], bkg2['SPECTRUM'].data['COUNTS']
fig, ax = plt.subplots()
ax.plot(x, y)
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_title('Background Spectrum from simulate_spectrum')
plt.show()

You will see the counts of the spectrum from instrument_simulator is about 30 times of simulate_spectrum.
The spectrum is the sum of cosmic X-ray background (CXB) and Galactic foreground and the Galactic foreground dominates the CXB.
Although simulate_spectrum uses a simplified model to simulate the CXB, it should not deviate so much from that of instrument_simulator.

The resultant figures are

[liuguanfu1120]

How to fix it

Line 952 of the file instrument.py of version 4.8.3 (the latest version) should be modified as follows if line 934 is deleted:

fov = None if bkgnd_area is None else np.sqrt(parse_value(bkgnd_area, "arcmin**2"))

[jzuhone]

Fixed in version 4.8.4, thanks for the report!