Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

xsapec abundance not available #15

Open
crhea93 opened this issue Sep 27, 2020 · 4 comments
Open

xsapec abundance not available #15

crhea93 opened this issue Sep 27, 2020 · 4 comments

Comments

@crhea93
Copy link

crhea93 commented Sep 27, 2020

Hello,
There does not appear to be any method for determining the abundance calculated from xspec's apec function. If I set the model to 'xsphabs*xsapec' I am only able to recover xsapec.kT and xsapec.norm, but not xsapec.abundanc despite the fact that I am able to set an initial abundance. Am I missing something?

@DougBurke
Copy link
Collaborator

That doesn't sound right. I'll have to dig into this but I have just had no free time recently (hence the delay in responding).

@DougBurke
Copy link
Collaborator

You are meant to be able to access the xsapec.abundanc values - e.g. https://deproject.readthedocs.io/en/latest/examples/tie.html shows an example with xsmekal, and it should work the same way. Can you show a fragment of code showing your call to get the abundance and the output?

@crhea93
Copy link
Author

crhea93 commented Oct 2, 2020

Thank you for your reply! I thought it would be simple too... Here is an excerpt from the code.

`dep = deproject.Deproject(radii=[float(x) for x in region_vals]* u.arcsec)

# load associated datasets (there are three extracted spectra)

for ann in range(3):
    for obsid in ObsIDs:

        # Take pi file from obsid/repro
        dep.load_pha(obsid+'/repro/'+file_name+'%s.pi' % (ann), annulus=ann)
# Set up fit parameters
dep.set_source('xsphabs*xsapec')
dep.ignore(None, 0.5)
dep.ignore(8, None)
dep.freeze("xsphabs.nh")
dep.set_par('xsapec.redshift', redshift)
dep.set_par('xsphabs.nh', n_H)
dep.subtract()  # Subtract associated background. Read in automatically earlier
onion = dep.fit()
onion_errs = dep.conf()
res = dep.get_fit_results()['xsapec.Abundanc']`

It is in the last line where I try to get the abundance. If I change the line to res = dep.get_fit_results()['xsapec.kT'] for example, I do recover the temperature.

Thank you :)

@DougBurke
Copy link
Collaborator

So I can call

In [19]: dep.get_par('xsapec.abundanc')
Getting xsapec_0.abundanc
Getting xsapec_1.abundanc
Getting xsapec_2.abundanc
Getting xsapec_3.abundanc
Getting xsapec_4.abundanc
Getting xsapec_5.abundanc
Getting xsapec_6.abundanc
Getting xsapec_7.abundanc
Getting xsapec_8.abundanc
Getting xsapec_9.abundanc
Getting xsapec_10.abundanc
Getting xsapec_11.abundanc
Getting xsapec_12.abundanc
Getting xsapec_13.abundanc
Getting xsapec_14.abundanc
Getting xsapec_15.abundanc
Getting xsapec_16.abundanc
Getting xsapec_17.abundanc
Getting xsapec_18.abundanc
Getting xsapec_19.abundanc
Out[19]:
array([0.86248535, 0.99977104, 1.29992861, 1.30085166, 0.97578328,
       1.31974525, 0.95164077, 2.14776034, 0.58024465, 0.67658168,
       ... a bunch of other results

So I can get the values from an APEC model. If I call get_fit_results I get

In [20]: res = dep.get_fit_results()

In [21]: res.keys()
Out[21]:
['annulus',
 'rlo_ang',
 'rhi_ang',
 'rlo_phys',
 'rhi_phys',
 'datasets',
 'succeeded',
 'statval',
 'dstatval',
 'numpoints',
 'dof',
 'qval',
 'rstat',
 'nfev',
 'xsapec.kT',
 'xsapec.Abundanc',
 'xsapec.norm',
 'density']

In [22]: res['xsapec.Abundanc']
Out[22]:
<Column name='xsapec.Abundanc' dtype='float64' length=20>
 0.8624853514583314
 0.9997710377716991
 1.2999286074185443
  1.300851660232163
  0.975783282897546
 1.3197452528254248
 0.9516407743831731
  2.147760341123772
 0.5802446487575588
                ...
 0.6222286676702481
 0.9751073783601172
 0.6288740875191963
 0.6146461546793092
 0.2015569914717479
  0.538681807671724
0.32363535367057805
0.37999792437248636

So it works for me. I'm not sure why it isn't for you.

One thing I notice is that you have multiple datasets per annulus, and I wonder if that has messed things up somehow. If you change you code to the following (ie only have one obsid per annulus) does this work?

for ann in range(3):
    dep.load_pha(ObsIDs[0] +'/repro/'+file_name+'%s.pi' % (ann), annulus=ann)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants