In [ ]:
from pathlib import Path
from glob import glob
import numpy as np
from sherpa.astro import ui
from sherpa.astro import datastack as ds
import sherpa
from astropy.io import fits
from astropy.table import Table
import astropy.units as u
import ChiantiPy.core as ch
In [ ]:
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
# Monkey patch a sherpa bug for loading multiple response files
# see https://cxc.harvard.edu/sherpa/threads/grating_hrcsletg/
import numpy
from sherpa.models.model import CompositeModel
from sherpa.astro.instrument import MultiResponseSumModel
def startup_monkey(self, cache):
pha = self.pha
if numpy.iterable(pha.mask):
pha.notice_response(True)
self.channel = pha.get_noticed_channels()
self.mask = pha.get_mask()
self._get_noticed_energy_list()
CompositeModel.startup(self, cache)
MultiResponseSumModel.startup = startup_monkey
In [ ]:
path = 'data/Chandra/tgcat/obs_6443_tgid_2459/'
ui.load_data(path + 'pha2.gz')
for sign, sherpaid in zip(['-', ''], [1, 2]):
for num in [1,2,3]:
ui.load_arf(sherpaid, path+"leg_{}{}.arf.gz".format(sign, num), num)
ui.load_rmf(sherpaid, path+"leg_{}{}.rmf.gz".format(sign, num), num)
ui.copy_data(1, '6443_leg-1')
ui.copy_data(2, '6443_leg+1')
In [ ]:
# LEG / ACIS
path = 'data/Chandra/tgcat/obs_13250_tgid_3371/'
ui.load_data(path + 'pha2.gz')
for sign, sherpaid in zip(['-', ''], [3, 4]):
ui.load_arf(sherpaid, path+"leg_{}1.arf.gz".format(sign))
ui.load_rmf(sherpaid, path+"leg_{}1.rmf.gz".format(sign))
ui.copy_data(3, '13250_leg-1')
ui.copy_data(4, '13250_leg+1')
In [ ]:
# HEG / ACIS
hegobsids = ['5', '7435', '7436', '7437', '7438']
hegdirs = ['obs_5_tgid_4988', 'obs_7435_tgid_4053', 'obs_7436_tgid_4054', 'obs_7437_tgid_4057', 'obs_7438_tgid_4055']
for h, hd in zip(hegobsids, hegdirs):
path = 'data/Chandra/tgcat/' + hd + '/'
ui.load_data(path + 'pha2.gz')
ui.load_arf(3, path+"heg_-1.arf.gz")
ui.load_rmf(3, path+"heg_-1.rmf.gz")
ui.load_arf(4, path+"heg_1.arf.gz")
ui.load_rmf(4, path+"heg_1.rmf.gz")
ui.load_arf(9, path+"meg_-1.arf.gz")
ui.load_rmf(9, path+"meg_-1.rmf.gz")
ui.load_arf(10, path+"meg_1.arf.gz")
ui.load_rmf(10, path+"meg_1.rmf.gz")
ui.copy_data(3, h + '_heg-1')
ui.copy_data(4, h + '_heg+1')
ui.copy_data(9, h + '_meg-1')
ui.copy_data(10, h + '_meg+1')
In [ ]:
# XMM
path = 'data/XMM/0112880201/pps/'
ui.load_data('XMM_R1', path + 'P0112880201R1S004SRSPEC1003.FTZ')
ui.load_data('XMM_R2', path + 'P0112880201R2S005SRSPEC1003.FTZ')
ui.load_rmf('XMM_R1', path + 'P0112880201R1S004RSPMAT1003.FTZ')
ui.load_rmf('XMM_R2', path + 'P0112880201R2S005RSPMAT1003.FTZ')
In [ ]:
ui.set_analysis('wave')
In [ ]:
ui.ignore(None, 21.4)
ui.ignore(22.3, None)
ui.notice(21.4, 22.3)
In [ ]:
ui.set_stat('cstat')
In [ ]:
bkg = ui.xsbbody(name='bkg')
line_r = ui.delta1d(name='r')
line_i = ui.delta1d(name='i')
line_f = ui.delta1d(name='f')
In [ ]:
bkg = ui.xsconstant(name='bkg')
line_r = ui.xsgaussian(name='r')
line_i = ui.xsgaussian(name='i')
line_f = ui.xsgaussian(name='f')
In [ ]:
line_r.Sigma.val = 0.00001
line_r.Sigma.frozen = True
line_r.LineE = 0.5740
line_r.LineE.frozen = True
line_i.Sigma.val = 0.00001
line_i.Sigma.frozen = True
line_i.LineE = 0.5686
line_i.LineE.frozen = True
line_f.Sigma.val = 0.00001
line_f.Sigma.frozen = True
line_f.LineE = 0.5610
line_f.LineE.frozen = True
In [ ]:
# Ne IX
line_r.LineE = 0.9220
line_i.LineE = 0.9148
line_f.LineE = 0.9051
ui.ignore(None, 13.40)
ui.ignore(13.75, None)
ui.notice(13.40, 13.75)
In [ ]:
idnames = [['5_meg+1', '5_meg-1'],
['XMM_R2'],
['6443_leg+1', '6443_leg-1'],
['7435_meg+1', '7435_meg-1'],
['7436_meg+1', '7436_meg-1'],
['7437_meg+1', '7437_meg-1'],
['7438_meg+1', '7438_meg-1'],
['13250_leg+1', '13250_leg-1'],
]
obsids = [o[0].split('_')[0] for o in idnames]
In [ ]:
def fitlines(names, model):
for n in names:
ui.set_source(n, model)
#ui.group_width(n, 1)
#ui.group_counts(n, 5)
ui.fit(*names)
ui.conf(*names)
cov = ui.get_conf_results()
ui.plot_fit(names[0])
ax = plt.gca()
ax.set_title(names[0])
return np.array([cov.parvals, cov.parmins, cov.parmaxes], dtype=np.floating), (plt.gcf(), ax)
In [ ]:
# Ne IX
line_r.LineE = 0.9220
line_i.LineE = 0.9148
line_f.LineE = 0.9051
ui.ignore(None, 13.40)
ui.ignore(13.75, None)
ui.notice(13.40, 13.75)
fits = []
plots = []
for n in idnames:
f, p = fitlines(n, bkg + line_r + line_i + line_f)
fits.append(f)
plots.append(p)
plt.show()
In [ ]:
farr = np.stack(fits)
farr.shape
In [ ]:
# Ne IX Lyb
line_r.LineE = (11.55 * u.Angstrom).to(u.keV, equivalencies=u.spectral())
ui.ignore(None, 11.45)
ui.ignore(11.65, None)
ui.notice(11.45, 11.65)
fits = []
plots = []
for n in idnames:
f, p = fitlines(n, bkg + line_r)
fits.append(f)
plots.append(p)
plt.show()
In [ ]:
farr2 = np.stack(fits)
farr2.shape
In [ ]:
alpha2beta = farr[:, 0, 1] / farr2[:, 0, 1]
sig_alpha2beta = np.sqrt((farr[:, 2, 1] / farr2[:, 0, 1])**2 + (farr[:, 0, 1] * farr2[:, 2, 1] / farr2[:, 0, 1]**2)**2)
In [ ]:
plt.errorbar(obsids, alpha2beta, yerr=sig_alpha2beta, fmt='o')
In [ ]:
f2i = farr[:, 0, 3] / farr[:, 0, 2]
sig_f2i = np.sqrt((farr[:, 2, 3] / farr[:, 0, 2])**2 + (farr[:, 0, 3] * farr[:, 2, 2] / farr[:, 0, 2]**2)**2)
In [ ]:
fi2r = (farr[:, 0, 3] + farr[:, 0, 2]) / farr[:, 0, 1]
sig_fi2r = np.sqrt((farr[:, 2, 3] / farr[:, 0, 1])**2 +
(farr[:, 2, 2] / farr[:, 0, 1])**2 +
((farr[:, 0, 3] + farr[:, 0, 2]) * farr[:, 2, 1] / farr[:, 0, 1]**2)**2)
In [ ]:
plt.errorbar(obsids, f2i, yerr=sig_f2i, fmt='o')
In [ ]:
plt.errorbar(obsids, fi2r, yerr=sig_fi2r, fmt='o')
In [ ]:
phabs = ui.xsphabs("phabs")
In [ ]:
phabs.nH = 1
print(phabs)
In [ ]:
lam = np.arange(25., 1., -.1)
en = 12.4/lam
phabs.nH = 1
plt.plot(lam, phabs(en), label='$1 \\times 10^{22}$')
phabs.nH = .3
plt.plot(lam, phabs(en), label='$3\\times 10^{21}$')
phabs.nH = .1
plt.plot(lam, phabs(en), label='$1 \\times 10^{21}$')
plt.legend()
In [ ]:
phabs.nH = .1
np.interp([11.55, 13.45], lam[::-1], phabs(en)[::-1])
In [ ]:
phabs.nH = .3
np.interp([11.55, 13.45], lam[::-1], phabs(en)[::-1])
In [ ]:
0.81/0.74, 0.54/.40
In [ ]:
print(phabs)
In [ ]:
logtemp = np.arange(6.1, 7.01, .1)
temp = 10**logtemp
ne9 = ch.ion('ne_9', temperature=temp, eDensity=1.e+9, em=1.e+27)
ne9.intensity()
In [ ]:
# rename this Lya -> Hea to avoid confusion
# Really want Intensity here not Emiss, but for lines of the same ion the result should be the same
e_lyb = ne9.Emiss['emiss'][(ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.3p 1P1.0')].flatten()
e_lya = ne9.Emiss['emiss'][(ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.2p 1P1.0')].flatten()
In [ ]:
((ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.3p 1P1.0')).sum()
In [ ]:
((ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.2p 1P1.0')).sum()
In [ ]:
plt.plot(logtemp, e_lya / e_lyb)
plt.xlabel('log T')
plt.ylabel('ratio Ly$\\alpha$/Ly$\\beta$')
In [ ]:
In [ ]:
e_i1 = ne9.Emiss['emiss'][(ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.2p 3P2.0')].flatten()
e_i2 = ne9.Emiss['emiss'][(ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.2p 3P1.0')].flatten()
e_f = ne9.Emiss['emiss'][(ne9.Emiss['lvl1'] == 1) & (ne9.Emiss['pretty2'] == '1s.2s 3S1.0')].flatten()
In [ ]:
plt.plot(logtemp, (e_i1 + e_i2 + e_f) / e_lya)
plt.xlabel('log T')
plt.ylabel('ratio (f+i)/r')
In [ ]:
xp = (e_i1 + e_i2 + e_f) / e_lya
sortind = np.argsort(xp)
xp = xp[sortind]
fp = logtemp[sortind]
t_from_g = np.interp(fi2r, xp, fp)
t_from_g_up = np.interp(fi2r + sig_fi2r, xp, fp) - t_from_g
t_from_g_do = np.interp(fi2r - sig_fi2r, xp, fp) - t_from_g
In [ ]:
t_from_g_up
In [ ]:
plt.errorbar(fi2r, alpha2beta, xerr=sig_fi2r, yerr=sig_alpha2beta, fmt='o')
In [ ]:
phabs = ui.xsphabs("phabs")
phabs.nH = .1
print(phabs)
In [ ]:
# Input is "edge of bins", return values are for bin center.
# So need to make a few narrow bins around the range I care.
absval = phabs(([13.46, 13.45, 13.44, 11.56, 11.55, 11.54, 5.]
* u.Angstrom).to(u.keV, equivalencies=u.spectral()).value)
absval
In [ ]:
abscoeffalpha = - 1e-21 * np.log(absval[2])
abscoeffbeta = - 1e-21 * np.log(absval[4])
abscoeffalpha, abscoeffbeta
In [ ]:
Gratio, Habobs = np.mgrid[.3:1.4:.01, 3:14:.01]
In [ ]:
logT_from_G = np.interp(Gratio, xp, fp)
Hab_from_logT = np.interp(logT_from_G, logtemp, e_lya / e_lyb)
N_H = np.log(Habobs/Hab_from_logT) / (abscoeffbeta - abscoeffalpha) / 1e21
In [ ]:
# Define colors so that they match figure 3 in Brickhouse et al. (2012)
# Other colors are taken from default matplotlib color cycle
obscolors = {'5': '#ff7f0e',
'XMM': '#9467bd',
'6443': '#8c564b',
'7435': '#d62728',
'7436': '#1f77b4',
'7437': '#2ca02c',
'7438': '#e377c2',
'13250': '#7f7f7f'}
def grat2logt(g):
return np.interp(g, xp, fp)
def logt2grat(logt):
return np.interp(logt, fp[::-1], xp[::-1])
fig, ax = plt.subplots(figsize=(4,3))
cs = ax.contourf(Gratio, Habobs, np.ma.masked_less_equal(N_H, 0), cmap='binary')
cslines = ax.contour(Gratio, Habobs, N_H, levels=[0], linestyles=['dotted'], colors=['k'], linewidths=[4])
#ax.clabel(cs, cs.levels, inline=True, fontsize=10)
cb = fig.colorbar(cs, ax=ax)
for i, o in enumerate(obsids):
eb = ax.errorbar(fi2r[i], alpha2beta[i] * 11.55 / 13.45, xerr=sig_fi2r[i], yerr=sig_alpha2beta[i],
fmt='o', label=o, color=obscolors[o])
secax = ax.secondary_xaxis('top', functions=(grat2logt, logt2grat))
secax.set_xlabel('Temperature [log T in K]')
ax.set_xlim(.4, 1.22)
ax.set_ylim(3, 13)
cb.set_label('$N_H$ [$10^{21}$ cm$^{-2}$]')
ax.set_xlabel('Observed line ratio $(f+i)/r$')
ax.set_ylabel('Observed line ratio He$\\alpha$/He$\\beta$')
ax.legend(loc='upper left')
fig.savefig('/Users/guenther/MITDropbox/my_proposals/Chandra/ULYSSES-TWHya/figs/Ne-var.pdf', bbox_inches='tight')
See https://cxc.cfa.harvard.edu/cal/Letg/LetgHrcEEFRAC/index.html for how to improve the S/N somewhat with non-standard settings. For the proposal I don't need that, fur a publication it might be worthwhile to explore the difference.
Note I don't see changes between chunks, but the wavelength is slightly off (need to fix 0-order location by hand), so I just got the count number from the total 150 ks and devided by five for the proposal.
In [ ]:
import ciao_contrib.runtool as rt
In [ ]:
# Everything could be done with CIAO tools, but I know astropy better, so I migh just as well use that
from astropy.io import fits
from astropy.table import Table
In [ ]:
# From the header for ObsID 6443
# I'm running this on an old computer, thus I don't want to open a large file just to get those two numbers
# when I'm running the notebook again.
TSTART = 272620141.4059900045
TSTOP = 272774659.6507499814
In [ ]:
for i in range(5):
rt.dmcopy.punlearn()
rt.dmcopy(infile='data/Chandra/6443/repro/hrcf06443_repro_evt2.fits[EVENTS][time={}:{}]'.format(TSTART + i * delta_t,
TSTART + (i+1) * delta_t),
outfile='data/Chandra/6443/repro/evt2_{}'.format(i), option="")
rt.dmappend(infile='data/Chandra/6443/repro/hrcf06443_repro_evt2.fits[region][subspace -time]',
outfile='data/Chandra/6443/repro/evt2_{}'.format(i))
In [ ]:
# Assuming DTCOR ppreviously calculated (i.e. for the entire exposure) is OK here, too.
# See https://cxc.cfa.harvard.edu/ciao/threads/spectra_letghrcs/ for how to redo that calculation.
# Not needed for proposal, because effect is in the percent range.
In [ ]:
for i in range(5):
rt.tgextract.punlearn()
rt.tgextract(infile='data/Chandra/6443/repro/evt2_{}'.format(i),
outfile='data/Chandra/6443/repro/pha2_{}'.format(i),
inregion_file='CALDB')
In [ ]:
path = 'data/Chandra/tgcat/obs_6443_tgid_2459/'
for i in range(5):
ui.load_data('data/Chandra/6443/repro/pha2_{}'.format(i))
for sign, sherpaid in zip(['-', ''], [1, 2]):
for num in [1,2,3]:
ui.load_arf(sherpaid, path+"leg_{}{}.arf.gz".format(sign, num), num)
ui.load_rmf(sherpaid, path+"leg_{}{}.rmf.gz".format(sign, num), num)
ui.copy_data(1, 't6443_{}_leg-1'.format(i))
ui.copy_data(2, 't6443_{}_leg+1'.format(i))
In [ ]:
ui.set_analysis('wave')
ui.ignore(None, 13.40)
ui.ignore(13.75, None)
ui.notice(13.40, 13.75)
ui.plot_data('6443_leg-1')
for i in range(5):
ui.plot_data('t6443_{}_leg-1'.format(i), overplot=True)
In [ ]:
ui.set_analysis('6443_leg-1', "energy", "counts", 0)
In [ ]:
ui.plot_data('6443_leg-1')
In [ ]:
pl = ui.get_data_plot('6443_leg-1')
In [ ]:
pl.y.sum()
In [ ]:
ui.set_analysis('wave')
#ui.ignore(None, 21.5)
#ui.ignore(22.3, None)
#ui.notice(21.5, 22.3)
ui.ignore(None, 13.40)
ui.ignore(13.75, None)
ui.notice(13.40, 13.75)
ui.set_analysis('6443_leg+1', "energy", "counts", 0)
pl = ui.get_data_plot('6443_leg+1')
ui.plot_data('6443_leg+1')
print(pl.y.sum())
In [ ]:
5.**0.5
In [ ]: