Download the notebook.
Astrosprint: X-ray tutorial¶
- Work through one example
- won't explain every step
- Lots of documentation available on the internet. For Chandra data, look here:
- Homepage: https://cxc.cfa.harvard.edu/ciao/
- examples similar to what I show now: https://cxc.cfa.harvard.edu/ciao/threads/index.html
- help for individual CIAO commands: https://cxc.cfa.harvard.edu/ciao/ahelp/index_alphabet.html
Using CIAO¶
Two ways of using CIAO:
- On the command line (that's what most the documenation looks like)
- through a Python interface (that's what I'll do here, because it's easy to show in a notebook)
In the following, commands that can be done either at the terminal or in Python are in two consequtive cells. Either copy the comment in the cell that says TERMINAL into the terminal OR run the Python version below but not both.
Choose a dataset to work with¶
pick an ObsID from the list¶
I've selected a set of different datasets that are all taken in the same mode and look at similar targets. That way, we can all follow the same steps and use similar models for analysis.
Select an ObsID from https://docs.google.com/spreadsheets/d/167MmraWsYMN7fujjxG1-WpfvtXYALxTCIeLDqSLbZH0/edit?usp=sharing and put your name in. It's OK if more than one person works on the same dataset.
In independent analysis, how do you find Chandra data?¶
Looks at Chandra data search page: https://cda.harvard.edu/chaser/ and filter by Ra/Dec, time, instrument etc.
Set up for Python users¶
from ciao_contrib.cda.data import download_chandra_obsids
from ciao_contrib import runtool as rt
Download and extract that data¶
The steps in here take some time (15 min on my laptop), so you want to run them only once. The results (raw data and extracted spectra) are saved on your hard drive. If stop and start the notebook agian, you could just skip those cells, I can just skip those steps in the notebook.
# Can give multiple ObsIDs in a list, but we do only one right now.
# Put in your dataset here
download_chandra_obsids([26133])
# True means the download worked
Download the required calibration files¶
This is ONLY needed if you did not install the caldb_main package. If you did install the full CALDB skip the next few cells!
CALDB is the Chandra calibration database. It holds all calibration information for all Chandra observations ever done. However, many files are only valid for a certain period, e.g. for a specific year. So, if the full CALDB is not installed, then we need to download the files relevant (e.g. for the year this observations was taken) now.
rt.download_obsid_caldb("26133")
Now, here is a bit of a disadvantage of using Python. We need to set the environment variables that tell CIAO where to find the CALDB. This task writes a nice script, to do that, but that script won't run in Python. So, need to do some handwork here. However, the path to where it all went is printed in the output above, so we can do that by hand.
import os
os.environ['CALDB'] = '/Users/guenther/mambaforge/envs/india_workshop/CALDB'
os.environ['CALDBCONFIG'] = '/Users/guenther/mambaforge/envs/india_workshop/CALDB/software/tools/caldb.config'
os.environ['CALDBALIAS'] = 'none'
Reprocess data¶
Always reprocess the data to apply the newest calibration. Also, this data is taken in VFAINT mode, which allows us to filter out background, that the default processing cannot remove from the event list.
# reprocess the data with the newest calibration files.
# Not really needed for data taken in the last year,
# but just do it here so we won't forget if we ever deal with older data
rt.chandra_repro("26133", # name of the directory
check_vf_pha=True) # VFAINT extra filtering
# This step makes some intermediate files
# If you want to run it a second time, you need to remove those files first.
Now, we need to select the region on the detector we want to extract. There are several ways to do that -and even ways to select the appropriate regions autmomatically- but the simplist way of doing it is to open the evt2 files in ds9 like so: ds9 26133/repro/*evt2* -log -cmap heat &. ds9 is included in the CIAO installation so if you activate that conda environment, you would be able to simply type ds9. In ds9, you can move around and click on the buttons with your mouse. You then look for the obvious brihgt source and look up the "physical coordinates" (= the pixels on the camera). In this case, the source is at 4096, 4088.5, so we put that in the infile parameter below with a radius of 2 pixels. Then, we look for a region with no sources, where we can measure the background. We want a big region (so we get a good statistical sample). I selected a circle centered on 3700, 4200, with a radius of 300.
I can explain this more in person / on zoom. Once you know what to look for it's really simple.
Use ds9 to look at your data ...¶
Start ds9 from command line or from notebook. If you start it from notebook, cells will be frozen, until you close ds9 again!
!ds9 26133/repro/*evt2* -log -cmap heat
... or find the bright sources using a tool¶
rt.dmcopy("26133/repro/acisf26133_repro_evt2.fits[energy=1000:5000][bin x=3000:5000:4,y=3000:5000:4]",
"image.fits", clobber=True)
rt.mkpsfmap("image.fits", "mypsfmap.fits", energy=1.49, ecf=0.393)
rt.wavdetect("image.fits", "source_list.fits", "source_cell.fits",
"srcimage.fits", "background.fits", psffile="mypsfmap.fits")
# We can also use "!" to run programs on the console instead of using the rt.dmllist syntax
!dmlist "source_list.fits[cols RA,dec,pos,net_counts]" data
Check the lightcurve¶
I would start from here when I run the notebook a second time and just want to tweak the fit a little or make a new plot.
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
rt.dmextract(infile=f"26133/repro/acisf26133_repro_evt2.fits[sky=circle(4108.9,4083.8,4)][bin time=::1000]",
outfile=f"lightcurve.fits", opt="ltc1", clobber=True)
# I installed astropy because I'm more used to it, but it's not included in CIAO by default.
# Simply "pip install astropy" AFTER activating your CIAO conda environment should do the trick.
from astropy.table import Table
lc = Table.read(f"lightcurve.fits", hdu=1)
lc
plt.errorbar(lc['TIME'], lc['COUNT_RATE'], yerr=lc['COUNT_RATE_ERR'])
The lightcurve has some variability, maybe that's even ust noise, but certainly not a major flare. So, we continue to observe it in one piece. If there was a flar, we might want to extract two spectra (one befor ethe flar and one during the flare), but here, we don't need to do any of that.
Extract the spectrum¶
# reset specextract to default settings because you never know
# what old version of settings might be lurking somewhere in your file system.
rt.specextract.punlearn()
rt.specextract(infile=f"26133/repro/acisf26133_repro_evt2.fits[sky=circle(4108.9,4083.8,4)]",
outroot=f'src',
bkgfile=f"26133/repro/acisf26133_repro_evt2.fits[sky=circle(3700,4200,300)]",
grouptype="NONE", binspec="NONE", weight=False, correctpsf=True, clobber=True
)
Fit the spectrum¶
# Set better defaults for log plots with a small (< 2 mag) range
mpl.rcParams['axes.formatter.min_exponent'] = 2
from sherpa.astro import ui
ui.load_data(26133, 'src.pi') # 26133 is the index of our data. I use the ObsID because we probably want
# to load data from different OBsIDs to compare at some point, but I could
# call it 1 or "asefsdjf" instead.
# Sherpa can load different datasets at the same time, so in the following
# we give this ID when we plot, fit etc, to say which dataset we use.
This is an example of the Sherpa magic. Underneath there are "data objects" (like a table, an array, etc.), but instead of remembering how all those are called, sherpa.ui keeps an internal list and we just type the ID for poltting, fitting, etc. We could use the full objects instead, but a simple ID is faster to type, so we do that here.
ui.plot_data(26133)
We see that
- There are very few counts at high energies.
- Some of the bins just don't have much data.
ui.group_counts(26133, 10)
ui.plot_data(26133)
Better. Still, Chandra is not sensitive > 8 keV or so, everything else above that is noise. So, we want to select just a useful range of data to work with.
ui.notice(None, None)
ui.ignore(None, .5)
ui.ignore(5., None)
ui.plot_data(26133)
We know this is a star and stars have collisionally ionized plasma almost in ionization equilibrium. The model to describe this type of plasma is called APEC (xsvapec in Sherpa) and we give it the name "a1" (using the non-stanard sherpa way that tried to hide the fact that you really are using Python here, so the syntax is a little odd).
ui.set_source(26133, ui.xsvapec.a1 + ui.xsvapec.a2)
ui.fit(26133)
ui.plot_fit(26133)
ax = plt.gca()
ax.set_xscale("log")
ax.set_yscale("log")
ax.get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
minor_thresholds=(2, .5)))
ax.tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
ui.conf(26133)
ui.calc_energy_flux(id=26133)
ui.plot_energy_flux(id=26133)
Result¶
We fit a single temperature plasma with a temperature of 1.3 +/0.04 keV and a flux of (3.4 +/ 0.4) * e-13 erg/s/cm^2.
There are additional things we could do, e.g. treat the background, but for this data, we don't really need to (the source is bright, so the background is negligible). There are also some details to take into account for errors and statistics (the data really is Poisson distributed, I just estimated the uncertainty on the flux by hand, the fit is not terribly good below 1 keV, ...) but this is good enough for a first impression and, depending on the exact question we want to answer, might already be good enough for science analysis.
Put your result in here: https://docs.google.com/spreadsheets/d/167MmraWsYMN7fujjxG1-WpfvtXYALxTCIeLDqSLbZH0/edit?usp=sharing
Download the notebook.