Open
Description
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.