<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Test Whether Clàrsach Works With Athena X-IFU Simulations
Let's test whether I can make it work on Athena X-IFU simulations. This will be complicated by the fact that sherpa doesn't seem to read the RMF file
Step1:  Let's load some data
Step2:  Let's load the data using Clàrsach
Step3:  Let's also load the ARF and RMF
Step4:  Let's make an empty model to divide out the responses
Step5:  This seems to be working not badly. Let's try the same with sherpa
Step6:  Do the ARF and RMF exist?
Step7:  There's no RMF, because the RMF for Athena does not seem to be a variable length field and is thus not read in. Oh well.
Let's check whether at least the ARF is the same
Step8:  Looks like this worked. Let's do the deconvolution with sherpa and look at the results
Step9:  Well, I didn't actually expect them to be the same, so all right. This might also be due to the fact that I don't understand everything about what ratio_plot is doing; at some point I need to talk to Victoria about that.
Modeling the Spectrum
I can't actually model the full spectrum right now, because I don't have the XSPEC models, thus I don't have the galactic absorption model. This'll have to wait to later. Let's take the second half of the spectrum and fit a model to it. For now, let's ignore the lines and just fit a basic power law
Step10:  We'll need to fix the x_0 parameter of the power law model to continue
Step11:  Let's define a Poisson log-likelihood
Step12:  Ok, cool, let's make a PoissonLikelihood object to use
Step13:  Let's fit this with a minimization algorithm
Step14:  Looks like it has accurately found the photon index of 2. Let's make a best-fit example model and plot the raw spectra
Step15:  Let's also plot the deconvolved version for fun
<ASSISTANT_TASK:>
Python Code:
%matplotlib notebook
import matplotlib.pyplot as plt
try:
    import seaborn as sns
except ImportError:
    print("Seaborn not installed. Oh well.")
    
import numpy as np
import astropy.io.fits as fits
import sherpa.astro.ui as ui
from clarsach.respond import RMF, ARF
Explanation: Test Whether Clàrsach Works With Athena X-IFU Simulations
Let's test whether I can make it work on Athena X-IFU simulations. This will be complicated by the fact that sherpa doesn't seem to read the RMF file:
End of explanation
datadir = "../data/athena/"
data_file = "26.pha"
rmf_file = "athena_xifu_rmf_highres_v20150609.rmf"
arf_file = "athena_xifu_sixte_1469_onaxis_v20150402.arf"
Explanation: Let's load some data:
End of explanation
hdulist = fits.open(datadir+data_file)
hdulist.info()
s = hdulist["SPECTRUM"]
s.columns
channel = s.data.field("CHANNEL")
counts = s.data.field("COUNTS")
plt.figure(figsize=(10,5))
plt.plot(channel, counts)
Explanation: Let's load the data using Clàrsach:
End of explanation
arf = ARF(datadir+arf_file)
rmf = RMF(datadir+rmf_file)
Explanation: Let's also load the ARF and RMF:
End of explanation
resp_model = np.ones_like(counts)
m_arf = arf.apply_arf(resp_model)
m_rmf = rmf.apply_rmf(m_arf)
c_deconv = counts/m_rmf
plt.figure(figsize=(10, 5))
plt.plot(channel, c_deconv)
plt.xscale("log")
Explanation: Let's make an empty model to divide out the responses:
End of explanation
ui.load_data("26", datadir+data_file)
d = ui.get_data("26")
arf_s = d.get_arf()
rmf_s = d.get_rmf()
Explanation: This seems to be working not badly. Let's try the same with sherpa:
End of explanation
print("ARF: " + str(arf_s))
print("RMF: " + str(rmf_s))
Explanation: Do the ARF and RMF exist?
End of explanation
assert np.all(arf_s.specresp == arf.specresp), "Clarsach ARF is different from Sherpa ARF"
Explanation: There's no RMF, because the RMF for Athena does not seem to be a variable length field and is thus not read in. Oh well.
Let's check whether at least the ARF is the same:
End of explanation
ui.set_source("26", ui.polynom1d.truespec)
c_deconv_s = ui.get_ratio_plot("26").y
e_deconv_s = ui.get_ratio_plot("26").x
plt.figure(figsize=(10,5))
plt.plot(e_deconv_s, c_deconv, label="Clarsach Deconvolution")
plt.plot(e_deconv_s, c_deconv_s, label="Sherpa Deconvolution")
plt.legend()
plt.yscale('log')
np.allclose(c_deconv, c_deconv_s)
Explanation: Looks like this worked. Let's do the deconvolution with sherpa and look at the results:
End of explanation
import astropy.modeling.models as models
from astropy.modeling.fitting import _fitter_to_model_params
from scipy.special import gammaln as scipy_gammaln
pl = models.PowerLaw1D()
Explanation: Well, I didn't actually expect them to be the same, so all right. This might also be due to the fact that I don't understand everything about what ratio_plot is doing; at some point I need to talk to Victoria about that.
Modeling the Spectrum
I can't actually model the full spectrum right now, because I don't have the XSPEC models, thus I don't have the galactic absorption model. This'll have to wait to later. Let's take the second half of the spectrum and fit a model to it. For now, let's ignore the lines and just fit a basic power law:
End of explanation
pl.x_0.fixed = True
Explanation: We'll need to fix the x_0 parameter of the power law model to continue:
End of explanation
class PoissonLikelihood(object):
    
    def __init__(self, x, y, model, arf=None, rmf=None, bounds=None):
        self.x = x
        self.y = y
        self.model = model
        self.arf = arf
        self.rmf = rmf
        
        if bounds is None:
            bounds = [self.x[0], self.x[-1]]
        
        min_idx = self.x.searchsorted(bounds[0])
        max_idx = self.x.searchsorted(bounds[1])
        
        self.idx = [min_idx, max_idx]
        
    def evaluate(self, pars):
        # store the new parameters in the model
        _fitter_to_model_params(self.model, pars)
        # evaluate the model at the positions x
        mean_model = self.model(self.x)
        # run the ARF and RMF calculations
        if arf is not None and rmf is not None:
            m_arf = arf.apply_arf(mean_model)
            ymodel = rmf.apply_rmf(m_arf)
        else:
            ymodel = mean_model
        
        # cut out the part of the spectrum that's of interest
        y = self.y[self.idx[0]:self.idx[1]]
        ymodel = ymodel[self.idx[0]:self.idx[1]]
        
        # compute the log-likelihood
        loglike = np.sum(-ymodel + y*np.log(ymodel) \
               - scipy_gammaln(y + 1.))
        if np.isfinite(loglike):
            return loglike
        else:
            return -1.e16
    def __call__(self, pars):
        l = -self.evaluate(pars)
        #print(l)
        return l
Explanation: Let's define a Poisson log-likelihood:
End of explanation
loglike = PoissonLikelihood(e_deconv_s, counts, pl, arf=arf, rmf=rmf, bounds=[1.0, 6.0])
loglike([1.0, 2.0])
Explanation: Ok, cool, let's make a PoissonLikelihood object to use:
End of explanation
from scipy.optimize import minimize
opt = minimize(loglike, [1.0, 1.0])
opt
Explanation: Let's fit this with a minimization algorithm:
End of explanation
_fitter_to_model_params(pl, opt.x)
mean_model = pl(loglike.x)
m_arf = arf.apply_arf(mean_model)
ymodel = rmf.apply_rmf(m_arf)
ymodel_small = ymodel[loglike.idx[0]:loglike.idx[1]]
y_small = loglike.y[loglike.idx[0]:loglike.idx[1]]
e_deconv_small = e_deconv_s[loglike.idx[0]:loglike.idx[1]]
print(np.mean(y_small-ymodel_small))
plt.figure()
plt.plot(e_deconv_small, y_small, label="Data")
plt.plot(e_deconv_small, ymodel_small, label="Model")
plt.legend()
Explanation: Looks like it has accurately found the photon index of 2. Let's make a best-fit example model and plot the raw spectra:
End of explanation
plt.figure(figsize=(10,5))
plt.plot(e_deconv_small, c_deconv[loglike.idx[0]:loglike.idx[1]], label="Data")
plt.plot(e_deconv_small, mean_model[loglike.idx[0]:loglike.idx[1]], label="Best-fit model")
plt.xlim(1.0, 6.0)
plt.legend()
Explanation: Let's also plot the deconvolved version for fun:
End of explanation
<END_TASK>