plot the data in $\rm counts~s^{-1}~cm^{-2}~keV^{-1}$

[rhuang-astro]

Hi,

Is it possible to plot the spectrum in $\rm counts~s^{-1}~cm^{-2}~keV^{-1}$?

The result of ui.plot_data(1) is $\rm counts~s^{-1}~keV^{-1}$.

It is straight forward to do so in xspec with setplot area.

Best, Rui

[DougBurke]

So, XSPEC's setplot area command is described as

After setplot area is entered, plot data and plot ldata will show the data divided by the response effective area for each particular channel. plot residuals will necessarily also be affected by this. Usual plotting is restored by setplot noarea. If data is associated with more than 1 response, the response effective area is calculated by simply summing the contributions from each response.

We have the information to do this, but at the moment it would have to be a separate routine, not part of Sherpa. We are planning to look at Sherpa plots so we can add this to the list.

For now, it'll take some time to come up with a solution to

a) check how XSPEC deals with grouping b) check how XSPEC handles the mapping between channel and energy (does it just use the EBOUNDS block of the RMF?)

[rhuang-astro]

Thank you!

[hamogu]

OK, you did not say way exactly the purpose of this plot is, so it's entirely possibly that all I'm going to say below does not apply to your use case. First, I'm going to explain why plots like this are often misleading so that people who find this issue later with Google understand; this may or may not apply to you. Then, I'll show you how to do what you ask for to answer the question in the title in a separate post.

There is a reason why a plot like this it's not easy: We believe that plots that show the X-ray spectrum divided by the effective area are often misleading.

Here are a few thoughts that show the problems:

• If you have binned data with large energy bins, which ARF value do you use? Do you average over a bin? What if the spectrum changes in the bin, e.g. it has a spectral line or a steep powerlaw? Do you use a weighted average, which gives more weight to e.g. the ARF at the left side of a bin? And if so, do you use the observed, noisy spectru for the weight or a model (and then: which model)?
• In X-ray astronomy, we often have a significant redistribution (RMF), so that photons are not detected at the right energy because of the limited energy resolution of the CCDs. So, a photon that has 1 keV is detected somewhere between 0.9 and 1.1 keV. Now imagine a detector that has a specific energy where the effective area is very low. That's not uncommon: See e.g. Fig 6.5 here: https://cxc.harvard.edu/proposer/POG/html/chap6.html and note the drop at 1.9 keV in the dotted line for the ACIS-I3 chip. That's due to absorption by one of the components that ACIS is made out of (on top of my hat, I think that one is silicon, but I would have to check). In an observed spectrum plotted with ui.plot_data(1) you don't see that. Because the ACIS energy resolution is lower, some of the photons with energies of e.g. 1.8 keV or 2 keV are also detected at 1.9 keV, so it all looks smooth. If you now take that smooth, observed spectrum and devide it by the ARF like XPSEC does, then you devide the spectrum at 1.9 by a small number and the result will look af if there is a 1.9 keV emission line. That's wrong!

[hamogu]

I see that @DougBurke outtyped me by 10 min... I hope my answer is still useful.

[hamogu]

That said, here is how you could do it. In this case, I've made the choice to user interpolation to find the value of the ARF in the middle of the bin (which might work for an ARF that is smooth, but then again, I have to make some simplification to do what you asked for):

import os
import numpy as np
import matplotlib.pyplot as plt
import sherpa
import sherpa.astro.ui as ui

%matplotlib inline

datafile = os.path.dirname(sherpa.__file__) + '/astro/datastack/tests/data/3c273cp.pi'
ui.load_pha(datafile, use_errors=True)
arf = ui.get_arf()
dataplot = ui.get_data_plot()

approx_arf = np.interp(dataplot.x, arf.get_x(), arf.get_y())
plt.errorbar(dataplot.x,
             dataplot.y / approx_arf,
             yerr=dataplot.yerr / approx_arf,
            marker='o')
plt.xlabel(r'Energy (keV)')
plt.ylabel(r'$\rm counts~s^{-1}~cm^{-2}~keV^{-1}$')

That will give you this plot:

And there you see the problems I was describing: The spectrum seems very, very bright at low energies - in fact so bright that you can't see what is going on at energy > 0.2 keV. However, that does not mean that the source really is that bright - it just means that the ARF is very low in this region and we divide the few photons be have by a very small number - leading to a very large datapoint in the plot.

If I zoom in more, then I see that the rest of the plot looks fine:

But, still, there are features in the ARF. Just as I said in my comment above, there is an apparent emission line at 1.9 keV, which is not real emission, it's just that the ARF has a lot of structure in this region, which together with the fast that photons are not always detected at the real energy can cause these apparent feature.

Here is a picture of the ARF for this dataset:

[hamogu]

So, personally, I quite like that it's hard to make this type of plot in Sherpa; making this type of plot forces you to think about the choices that go into it.

However, if users need it, we should add a function to make it easier to plot that.

@rhuang-astro It would be great if you could let us know your thoughts about this? After all the discussion, do you feel we should make a function to make it easier to make a plot like this? Is there anything else that is missing in Sherpa plotting?

[rhuang-astro]

I understand your concern.

We deal with spectrum in this way: model ( $\rm ph~s^{-1}~cm^{-2}~keV^{-1}$ ) $\times$ ARF ( $\rm cm^2$ ) $\times$ RMF ( $\rm keV ~ channel^{-1}$ ) = data ( $\rm count ~ s^{-1} ~ channel^{-1}$). It's not easy to restore the $\rm cm^2 ~ keV$ directly since the RMF is a matrix. In principle, we should plot the spectrum in $\rm count ~ s^{-1} ~ channel^{-1}$, that's how it just binned. However, we still often use $\rm count ~ s^{-1} ~ keV^{-1}$.
This is the trade-off between correctness and intuition. Such desire for plotting (e.g. plot eeufspec ) is only for plotting, we always fit the data in the correct way.

[rhuang-astro]

What a high efficiency @hamogu , at least the code can be used in some cases.

[rhuang-astro]

Sometimes people also like ARF ( $\rm cm^2$ ) $\times$ RMF ( $\rm keV ~ channel^{-1}$ ) = RSP ( $\rm cm^{2} ~ keV ~ channel^{-1}$ ). If we have already turn the $\rm count ~ s^{-1} ~ channel^{-1}$ to $\rm count ~ s^{-1} ~ keV^{-1}$, it should be the same for $\rm count ~ s^{-1} ~ cm^{-2} ~ keV^{-1}$.
So let me guess, we just use the E-C relation stored in EBOUNDS to make the conversion, so should the effective area.

[DougBurke]

For reference, here are the XSPEC plots I get for our 3c273.pi dataset:

% xspec

        XSPEC version: 12.12.1
    Build Date/Time: Thu Apr 28 18:10:42 2022

XSPEC12>cpd /xw
XSPEC12>data 3c273.pi

1 spectrum  in use

Spectral Data File: 3c273.pi  Spectrum 1
Net count rate (cts/s) for Spectrum:1  1.833e-02 +/- 7.054e-04 (96.0 % total)
 Assigned to Data Group 1 and Plot Group 1
  Noticed Channels:  1-46
  Telescope: CHANDRA Instrument: ACIS  Channel Type: PI
  Exposure Time: 3.856e+04 sec
 Using fit statistic: chi
 Using Background File                3c273_bg.pi
  Background Exposure Time: 3.856e+04 sec
 Using Response (RMF) File            3c273.rmf for Source 1
 Using Auxiliary Response (ARF) File  3c273.arf

XSPEC12>setplot energy
XSPEC12>ignore **-**
    46 channels (1,46) ignored in spectrum #     1

XSPEC12>notice 0.5-7.0
    42 channels (4-45) noticed in spectrum #     1

XSPEC12>plot data
XSPEC12>iplot data
PLT> wdata before.txt
PLT> quit

and then

XSPEC12>setplot area
XSPEC> plot data
XSPEC12>iplot data
PLT> wdata after.txt
PLT> quit

and the output files are

before.txt

% cat before.txt 
READ SERR 1 2
@before.pco
!
0.518299997 5.11000007E-2 3.77161545E-3 9.83261736E-4
0.605900049 3.65000069E-2 5.63547295E-3 1.42165762E-3
0.671599984 2.91999876E-2 6.66023651E-3 1.71966571E-3
0.7227 2.18999982E-2 1.06563745E-2 2.51173158E-3
0.766499996 2.18999982E-2 1.06563745E-2 2.51173158E-3
0.802999973 1.46000087E-2 1.32006463E-2 3.4414141E-3
0.83950001 2.18999982E-2 1.05764987E-2 2.51300144E-3
0.876000047 1.46000087E-2 1.332046E-2 3.43932793E-3
0.919800043 2.91999876E-2 6.48051593E-3 1.72279321E-3
0.970899999 2.18999982E-2 1.116852E-2 2.58179498E-3
1.01469994 2.18999982E-2 8.80043674E-3 2.29427777E-3
1.06579995 2.92000175E-2 6.60032313E-3 1.72070705E-3
1.11689997 2.18999982E-2 1.00643542E-2 2.44096434E-3
1.16799998 2.91999578E-2 7.10425945E-3 1.77606486E-3
1.24099994 4.3800056E-2 4.65628458E-3 1.18538644E-3
1.34319997 5.83999753E-2 3.71422712E-3 9.16341611E-4
1.43809998 3.64999771E-2 5.32819051E-3 1.37573283E-3
1.54029989 6.56999946E-2 3.66958953E-3 8.61421635E-4
1.64980006 4.3800056E-2 4.40021232E-3 1.14713726E-3
1.75200009 5.83999753E-2 3.46226594E-3 8.89545772E-4
1.85420001 4.37999964E-2 4.40021837E-3 1.14713889E-3
1.91989994 2.18999982E-2 9.47233289E-3 2.36808322E-3
1.98559999 4.3800056E-2 4.91235685E-3 1.22243934E-3
2.05859995 2.91999578E-2 6.6602435E-3 1.71966746E-3
2.13890004 5.11000156E-2 4.46431851E-3 1.07808772E-3
2.23379993 4.37999964E-2 4.73616645E-3 1.18404161E-3
2.33599997 5.83999157E-2 3.33012175E-3 8.5983373E-4
2.48929977 9.4900012E-2 1.97557104E-3 5.30410325E-4
2.64989996 6.57000542E-2 3.1574415E-3 7.89360376E-4
2.78859997 7.29999542E-2 2.59220693E-3 6.89117471E-4
2.97109985 0.109500051 1.72813609E-3 4.59411123E-4
3.23390007 0.153300047 1.22297229E-3 3.28349182E-4
3.47480011 8.75999928E-2 2.34811427E-3 5.92357537E-4
3.67920017 0.116799951 1.73113297E-3 4.44772886E-4
3.91280007 0.11680007 1.59017392E-3 4.31218359E-4
4.13910007 0.109499931 1.69618754E-3 4.59966803E-4
4.48219967 0.233599901 8.50589713E-4 2.22638453E-4
4.86909962 0.153300047 1.28472503E-3 3.39450198E-4
5.19759989 0.175199986 1.06011669E-3 2.87479081E-4
5.63560009 0.262799978 7.49423169E-4 1.98012684E-4
6.2342 0.335800171 5.42685448E-4 1.50169813E-4
8.2198 1.64980006 9.77345699E-5 3.07854498E-5

after

% cat after.txt 
READ SERR 1 2
@after.pco
!
0.518299997 5.11000007E-2 3.8153425E-4 9.94661459E-5
0.605900049 3.65000069E-2 3.98739387E-4 1.00589765E-4
0.671599984 2.91999876E-2 3.17908707E-4 8.20836722E-5
0.7227 2.18999982E-2 3.9939268E-4 9.41377584E-5
0.766499996 2.18999982E-2 3.3560631E-4 7.91031634E-5
0.802999973 1.46000087E-2 3.67285858E-4 9.57515804E-5
0.83950001 2.18999982E-2 2.64582312E-4 6.2865387E-5
0.876000047 1.46000087E-2 3.04721791E-4 7.86788296E-5
0.919800043 2.91999876E-2 1.36434726E-4 3.62700812E-5
0.970899999 2.18999982E-2 2.18927351E-4 5.06088109E-5
1.01469994 2.18999982E-2 1.66690108E-4 4.34561844E-5
1.06579995 2.92000175E-2 1.22380938E-4 3.19047649E-5
1.11689997 2.18999982E-2 1.86378646E-4 4.52034619E-5
1.16799998 2.91999578E-2 1.33637135E-4 3.34092838E-5
1.24099994 4.3800056E-2 9.06178975E-5 2.30693004E-5
1.34319997 5.83999753E-2 7.51089319E-5 1.85302179E-5
1.43809998 3.64999771E-2 1.0900855E-4 2.8145887E-5
1.54029989 6.56999946E-2 7.44015342E-5 1.74654651E-5
1.64980006 4.3800056E-2 8.03869552E-5 2.09569134E-5
1.75200009 5.83999753E-2 4.97635519E-5 1.27855446E-5
1.85420001 4.37999964E-2 5.17818589E-5 1.34995526E-5
1.91989994 2.18999982E-2 9.41539765E-5 2.35384941E-5
1.98559999 4.3800056E-2 4.42773235E-5 1.10184064E-5
2.05859995 2.91999578E-2 6.89081498E-5 1.77920083E-5
2.13890004 5.11000156E-2 5.19338282E-5 1.25414936E-5
2.23379993 4.37999964E-2 4.66068268E-5 1.16517067E-5
2.33599997 5.83999157E-2 3.80757774E-5 9.83112386E-6
2.48929977 9.4900012E-2 2.67143605E-5 7.17239345E-6
2.64989996 6.57000542E-2 4.10687717E-5 1.02671929E-5
2.78859997 7.29999542E-2 3.24075736E-5 8.61529406E-6
2.97109985 0.109500051 2.09977534E-5 5.58208467E-6
3.23390007 0.153300047 1.32659798E-5 3.56171108E-6
3.47480011 8.75999928E-2 2.29117395E-5 5.77993205E-6
3.67920017 0.116799951 1.54334848E-5 3.96526184E-6
3.91280007 0.11680007 1.28723614E-5 3.49068637E-6
4.13910007 0.109499931 1.27448902E-5 3.4561192E-6
4.48219967 0.233599901 5.92840343E-6 1.55173586E-6
4.86909962 0.153300047 8.77426919E-6 2.31833837E-6
5.19759989 0.175199986 7.61328101E-6 2.06454547E-6
5.63560009 0.262799978 6.19864295E-6 1.63780612E-6
6.2342 0.335800171 5.68464657E-6 1.5730335E-6
8.2198 1.64980006 3.36574021E-6 1.0601758E-6

I can get the following which is close, but not quite what XSPEC gives (the error bars aren't quite the same, but that's not what I'm looking at here):

In [227]: for xx, yy in zip(plot.x, plot.y):
     ...:     print(f"{xx:.5f} {yy:.5e}")
     ...: 
0.51830 3.92170e-04
0.60590 4.02419e-04
0.67160 3.21622e-04
0.72270 3.97637e-04
0.76650 3.42294e-04
0.80300 3.70239e-04
0.83950 2.66789e-04
0.87600 3.03820e-04
0.91980 1.41679e-04
0.97090 2.20512e-04
1.01470 1.68118e-04
1.06580 1.23484e-04
1.11690 1.86535e-04
1.16800 1.33266e-04
1.24100 9.22264e-05
1.34320 7.64389e-05
1.43810 1.08926e-04
1.54030 7.58734e-05
1.64980 8.20311e-05
1.75200 4.94513e-05
1.85420 5.23425e-05
1.91990 9.42490e-05
1.98560 4.54554e-05
2.05860 6.80583e-05
2.13890 5.22854e-05
2.23380 4.65749e-05
2.33600 3.72515e-05
2.48930 2.78854e-05
2.64990 4.13662e-05
2.78860 3.34691e-05
2.97110 2.13197e-05
3.23390 1.37639e-05
3.47480 2.30592e-05
3.67920 1.55699e-05
3.91280 1.33310e-05
4.13910 1.33650e-05
4.48220 6.19793e-06
4.86910 9.24018e-06
5.19760 8.08975e-06
5.63560 6.54390e-06
6.23420 6.31209e-06
8.21980 3.31392e-06

[DougBurke]

You can try this

• VERY EXPERIMENTAL
• NOT-WELL-TESTED
• NOT-GUARANTEED
• NOT-TESTED-AGAINST-XSPEC

bit of code (unfortunately GH don't make it possible to attach a .py file hence the odd name).

areaplot.py.txt

If saved as areaplot.py (the name doesn't matter) then you can try

% sherpa areaplot.py

and, as well as a warning frmo traitlets.py which we can't do much about, you should also see

To turn on 'setplot noarea' behaviour, call

> set_hack_area(True)

and to turn it off,

> set_hack_area(False)

NOTE: This changes plot_data, plot_model, and plot_fit.
      It does *not* change plot_resid.

WARNING
  This is an approximation, and it is not a well-tested
  approximation at that (the approximation for the error
  term is 'even worse' than for the data values).

The warning messages are controlled by

> set_warn_area(False)
> set_warn_area(True)

and have been left 'on' just to help you see when the
code is in use.

This should change plot_data, plot_model, plot_fit calls when you've called set_hack_area(True).

The code will note when it's doing things, because I want to make sure it isn't accidentally changing things it shouldn't (I'd be interested to know if you find any like this). You can turn this off with set_warn_area(False) but I suggest leaving it on until you're happy.

Have I mentioned this is COMPLETELY UNTESTED. Although for the one test case I did

> set_hack_area(True)
> plot_fit(ylog=True)
> plot_source(overlpot=True)

did show "okay" agreement (I am not qualifying what I mean by okay).

You should see similar results to the code that @hamogu shows except for areas where the ARF changes significantly across a group; in this case I'd expect this version to be "better behaved", but it is still not the "true" value, as we've mentioned above.

[rhuang-astro]

This is awesome! Thank you very much.

[DougBurke]

Here is a tweaked version (to better support data which has a RSP, so no ARF, just a RMF) - if the previous code was working for you then this should not change it (but I can't guarantee things!).

areaplot.py.txt