{ "cells": [ { "cell_type": "markdown", "id": "7a454ff3", "metadata": {}, "source": [ "# PlotXspec example\n", "\n", "## Basic PyXspec usage & additional plotting\n", "\n", "This notebook is aimed at illustrating the _basic_ usage of PyXspec, as well as a few supporting methods related to plotting. The supporting methods are part of the [PlotXspec class available for download here](https://github.com/JohannesBuchner/BXA/tree/master/examples/xspec/bayesian-workflow). The notebook will cover how to load spectra, define a model, fit it to the data, and plot the results. PyXspec makes all Xspec functionality accessible through Python, so has a much broader range of options than covered here. For further details on PyXspec usage, please see the [documentation](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html).\n", "\n", "In the following, we will assume a basic familiarity with both Python and Xspec usage. Users will need to have a working HEASOFT installation, as well as a working Python3 installation. The examples focus on the 'out-of-the-box' usage of the PlotXSpec public methods. Some users may find it more helpful to use these methods as a template for their own plotting functions; the PlotXspec class includes several private methods that could also be of use as a guide (such as the ```_get_xspec_data()``` method).\n", "\n", "_Please note_: this notebook relates to PyXspec only. For examples of using BXA and the PlotBXA class, please see the notebook [_tutorial_usage_plotbxa_](https://johannesbuchner.github.io/BXA/tutorial_usage_plotbxa.html). \n", "\n", "#### The notebook is structured as follows:\n", "1) Basic loading and plotting data using PyXspec\n", "2) Further plotting examples\n", "3) A slightly more complicated example (multiple spectra & simultaneous fitting)\n", "\n", "---\n", "This notebook was tested with:\n", "- HEASOFT 6.33.1\n", "- Python 3.11.4" ] }, { "cell_type": "markdown", "id": "376b2522", "metadata": {}, "source": [ "## 1) Loading and plotting data using PyXspec" ] }, { "cell_type": "markdown", "id": "9a77f139", "metadata": {}, "source": [ "We will start with a simple example: loading a single spectrum into Xspec and fitting a model to it. The example data are a NICER spectrum for a bright extragalactic nuclear ignition and the associated background spectrum. The data have been rebinned to a minimum bin count of 20, which means we can use $\\chi^2$ statistics." ] }, { "cell_type": "code", "execution_count": null, "id": "09dd427b", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Loading the required modules\n", "import xspec\n", "import os,sys\n", "\n", "# load the PlotXspec class and create an instance\n", "sys.path.append(os.getcwd()+'/..')\n", "from plot_xspec import PlotXspec\n", "px = PlotXspec()" ] }, { "cell_type": "markdown", "id": "56349df3", "metadata": {}, "source": [ "### A quick aside on getting more information on class methods in Python\n", "\n", "If at any point you are unclear about what a method is supposed to do, or what arguments should be passed to it, you can make use of Python's built-in help functionality.\n", "\n", "This will list all public methods included in the class:\n", "```python\n", "help(px)\n", "```\n", "\n", "And this will show the 'docstring' for a particular method:\n", "```python\n", "help(px.plot_model_and_data)\n", "```\n", "\n", "**The public methods in PlotXspec all have docstrings: these explain in more detail what the methods are supposed to and what additional arguments can be passed to make the best use of them.**" ] }, { "cell_type": "markdown", "id": "f7b6b1b0", "metadata": {}, "source": [ "### A) Loading" ] }, { "cell_type": "code", "execution_count": null, "id": "1160818c", "metadata": {}, "outputs": [], "source": [ "# Global settings for XSpec\n", "xspec.Xset.abund = 'wilm' # Wilms et al. '00'\n", "xspec.Xset.xsect = 'vern' # Verner et al. '96'\n", "xspec.Xset.cosmo = '70 0 0.73' # Flat LambdaCDM" ] }, { "cell_type": "code", "execution_count": null, "id": "5d99a2eb", "metadata": {}, "outputs": [], "source": [ "# Load the data\n", "#\n", "# In this case we will explicitly set the background, RMF, and ARF in XSPEC\n", "# If all files are linked correctly, this is not necessary as XSPEC loads\n", "# the required files automatically\n", "\n", "xspec.AllData.clear() # <= just to be sure there is nothing in the way\n", "\n", "olddir = os.getcwd()\n", "os.chdir('example_data/athena/')\n", "\n", "sp = xspec.Spectrum('example-file.fak')\n", "\n", "os.chdir(olddir)\n", "\n", "# NOTE: in our example case, the response file is linked directly\n", "# from the spetral file. It is, of course, also possible to\n", "# link these files manually. In PyXspec backround, RMF and \n", "# ARF are all atributes of the 'Spectrum' class instance.\n", "# \n", "# For example (setting background, RMF, and ARF explicitly):\n", "#\n", "#sp.background = 'path/to/bg_file.bak'\n", "#sp.response = 'path/to/response_file.rmf'\n", "#sp.response.arf = 'path/to/auxiliary_file.arf'" ] }, { "cell_type": "code", "execution_count": null, "id": "2e62796d", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# The loading above will have caused a few warning messages, but we can check whether\n", "# everything is now in place\n", "xspec.AllData.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "a5a2d034", "metadata": {}, "outputs": [], "source": [ "# Set the appropriate spectral range\n", "sp.notice('all')\n", "xspec.AllData.ignore('bad')\n", "sp.ignore('**-1.0 10.0-**')" ] }, { "cell_type": "code", "execution_count": null, "id": "23e82947", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# A final look to make sure everything is in place\n", "xspec.AllData.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "16315838", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "540525c2", "metadata": {}, "source": [ "#################################################################\n", "#### NOTE:\n", "From this point on we will use Python methods to visualise the data, as implemented in PlotXspec. However, XSpec's normal visualisation can of course also still be used instead. For example, to make use of one of XSpec's plotting windows, enter\n", "\n", "```python\n", "xspec.Plot.device = \"/xs\"\n", "```\n", "and then pass any of XSpec's plotting arguments as strings to xspec.Plot, e.g.\n", "```python\n", "xspec.Plot('ldata')\n", "```\n", "#################################################################" ] }, { "cell_type": "code", "execution_count": null, "id": "efa55c61", "metadata": {}, "outputs": [], "source": [ "# A first look at the data\n", "px.first_look(ymin=-0.05,ymax=2.5,ylog=False, #<= several keyword arguments can be used, \n", " rebinsig=5,rebinbnum=40) # to manipulate the plot" ] }, { "cell_type": "markdown", "id": "f91d7aac", "metadata": {}, "source": [ "### B) Define and fit a model" ] }, { "cell_type": "markdown", "id": "61e58b5d", "metadata": {}, "source": [ "We will now define a simple model and fit it to the data. For the model, we will use redshifted blackbody, modified by Galactic absorption.\n", "\n", "There are many ways to set the parameters in XSpec, and we show only the simplest example here. For a slightly more complex usage, please see section 3) of this notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "687e0885", "metadata": {}, "outputs": [], "source": [ "# The simplest way to define a model in PyXspec\n", "#\n", "# We define the model components and then take a look\n", "# at the parameters we will need to set\n", "\n", "xspec.AllModels.clear() # <= just to sure, remove any possible existing definitions\n", "mod = xspec.Model(\"wabs*pow+gauss\") # <= the same model as used to create the spectrum" ] }, { "cell_type": "code", "execution_count": null, "id": "5e9a3b6f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Now we set each of the parameters to an initial value.\n", "# We can access the models and the parameters by name:\n", "mod.wabs.nH = 1e-2\n", "mod.powerlaw.PhoIndex = 1.7\n", "mod.powerlaw.norm = 1e-4\n", "mod.gaussian.LineE = 6.4\n", "mod.gaussian.Sigma = 0.1\n", "mod.gaussian.norm = 2\n", "\n", "# If we want to freeze one of the parameters using \n", "# this syntax, we could specify, e.g.:\n", "#\n", "# mod.gaussian.LineE.frozen = True" ] }, { "cell_type": "code", "execution_count": null, "id": "22d87335", "metadata": {}, "outputs": [], "source": [ "# Let us have another look at the model, to make sure\n", "# all parameters are set\n", "xspec.AllModels.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "fa53b7c7", "metadata": {}, "outputs": [], "source": [ "# Now we can run the fit\n", "\n", "# some fitting meta-parameters\n", "xspec.Fit.statMethod = 'chi'\n", "xspec.Fit.nIterations = 1000\n", "# and run...\n", "xspec.Fit.perform()" ] }, { "cell_type": "code", "execution_count": null, "id": "80b1cace", "metadata": {}, "outputs": [], "source": [ "# With the fit completed, we can access the usual functionalities in XSpec\n", "xspec.AllModels.calcFlux(\"1.0,10.0\")" ] }, { "cell_type": "code", "execution_count": null, "id": "8bb4e8c9", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Should we want to look at the scaling factor from rate to flux, we could do the following\n", "print(f'scaling rate (1-10keV) = {sp.flux[0]/sp.rate[0]:.3e}')" ] }, { "cell_type": "code", "execution_count": null, "id": "824b245d", "metadata": {}, "outputs": [], "source": [ "# And we can look at the resulting fit\n", "px.plot_model_and_data(ymin=1e-4,rebinsig=5,rebinbnum=40)" ] }, { "cell_type": "code", "execution_count": null, "id": "23e1d4d4", "metadata": {}, "outputs": [], "source": [ "# If we, e.g., would like to look at the unfolded spectrum,\n", "# we can specify the 'plottype' keyword argument.\n", "px.plot_model_and_data(ymin=9e-7,rebinsig=5,rebinbnum=40,plottype='unfolded')" ] }, { "cell_type": "markdown", "id": "fa2f34a0", "metadata": {}, "source": [ "##################\n", "\n", "PlotXspec also includes two methods to print a quick overview of the fitting results to screen: print_model_results() & print_errors(). The print_errors method requires the user to specify the $\\chi^2$ level for which the errors should be calculated. Of course, if the overall $\\chi^2$ statistic for the fit is too high, no error can be calculated (as XSpec raises an exception)." ] }, { "cell_type": "code", "execution_count": null, "id": "7f60cc28", "metadata": {}, "outputs": [], "source": [ "px.print_model_results()" ] }, { "cell_type": "code", "execution_count": null, "id": "10f72ff3", "metadata": {}, "outputs": [], "source": [ "px.print_errors(1.0)" ] }, { "cell_type": "markdown", "id": "a1a40c86", "metadata": {}, "source": [ "## 2) Further Plotting Examples" ] }, { "cell_type": "markdown", "id": "082cd94c", "metadata": {}, "source": [ "In some cases, it can be usefule to explore the parameter space of a model, to get a better understanding of the distribution of $\\chi^2$, and therefore of the uncertainties. We can do this with steppar. The command for steppar is passed as a string and is the same as it would be in XSpec (we can see what the indices are for the model parameters of interest using e.g. the output of print_model_results or xspec.AllModels.show() )." ] }, { "cell_type": "code", "execution_count": null, "id": "28358624", "metadata": {}, "outputs": [], "source": [ "# when running in Jupyer notebooks, it can be helpful to reduce the log-output here (especially for steppar)\n", "# -- BXA's 'XSilence' method is helpful for this too --\n", "xchat = xspec.Xset.chatter, xspec.Xset.logChatter\n", "xspec.Xset.chatter = 0\n", "xspec.Xset.logChatter = 0" ] }, { "cell_type": "code", "execution_count": null, "id": "5b4d4706", "metadata": {}, "outputs": [], "source": [ "# run steppar and plot. Here we vary parameters 1 (nH) and 2 (Gamma)\n", "xspec.Fit.steppar('1 7.5 15.5 50 2 1.15 2.9 50')\n", "px.plot_chisq_contours()" ] }, { "cell_type": "markdown", "id": "08664418", "metadata": {}, "source": [ "If it is only a single parameter we are interested in, we should consider the 1-dimensional $\\chi^2$ distribution. By passing the best-fit value of a particular parameter (in the example below: kT), we can get an estimate on the error in this way." ] }, { "cell_type": "code", "execution_count": null, "id": "46c47230", "metadata": {}, "outputs": [], "source": [ "par = mod.powerlaw.PhoIndex #<= the parameter of interest\n", "xspec.Fit.steppar('2 1.2 2.9 1000')\n", "px.calc_error_from_1Dchisq(par,level=1)" ] }, { "cell_type": "code", "execution_count": null, "id": "03f452a7", "metadata": {}, "outputs": [], "source": [ "# set XSpec's outbut levels back to their original values\n", "xspec.Xset.chatter, xspec.Xset.logChatter = xchat" ] }, { "cell_type": "code", "execution_count": null, "id": "c339891b", "metadata": {}, "outputs": [], "source": [ "# Note that the error value is above is different from the initial estimate shown in the\n", "# fitting results overview (see output below). This is because the method used above \n", "# provided a different (more refined) approach for the calculation of the uncertainties\n", "xspec.AllModels.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "ccd298d1", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a0022a2a", "metadata": {}, "source": [ "## 3) More detailed fitting example (XMM-Newton: pn, MOS1, and MOS2)" ] }, { "cell_type": "markdown", "id": "35bcb75e", "metadata": {}, "source": [ "When fitting multiple datasets simultaneously, we need to slightly adjust our fitting procedure. In the following we will load XMM-Newton data from all three EPIC instruments. To account for a possible offset in the flux level between the spectra, we will multiply our model with a constant factor, which will be allowed to vary among the spectra. All other model parameters will remain tied.\n", "\n", "The example data used here are for an AGN showing intrinsic obscuration and a weak Iron K$\\alpha$ line. We will fit these data with a redshifted power-law model, modified with both intrinsic and Galactic absorption." ] }, { "cell_type": "code", "execution_count": null, "id": "14a5a4e2", "metadata": {}, "outputs": [], "source": [ "# Load the data\n", "xspec.AllData.clear()\n", "xspec.AllModels.clear()\n", "\n", "olddir = os.getcwd()\n", "os.chdir('example_data/xmm/')\n", "\n", "epicfn = 'epic_pn_agn.fak'\n", "mos1fn = 'epic_mos1_agn.fak'\n", "mos2fn = 'epic_mos2_agn.fak'\n", "\n", "# load into different data groups\n", "xspec.AllData(f\"1:1 {epicfn} 2:2 {mos1fn} 3:3 {mos2fn}\")\n", "\n", "os.chdir(olddir)\n", "\n", "# Spectral range\n", "xspec.AllData.notice('all')\n", "xspec.AllData.ignore('bad')\n", "for ii in range(1,4):\n", " xspec.AllData(ii).ignore('**-0.2 10.0-**')" ] }, { "cell_type": "code", "execution_count": null, "id": "f712a553", "metadata": {}, "outputs": [], "source": [ "# A first look at the data\n", "px.first_look(ymin=1e-4,ymax=1e-1,rebinsig=5,rebinbnum=20)" ] }, { "cell_type": "code", "execution_count": null, "id": "3a2763bb", "metadata": {}, "outputs": [], "source": [ "# Set all necessary parameters in XSpec & define the model\n", "#\n", "# In the simulated data, there is an offset in flux between the pn and the\n", "# MOS1 & MOS2 spectra. We will account for this using multiplicative\n", "# constant in the model. In the model definition, we need to make sure only\n", "# the 'constant' factor is allowed to vary for the MOS1 & MOS2 spectra, \n", "# w.r.t. the pn spectrum. All other parameters should remain tied together.\n", "\n", "## initial settings\n", "xspec.Xset.abund = 'wilm' # Wilms et al. '00'\n", "xspec.Xset.xsect = 'vern' # Verner et al. '96'\n", "xspec.Xset.cosmo = '70 0 0.73'\n", "\n", "z_opt = 0.015 # assume this is known\n", "\n", "xspec.Fit.statMethod = 'chi'\n", "\n", "## define model\n", "xspec.AllModels += (\"constant(zTBabs*(zpowerlw))\")\n", "mod = xspec.AllModels(1) # <= we explicitly set mod to the first of the models \n", " # (XSspec has now loaded three, one for each spectrum)\n", "\n", "# Let us have a first look\n", "# xspec.AllModels.show() # <= only necessary if the model is not automatically printed below" ] }, { "cell_type": "code", "execution_count": null, "id": "0e0f5c76", "metadata": {}, "outputs": [], "source": [ "# Silence XSpec for the moment...\n", "xchat = xspec.Xset.chatter, xspec.Xset.logChatter\n", "xspec.Xset.chatter = 0\n", "xspec.Xset.logChatter = 0\n", "\n", "# Set parameter definitions\n", "mod.constant.factor.values = (1,-1) #<= note the different notation here; \n", " # 1 is the initial value, -1 means frozen\n", "mod.zTBabs.Redshift.values = (z_opt, -1.)\n", "mod.zTBabs.nH.values = (1., 0.01, 1e-3, 1e-3, 10., 10.) #<= standard XSpec syntax for parameter definition\n", "\n", "mod.zpowerlw.PhoIndex.values = (1.8, 0.1, -0.5, -0.5, 5., 5.)\n", "mod.zpowerlw.Redshift.values = (z_opt, -1)\n", "mod.zpowerlw.norm.values = (1.e-1, 0.01, 1.e-5, 1.e-5, 1e1, 1e1)\n", "\n", "# Below we untie the model parameters for the MOS data. As noted,\n", "# there are in fact three models (one per data group). In PyXspec\n", "# we can access each model by its index in the AllModels object.\n", "# Below, we first untie the contant parameter from its counterpart\n", "# in the first model and then unfreeze it, by giving it a possible\n", "# range of values.\n", "xspec.AllModels(2).constant.factor.untie()\n", "xspec.AllModels(2).constant.factor.values = (1, 0.01, 0.5, 0.5, 1.5, 1.5)\n", "xspec.AllModels(3).constant.factor.untie()\n", "xspec.AllModels(3).constant.factor.values = (1, 0.01, 0.5, 0.5, 1.5, 1.5)\n", "\n", "# ... and show XSpec output again\n", "xspec.Xset.chatter, xspec.Xset.logChatter = xchat\n", "\n", "# Let's have another look at the models now\n", "xspec.AllModels.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "c5f7c465", "metadata": {}, "outputs": [], "source": [ "# now we can run the fit\n", "xspec.Fit.perform()" ] }, { "cell_type": "code", "execution_count": null, "id": "c0d4143d", "metadata": {}, "outputs": [], "source": [ "# ... and plot the results\n", "px.plot_model_and_data(ymin=1e-4,ymax=1,rebinsig=5,rebinbnum=40)" ] }, { "cell_type": "code", "execution_count": null, "id": "1dd35c6f", "metadata": {}, "outputs": [], "source": [ "# or if we want to look at e.g. only the MOS1 spectrum\n", "px.plot_model_and_data(ymin=1e-4,ymax=1,idsp=2,rebinsig=5)" ] }, { "cell_type": "code", "execution_count": null, "id": "e118dcd5", "metadata": {}, "outputs": [], "source": [ "# A quick overview of the fit results\n", "px.print_model_results()" ] }, { "cell_type": "code", "execution_count": null, "id": "9564e70e", "metadata": {}, "outputs": [], "source": [ "# ... and of the errors\n", "px.print_errors(1.0)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }