One task for Athena X-IFU will be measuring the Warm-Hot Intergalactic Medium (WHIM) via absorption spectroscopy. While WHIM itself does not emit X-rays, it will produce absorption lines in the X-ray spectra of bright Gamma-ray Burst (GRB) afterglows. This is important for solving the missing baryon problem.
How well Athena will be able to measure the WHIM depends on how well it is going to be able to identify lines in GRB afterglows. To this purpose, the Athena X-IFU WHIM Data Challenge was called into existence. The point here is figure out the best way to identify lines and their redshift, and find the method that has the highest detection probability for the weak lines expected from WHIM.
Let's first actually look at the data, so that we can set up an appropriate model.
In [1]:
%matplotlib notebook
import matplotlib.pyplot as plt
import seaborn as sns
import glob
import numpy as np
import pandas as pd
import sherpa.astro.ui
import astropy.io.fits as fits
datadir = "../data/"
Before we do anything clever, I want to have a look at the FITS files so that I know what's going on:
In [2]:
import sherpa
In [3]:
sherpa.__version__
Out[3]:
In [4]:
hdulist = fits.open("../data/47.pha")
hdulist
Out[4]:
In [5]:
hdr = hdulist[0].header
In [6]:
hdr
Out[6]:
What's in the data extension?
In [7]:
hdulist[1].header
Out[7]:
In [8]:
hdulist[1].columns
Out[8]:
Hmm, so it doesn't look like I have energy information for each bin. Maybe the ARF/RMF files do?
In [9]:
arf_hdulist = fits.open(datadir+"athena_xifu_sixte_1469_onaxis_v20150402.arf")
In [10]:
arf_hdulist.info()
In [11]:
e_low = arf_hdulist[1].data.field("ENERG_LO")
e_high = arf_hdulist[1].data.field("ENERG_HI")
What about the RMF?
In [12]:
rmf_hdulist = fits.open(datadir+"athena_xifu_rmf_highres_v20150609.rmf")
In [13]:
rmf_hdulist.info()
In [14]:
rmf_hdulist[1].header
Out[14]:
In [15]:
rmf_data = rmf_hdulist[2].data
In [16]:
rmf_data.columns
Out[16]:
In [17]:
datadir = "../data/"
datafiles = glob.glob(datadir+"*.pha")
print(datafiles)
In [18]:
sherpa.astro.ui.load_data("../data/47.pha")
In [19]:
d = sherpa.astro.ui.get_data()
In [20]:
arf = d.get_arf()
In [21]:
arf
Out[21]:
In [22]:
rmf = d.get_rmf()
In [23]:
rmf
In [24]:
for i,f in enumerate(datafiles):
sherpa.astro.ui.load_data(id="%i"%(i+1), filename=datadir+f)
In [25]:
d1 = sherpa.astro.ui.get_data("1")
In [26]:
energy = d1.get_x()
In [27]:
plt.figure()
plt.plot(energy, d1.counts, label="HEG P1", alpha=0.9,
c=sns.color_palette()[0], linestyle="steps-mid")
plt.xlim(0.5, 0.68)
Out[27]:
In the original document, the line search strategy is between 250 and 584 eV for the OVII line and 250 to 754 eV for the OVII and OVIII line. They use Gaussian statistics, i.e. rebin the data, which we are not going to do.
Question: Do I need to worry about over-sampling?
The authors then do a blind fit to the continuum, followed by a blind line search in the range between $0 \leq z \leq 0.1$, considering the simple case of a WHIM line in the nearby IGM. The search is performed using Protassov et al (2002), searching at every energy for a line and producing a posterior predictive p-value for each line energy in its rest frame $E_0$ shifted to $(1+z)E_0$ while oversampling the detector resolution. They produce a new fit for each energy and test for the presence of a line using a significance test.
Question: What are the priors used in the original approach?
We're going to stick this into ShiftyLines and see what happens! But first, let's look at the relevant energy ranges:
In [28]:
len(e_low)
Out[28]:
In [29]:
x_min = 0.25
x_max = 0.754
d = sherpa.astro.ui.get_data("47")
plt.figure()
plt.plot(e_low, d.counts, alpha=0.9,
c=sns.color_palette()[0], linestyle="steps-mid")
plt.xlim(x_min, x_max)
Out[29]:
I have no idea whether there's a line in there. I know that there are different simulations (1) for a single line search versus two line search, (2) different equivalent widths, and (3) different starting fluxes of the GRB. And apparently also for different exposure times.
Let's plot all spectra, just so I know what I'm dealing with, including the two lines of interest:
In [30]:
0.574*(1.0-0.1)
Out[30]:
In [31]:
o_lines = [0.574, 0.654]
In [32]:
plt.figure()
for i,f in enumerate(datafiles):
d = sherpa.astro.ui.get_data("%i"%(i+1))
plt.plot(e_low, d.counts, label="HEG P1", alpha=0.9,
c=sns.color_palette()[0], linestyle="steps-mid")
for l in o_lines:
plt.vlines(l/(1.0 + 0.07), 0, 7000, lw=2, color="black")
plt.xlim(x_min, x_max)
Out[32]:
In [33]:
sample = np.loadtxt("../../atrytone/data/47_posterior_sample.txt")
The $log(Z) = -2580.20885111$.
Let's plot the number of Doppler shifts:
In [34]:
plt.figure()
plt.hist(sample[:,10], bins=50);
In [35]:
plt.figure()
plt.hist(sample[:,11], bins=30)
Out[35]:
In [36]:
x_min = 0.25
x_max = 0.75
In [111]:
min_ind = e_low.searchsorted(x_min)
max_ind = e_low.searchsorted(x_max)
e_low_small = e_low[min_ind:max_ind]
e_high_small = e_high[min_ind:max_ind]
In [38]:
# the data
d = sherpa.astro.ui.get_data("1")
In [39]:
sample = np.loadtxt("../../atrytone/data/01_posterior_sample.txt")
print("There are %i samples in the posterior."%len(sample))
In [40]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
ax1.hist(sample[:,10], bins=50);
ax1.set_title("Number of Doppler Shifts")
ax2.hist(sample[sample[:,11] != 0.0,11], bins=30)
ax2.set_xlabel("Doppler Shift")
ax2.set_ylabel("p(Doppler Shift)")
plt.figure(figsize=(8,4))
plt.plot(e_low, d.counts,
c="black", linestyle="steps-mid")
for s in sample[-20:]:
plt.plot(e_low_small, s[-len(e_low_small)-5:-5], lw=1,
c=sns.color_palette()[0], alpha=0.3)
plt.xlim(x_min, x_max)
plt.ylim(0, 400)
Out[40]:
What's the probability of having a Doppler shift in this source?
In [41]:
len(sample[sample[:,10]!= 0.0,11])/len(sample[:,10])
Out[41]:
What are the amplitudes?
In [42]:
istart = 12
ndoppler = 1
nd1 = len(sample[sample[:,10]!=0.0, 10])
amplitudes = np.zeros((nd1, len(o_lines)))
for i in range(len(o_lines)):
amplitudes[:,i] = sample[sample[:,istart+i*ndoppler]!=0.0, istart+i*ndoppler]
In [43]:
fig, axes = plt.subplots(1,2, figsize=(8,4))
axes = np.hstack(axes)
for i in range(len(o_lines)):
a = amplitudes[:,i]
axes[i].hist(a, bins=30)
axes[i].set_title(o_lines[i])
What are the means and standard deviations of the amplitudes?
In [44]:
for i in range(amplitudes.shape[1]):
print("Mean for amplitude for line at %.5f is %.5f +/- %.5f"%(o_lines[i],
np.mean(amplitudes[:,i]),
np.std(amplitudes[:,i])))
We have to compute the equivalent width of the two lines. For that, we're going to need (1) The Doppler shift of the line, (2) the shifted mean of the line, (3) the width and amplitude of the line, (4) the background
In [45]:
ndoppler = 1
In [46]:
import scipy.special
def gaussian_cdf(x, x0, gamma):
c = 0.5*(1. + scipy.special.erf((x-x0)/(gamma*np.sqrt(2.))))
return c
This is the factor for converting keV into Angstrom
In [47]:
cfac = 12.3984191
In [135]:
import astropy.units as u
In [137]:
k = 42. * u.keV
In [140]:
k.to(u.Angstrom, equivalencies=u.spectral())
Out[140]:
In [207]:
def equivalent_widths(e_low_small, e_high_small, s, o_lines, cfac = 12.3984191):
e_mid_small = e_low_small + (e_high_small - e_low_small)/2.0
e_low_small = e_low_small * u.keV
e_mid_small = e_mid_small * u.keV
e_high_small = e_high_small * u.keV
e_low_ang = e_low_small.to(u.Angstrom, equivalencies=u.spectral())
e_mid_ang = e_mid_small.to(u.Angstrom, equivalencies=u.spectral())
e_high_ang = e_high_small.to(u.Angstrom, equivalencies=u.spectral())
#e_low_ang = np.array([e.to(u.Angstrom, equivalencies=u.spectral()) for e in e_low_small])
#e_mid_ang = np.array([e.to(u.Angstrom, equivalencies=u.spectral()) for e in e_mid_small])
#e_high_ang = np.array([e.to(u.Angstrom, equivalencies=u.spectral()) for e in e_high_small])
# constant background
bkg = s[0]
# threshold parameter for positive/negative sign
thr = s[9]
# number of doppler shifts in this model
nd = s[10]
# value of the doppler shift (will be zero if nd is zero)
dop = s[11]
ew = np.zeros(len(o_lines))
if nd > 0:
lpos_kev = np.zeros(len(o_lines))
lpos_ang = np.zeros(len(o_lines))
amps = np.zeros(len(o_lines))
widths = np.zeros(len(o_lines))
widths_ang = np.zeros(len(o_lines))
sign = np.zeros(len(o_lines))
test = np.zeros(len(o_lines))
lpos_kev = o_lines * u.keV
lpos_ang = lpos_kev.to(u.Angstrom, equivalencies=u.spectral()) * (1. + dop)
for i in range(len(o_lines)):
#lpos[i] = o_lines[i]*1./(1.0+dop)
#lpos_kev[i] = o_lines[i] * u.keV
#lpos_ang[i] = lpost_kev[i].to(u.Angstrom, equivalencies=u.spectral())
amps[i] = np.exp(s[istart+i*ndoppler])
widths[i] = np.exp(s[istart+i*ndoppler+2])
sign[i] = s[istart+i*ndoppler+4]
widths = widths * u.keV
widths_ang = widths.to(u.keV, equivalencies=u.spectral())
for i in range(len(o_lines)):
bkg_model = np.zeros_like(e_mid_ang.value) + bkg
line_model = amps[i]*scipy.stats.norm(lpos_ang[i].value, widths_ang[i].value).pdf(e_mid_ang.value)
total_model = bkg_model + line_model
ew[i] = np.sum((-e_high_ang.value + e_low_ang.value) * np.abs(total_model - bkg_model)/bkg_model)
return ew
# ew[i] = (amps[i]/bkg)*gaussian_cdf(x_max, lpos[i], widths[i])
In [208]:
def equivalent_widths_old(s, o_lines, cfac = 12.3984191, x_max=0.68):
# background flux
bkg = s[0]
# threshold parameter for positive/negative sign
thr = s[9]
# number of doppler shifts in this model
nd = s[10]
# value of the doppler shift (will be zero if nd is zero)
dop = s[11]
ew = np.zeros(len(o_lines))
if nd > 0:
lpos = np.zeros(len(o_lines))
amps = np.zeros(len(o_lines))
widths = np.zeros(len(o_lines))
sign = np.zeros(len(o_lines))
test = np.zeros(len(o_lines))
for i in range(len(o_lines)):
lpos[i] = o_lines[i]*1./(1.0+dop)
amps[i] = np.exp(s[istart+i*ndoppler])
widths[i] = np.exp(s[istart+i*ndoppler+2])
sign[i] = s[istart+i*ndoppler+4]
lpos_angstrom = cfac/lpos
widths_angstrom = cfac/widths
for i in range(len(o_lines)):
ew[i] = (amps[i]/bkg)*gaussian_cdf(x_max, lpos[i], widths[i])
return ew
In [209]:
print(sample[2,10:20])
equivalent_widths(e_low_small, e_high_small, sample[2], o_lines)
Out[209]:
In [210]:
ew = np.array([equivalent_widths(e_low_small, e_high_small, s, o_lines) for s in sample])
ew_old = np.array([equivalent_widths_old(s, o_lines, cfac = 12.3984191, x_max=0.68) for s in sample])
In [211]:
ew
Out[211]:
In [212]:
fig, axes = plt.subplots(1, 2, figsize=(8,4))
for i, l in enumerate(o_lines):
axes[i].hist(np.log(ew[ew[:,i]!= 0.0,i]), bins=30)
print(np.median(ew[ew[:,i]!= 0.0,i]))
axes[i].set_title(l)
In [213]:
for i in range(len(o_lines)):
print("Equivalent width for line %.3f: %.5f +/ %.6f" %(o_lines[i],
np.median(ew[ew[:,i]!= 0.0,i]),
np.std(ew[ew[:,i]!= 0.0,i])))
print("Equivalent width for line %.3f: %.5f +/ %.6f" %(o_lines[i],
np.median(ew_old[ew_old[:,i]!= 0.0,i]),
np.std(ew_old[ew_old[:,i]!= 0.0,i])))
Let's make a function of this:
In [214]:
def plot_ew(sample, o_lines):
# calculate equivalent widths
ew = np.array([equivalent_widths_old(s, o_lines) for s in sample])
fig, axes = plt.subplots(1, 2, figsize=(8,4))
for i, l in enumerate(o_lines):
axes[i].hist(np.log(ew[ew[:,i]!= 0.0,i]), bins=30)
print(np.median(ew[ew[:,i]!= 0.0,i]))
axes[i].set_title(l)
for i in range(len(o_lines)):
print("Equivalent width for line %.3f: %.5f +/ %.6f" %(o_lines[i],
np.median(ew[ew[:,i]!= 0.0,i]),
np.std(ew[ew[:,i]!= 0.0,i])))
return fig, axes
In [215]:
_, _ = plot_ew(sample, o_lines)
Let's make a function of all of it:
In [219]:
def plot_all(nsim, datadir="../data/"):
# the rest energies of the two lines
o_lines = [0.574, 0.654]
# make a string for the correct data file to load
if nsim < 10:
dataid = "0" + str(nsim)
else:
dataid = str(nsim)
datafile = datadir + dataid + ".pha"
# use sherpa to load the data
sherpa.astro.ui.load_data(id="%i"%nsim, filename=datafile)
# get the data
d = sherpa.astro.ui.get_data("%i"%nsim)
# get the posterior sample
sample = np.loadtxt(datadir+"%s_posterior_sample.txt"%dataid)
print("There are %i samples in the posterior."%len(sample))
# compute the fraction of samples that have a Doppler shift
p_doppler = len(sample[sample[:,10]!= 0.0,11])/len(sample[:,10])
doppler_mean = np.mean(sample[sample[:,10]!= 0.0,11])
doppler_std = np.std(sample[sample[:,10]!= 0.0,11])
print("The probability that there is a signal in the data is %.4f"%p_doppler)
print("The estimate for the Doppler shift is %.3f +/- %.3f"%(doppler_mean, doppler_std))
# plot the number of Doppler shifts and the posterior distribution
# of the Doppler shift itself if it's not zero (i.e. doesn't exist)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))
ax1.hist(sample[:,10], bins=50);
ax1.set_title("Number of Doppler Shifts")
ax2.hist(sample[sample[:,11] != 0.0,11], bins=30)
ax2.set_xlabel("Doppler Shift")
ax2.set_ylabel("p(Doppler Shift)")
# plot posterior samples of the light curve
plt.figure(figsize=(8,4))
plt.plot(e_low, d.counts,
c="black", linestyle="steps-mid")
for s in sample:
plt.plot(e_low_small, s[-len(e_low_small)-5:-5], lw=1,
c=sns.color_palette()[0], alpha=0.3)
# set some sane y axis limits.
s = sample[0, -len(e_low_small)-5:-5]
y_min = 0
y_max = np.max(s) + 20.0
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
# let's also plot the posteriors of the amplitudes
# these are the needed indices
istart = 12
ndoppler = 1
# the number of samples where there is a Doppler shift
nd1 = len(sample[sample[:,10]!=0.0, 10])
# empty array for the amplitudes
amplitudes = np.zeros((nd1, len(o_lines)))
# get the amplitudes out of the posterior sample
for i in range(len(o_lines)):
amplitudes[:,i] = sample[sample[:,istart+i*ndoppler]!=0.0, istart+i*ndoppler]
# make a plot of the amplitudes
fig, axes = plt.subplots(1,2, figsize=(8,4))
axes = np.hstack(axes)
for i in range(len(o_lines)):
a = amplitudes[:,i]
axes[i].hist(a, bins=30)
axes[i].set_title(o_lines[i])
# print the mean and standard deviation of the amplitudes.
for i in range(amplitudes.shape[1]):
print("Mean for amplitude for line at %.5f is %.5f +/- %.5f"%(o_lines[i],
np.mean(amplitudes[:,i]),
np.std(amplitudes[:,i])))
# compute the equivalent widths
ew = np.array([equivalent_widths_old(s, o_lines) for s in sample])
# plot the equivalent widths
_, _ = plot_ew(sample, o_lines)
return
In [220]:
plot_all(1)
In [221]:
plot_all(2)
In [222]:
plot_all(3)
In [223]:
plot_all(4)
In [224]:
plot_all(5)
In [226]:
plot_all(6)
In [227]:
plot_all(7)