Spectra with CountSpectrum from event files are different from loading the same data with Xspec/astropy using pha files

[matteolucchini1]

Description of the Bug

Per the title, the spectra loaded with the two different methods look slightly different. The script below uses the same obsid as the documentation, and compares the (detector space) spectra produced by 1) Stingray 2) pyXspec 3) loading in the ungrouped data with Astropy and 4) loading in the grouped data with Astropy. For the last two, I loaded the response matrix using nDspec, but that should not affect the results (that part of the code just wraps around Astropy).

The obsid is identical to the obsid used in the documentation, but I re-reduced the observation on my machine and used my own event file. The spectrum/response files were produces taking my own event file and running it through nicerl3 as standard.

Steps/Code to Replicate the Bug

import os
import sys
import warnings
import copy

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.pylab as pl
from matplotlib import rc, rcParams
rc('text',usetex=True)
rc('font',**{'family':'serif','serif':['Computer Modern']})
fi = 22
plt.rcParams.update({'font.size': fi-5})

colorscale = pl.cm.PuRd(np.linspace(0.,1.,5))

from stingray import EventList, CountSpectrum

sys.path.append('/home/matteo/Software/nDspec/src/')
from ndspec.Response import ResponseMatrix

#need to load a response to convert to energy or I can't compare to xspec
obsid = str(106)
path = os.getcwd() + "/../EventFiles/"

rmfpath = path+"1200120106_rmf.pha"
nicer_matrix = ResponseMatrix(rmfpath)
arfpath = path+"1200120106_arf.pha"
nicer_matrix.load_arf(arfpath)

#load the event file and extract the spectrum in stingray
fname = path + "ni1200120" + obsid+ "_0mpu7_cl.evt"
events = EventList.read(fname, "hea", additional_columns=["DET_ID"],rmf_file=rmfpath)
energy_spec = np.geomspace(0.5, 10, 51)
countsp = CountSpectrum(events, energy_spec=energy_spec)
bin_width = np.diff(energy_spec)

#load the spectrum in xspec - neglect the background, although for this obsid it's irrelevant 
from xspec import *
AllData.clear()
s = Spectrum(path+"1200120106_rebin.pha")
s.response = path+"1200120106_rmf.pha"
s.response.arf = path+"1200120106_arf.pha"
s.ignore("**-0.5 10.0-**")
Plot.xAxis = "keV"
Plot("data")
xVals = Plot.x()
yVals = Plot.y()
yErrs = Plot.yErr()

#load the data directly from the .pha with astropy, grabbing both grouped and ungrouped spectra 
#now to load the same with astropy:
from astropy.io import fits
spectrum_path = path + "1200120106_rebin.pha"
with fits.open(spectrum_path,filemap=False) as spectrum:
    extnames = np.array([h.name for h in spectrum])
    exposure = spectrum['PRIMARY'].header['EXPOSURE']
    spectrum_data = spectrum['SPECTRUM'].data
    channels = spectrum_data['CHANNEL']
    counts = spectrum_data['COUNTS']
    sys_err = spectrum_data['SYS_ERR']
    grouping_data = spectrum_data['GROUPING']
    #load the rebinned spectrum 
    group_start = np.where(grouping_data==1)[0]
    total_groups = len(group_start)
    counts_per_group = np.zeros(total_groups,dtype=int)
    bin_bounds_lo = np.zeros(total_groups)
    bin_bounds_hi = np.zeros(total_groups)
    for i in range(total_groups-1):
        counts_per_group[i] = np.sum(counts[group_start[i]:group_start[i+1]])
        bin_bounds_lo[i] = nicer_matrix.emin[group_start[i]]
        bin_bounds_hi[i] = nicer_matrix.emin[group_start[i+1]]
    counts_error = np.sqrt(counts_per_group)
#define bin centers for both spectra 
bound_midpoint = 0.5*(bin_bounds_hi+bin_bounds_lo)
bin_diff = bin_bounds_hi - bin_bounds_lo

baseline_midpoint = 0.5*(nicer_matrix.emin+nicer_matrix.emax)
baseline_width = nicer_matrix.emax-nicer_matrix.emin
baseline_err = np.sqrt(counts)

#plot all 4 and compare 
#note that the files contain counts, so we need to divide by bin width as well 
#as exposure to get the default xspec units of counts/s/keV
plt.figure(figsize=(9,6))

plt.errorbar(xVals,yVals,label='xspec',drawstyle="steps-mid",marker='d')
plt.errorbar(baseline_midpoint,counts/exposure/baseline_width,
             label='astropy',drawstyle="steps-mid",marker='x')
plt.errorbar(bound_midpoint,counts_per_group/exposure/bin_diff,
             label='astropy rebin',drawstyle="steps-mid",marker='*')
plt.errorbar(countsp.energy,countsp.spectrum/bin_width/exposure,
             label='stingray',drawstyle="steps-mid",marker='o')

plt.loglog()
plt.legend(loc='best')
plt.xlabel("Energy bounds (keV)")
plt.ylabel("Counts/s/keV")
plt.xlim([0.5,10.])
plt.ylim([5,2e4])
plt.savefig("Spectrum_comparison.pdf")

Expected Results

In the plot produced at the end, the four spectra should be identical.

Actual Results

Instead, the spectra read with xspec/astropy look identical, but the one loaded with Stingray is offset slightly in some bins and/or shows some features.

[matteobachetti]

@matteolucchini1 thanks for opening the Issue. On Slack we already eliminated some hypotheses, but better continuing the discussion here.
Is it possible to produce the count spectrum with the exact same energy ranges as the pha one and plot a residual (difference and ratio)? This gives the possibility to quantify the difference on a bin-by-bin scale.

[matteolucchini1]

Sure! Following from the code above:

stingray_bins = nicer_rebin.emin[nicer_rebin.emin>0.5]
stingray_bins = stingray_bins[stingray_bins<10.]
stingray_width = np.diff(stingray_bins)

countsp = CountSpectrum(events, energy_spec=stingray_bins)

print(np.allclose(countsp.energy,xVals))

fig, (ax1,ax2,ax3) = plt.subplots(3,1,figsize=(9.,12.),sharex=True,gridspec_kw={'height_ratios': [2, 1, 1]})

ax1.errorbar(countsp.energy,countsp.spectrum/exposure/stingray_width,label='stingray',drawstyle="steps-mid",marker='d')
ax1.errorbar(xVals,yVals,label='xspec',drawstyle="steps-mid",marker='o')
ax1.loglog()
ax1.set_xlabel("Energy bounds (keV)")
ax1.set_ylabel("Counts/s/keV")
ax1.legend(loc='best')
ax1.set_xlim([0.5,10.])

ax2.errorbar(xVals,countsp.spectrum/exposure/stingray_width/yVals)
ax2.set_xscale("log")
ax2.set_xlabel("Energy bounds (keV)")
ax2.set_ylabel("Ratio")
#plt.xlim([0.5,10.])

ax3.errorbar(xVals,countsp.spectrum/exposure/stingray_width-yVals)
ax3.loglog()
ax3.set_xlabel("Energy bounds (keV)")
ax3.set_ylabel("Difference")
ax3.set_xlim([0.5,10.])

fig.tight_layout()
plt.savefig("Stingray_vs_xspec.png")

Which results in the plot attached - there is a systematic offset between the two, as well as some energy-dependent noise. What a mess.
Stingray_vs_xspec.pdf