Add flux spectra #15
Add flux spectra #15
Conversation
…into add_flux_spectra
…ctionality to plot method in XSpectrum class
|
|
||
| # make a flat spectrum so that I can integrate | ||
| # ARF and RMF only | ||
| flat_model = np.ones_like(spec.counts) * (spec.rmf.energ_hi - |
There was a problem hiding this comment.
This is what I'm working on in xpysis.
I am certain that you should not multiply by (energ_hi - energ_lo). When I do this with the test data on Mrk 421, it gives me something that I know is an order of magnitude off from the true values
There was a problem hiding this comment.
To see how this method for unfolding data performs, see this notebook:
https://github.com/eblur/xpysis/blob/master/notebooks/Example1_fit_HETG_powerlaw.ipynb
There was a problem hiding this comment.
Regarding the multiplication by (energ_hi - energ_lo): I don't know how to square that circle, because I have a case where the exact opposite is true: if I don't multiply by the bin width, I end up with something that is no longer a power law, because my power law and the fact that the energy bins are logarithmic messes up the slope.
I don't know how to make sense of that.
There was a problem hiding this comment.
Oh my. Is it instrument dependent? What data set are you using where you do need to multiply by the bin width?
There was a problem hiding this comment.
Well, right now simulated Chandra/HETG data, actually, so this is the same instrument.
Here's a test I just did: I took the simulated Chandra/HETG spectrum (the one that's just a power law), and applied both methods to it (once without multiplying by the energy bin, and once with), and then plotted that with a comparison to the code that the sherpa manual suggests using.
Here's the plot:
I can see that this is about an order of magnitude off from what sherpa calculates (and sherpa doesn't have those weird features), but the one where I don't multiply by the bin width is off, too, and it has a different slope. So, if sherpa is right, currently neither of our methods is working correctly. If sherpa is wrong, then I don't know at all anymore what's going on.
Victoria and I tested this by eye when I was at ESTEC: the code I used (which is similar to what I implemented in this PR) seemed to give the right answer (where "right" == "the same as ISIS").
I am super confused.
There was a problem hiding this comment.
Those features (e.g. the bump around 2 keV) should definitely not be there. Sherpa is right in this case.
If the fitting is still working (i.e. giving us a power law with photon index = 2), then I'm inclined to say it's our method of unfolding the data. If the fitting is not working any more, then it's something wrong with the way we are applying the response.
Thank you for testing this, as I have never been able to install sherpa properly and can only use ISIS as a comparison.
There was a problem hiding this comment.
By the way, did you notice PR #16 ? If you have not yet reviewed this change, this could be the problem.
(but even after fixing it in my own code, I still get those weird features around 2 keV, which are clearly artifacts from the ARF)
There was a problem hiding this comment.
Well, so I'm pretty sure it's the unfolding, not the application of the ARF/RMF, because I tested those on simulated spectra with a direct comparison to sherpa, and that gives me correct results for all instruments I tested so far.
Here's a slide I found in a etalk by David Huenemoerder that deals with this:
I find that incredibly confusing, though. The words sound like we should be integrating: if S(E) = 1, then we want the counts per unit flux, then we should be integrating S(E) in order to get something that's in the same units as the data. But the equation makes it look like we shouldn't be integrating. He does mention, though, that you can find spurious features. I've been trying to figure out how exactly sherpa computes this plot, but so far, I've not managed to track it down within the source code.
Do you think it would help if we talked this through per telecon rather that GitHub issues next week?
There was a problem hiding this comment.
Yes, let's make some time to skype about this. I think it's easier in person. Plus we may as well talk about some other considered changes. I am going home soon, so next week is good. I'll email you on Monday.
There was a problem hiding this comment.
Okay, I will go home soon, too, but first, two observations:
(1) The weird features in the spectrum are caused by applying frac_expo. I don't know why, but when I don't apply frac_expo, they go away.
(2) It looks like there should definitely not be a np.sum in the expression for the RMF calculation. The reference implementation that I used to compare with some of the simulated data sets (NICER, HEXTE, Chandra etc) doesn't have a sum.
If those two conditions are fulfilled, I get the same spectrum as sherpa does.
| spec.rmf.energ_lo) | ||
|
|
||
| # now apply ARF to the flat model | ||
| m_arf = spec.arf.apply_arf(flat_model) |
There was a problem hiding this comment.
Use spec.apply_resp(flat_model) here instead of apply_arf and apply_resp.
apply_resp takes into account fracexpo. It is not specific to Chandra. For ARF files where the "FRACEXPO" column does not exist, fracexpo is set to 1.0 and doesn't really do anything.
|
|
||
| # divide the observed counts by the flat model with the ARF/RMF applied | ||
| if counts is None: | ||
| flux_spec = spec.counts / m_rmf / spec.exposure |
There was a problem hiding this comment.
If apply_resp is used, you will not need the spec.exposure term here.
I keep getting asked to give people spectra in flux units rather than counts, so I stuck it into Clarsach.
Any comments welcome.