From 2624e24576d74de015a631e48028db2fd8f72147 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 12 Jun 2018 15:56:42 +0200 Subject: [PATCH] Add Xspec howto tutorial to User documentation (#2404). --- ChangeLog | 3 ++- NEWS | 5 ++++- doc/source/admin/release_history/1.6.rst | 4 +++- doc/source/users/tutorials/howto/howto_xspec.rst | 214 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/tutorials/howto/index.rst | 1 + 5 files changed, 224 insertions(+), 3 deletions(-) create mode 100644 doc/source/users/tutorials/howto/howto_xspec.rst diff --git a/ChangeLog b/ChangeLog index e7c445c..dc838ec 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2018-06-11 +2018-06-12 * Version 1.6.0 released ======================== @@ -7,6 +7,7 @@ Remove ctools::is_valid_filename() method Fix bug when using precomputed model in csresspec (#2448) Keep existing library path in installation script on Mac OS X (#2420) + csphagen produces Xspec compliant data (#2404) Add show_rmf.py example script Allow specifying of input sky map in ctskymap (#2311) Add ras() and decs() methods to csobsinfo (#2313) diff --git a/NEWS b/NEWS index c73dd3a..6e3b5cb 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.6.0 -7 June 2018 +12 June 2018 Introduction @@ -222,6 +222,9 @@ The true number of energy bins is now computed using making sure in addition that the result is at least 1. This provides better control over the number of energy bins (#2490). +The script now produces data that are compliant with Xspec. An Xspec howto +tutorial was added to the user documentation (#2404). + csphasecrv - Generate phase curve --------------------------------- diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 1dc6e1e..1828f86 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -10,7 +10,7 @@ ctools 1.6 is a major release that adds significant functionality. In particular, this release provides: -* ... +* Support for Xspec spectral analysis Bug fixes @@ -27,6 +27,8 @@ Improvements * [`2490 `_] - Change formula for number of true energy bins in :ref:`csphagen` +* [`2404 `_] - + :ref:`csphagen` produces Xspec compliant data * Add ``show_rmf.py`` example script * [`2313 `_] - Add `ras()` and `decs()` methods to :ref:`csobsinfo` diff --git a/doc/source/users/tutorials/howto/howto_xspec.rst b/doc/source/users/tutorials/howto/howto_xspec.rst new file mode 100644 index 0000000..001dc9d --- /dev/null +++ b/doc/source/users/tutorials/howto/howto_xspec.rst @@ -0,0 +1,214 @@ +.. _howto_xspec: + +How to perform spectral fitting using Xspec? +-------------------------------------------- + + .. admonition:: What you will learn + + You will learn how to **use Xspec to perform a spectral fitting of On/Off + data**. + + .. warning:: You can reproduce this tutorial only if you have + `Xspec `_ installed. + You may + `download Xspec here `_. + +If you are familiar with X-ray data analysis, or if you want to use a specific +spectral model that is available in the +`Xspec `_ +package, you can use that package to fit CTA data that were prepared for an +On/Off analysis. + +To illustrate how you can do this, let's simulate an observation of the Crab +nebula with a pointing direction that is offset by 1 deg: + +.. code-block:: bash + + $ ctobssim + RA of pointing (degrees) (0-360) [83.63] + Dec of pointing (degrees) (-90-90) [22.01] 21.01 + Radius of FOV (degrees) (0-180) [5] + Start time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:00:00] + Stop time (UTC string, JD, MJD or MET in seconds) [2020-01-01T00:30:00] + Lower energy limit (TeV) [0.1] + Upper energy limit (TeV) [100.0] + Calibration database [prod2] + Instrument response function [South_0.5h] + Input model definition XML file [$CTOOLS/share/models/crab.xml] + Output event data file or observation definition XML file [events.fits] + +No use the :ref:`csphagen` script to generate the On/Off data using a On +region radius of 0.2 deg and the reflected region model to determine the +Off counts spectrum: + +.. code-block:: bash + + $ csphagen + Input event list or observation definition XML file [obs.xml] events.fits + Calibration database [prod2] + Instrument response function [South_0.5h] + Start value for first energy bin in TeV [0.1] + Stop value for last energy bin in TeV [100] + Radius of source region circle (deg) (0-180) [0.2] + Input model definition XML file (if NONE, use point source) [NONE] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Number of energy bins [20] + Stack multiple observations into single PHA, ARF and RMF files? [no] + Output observation definition XML file [onoff_obs.xml] + Method for background estimation (REFLECTED|CUSTOM) [REFLECTED] + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] + Right Ascension of source region centre (deg) (0-360) [83.63] + Declination of source region centre (deg) (-90-90) [22.01] + +Now you can start +`Xspec `_ +for model fitting. + +.. code-block:: bash + + $ xspec + Creating a $HOME/.xspec directory for you + + XSPEC version: 12.10.0c + Build Date/Time: Tue Jun 12 12:46:58 2018 + + XSPEC12> + +First load the On file ``onoff_pha_on.fits`` using the ``data`` command: + +.. code-block:: bash + + XSPEC12>data onoff_pha_on.fits + ***Warning: Data file onoff_pha_on.fits has both POISSERR key set to 'true' and a STAT_ERR column. + XSPEC will assume Poisson errors. + ***Warning: Background file onoff_pha_off.fits has both POISSERR key set to 'true' and a STAT_ERR column. + XSPEC will use assume Poisson errors. + + 1 spectrum in use + + Spectral Data File: onoff_pha_on.fits Spectrum 1 + Net count rate (cts/s) for Spectrum:1 1.585e+00 +/- 3.099e-02 (94.0 % total) + Assigned to Data Group 1 and Plot Group 1 + Noticed Channels: 1-20 + Telescope: unknown Instrument: unknown Channel Type: PI + Exposure Time: 1764 sec + Using fit statistic: chi + Using test statistic: chi + Using Background File onoff_pha_off.fits + Background Exposure Time: 1764 sec + Using Response (RMF) File onoff_rmf.fits for Source 1 + Using Auxiliary Response (ARF) File onoff_arf.fits + +This will load the On ``PHA`` file as well as the Off ``PHA``file, the ``ARF`` +file and the ``RMF`` file. While you specified the On file on the command line, +the file names of the other files were extracted from ``backfile``, ``ancrfile`` +and ``respfile`` header keywords in the On file that were setup by the +:ref:`csphagen` script. + +Now set the fit statistic to ``cstat`` using + +.. code-block:: bash + + XSPEC12>statistic cstat + Default fit statistic is set to: C-Statistic + This will apply to all current and newly loaded spectra. + +and set the spectral model as well as the initial parameters of the spectral +model using + +.. code-block:: bash + + XSPEC12>model powerlaw + + Input parameter value, delta, min, bot, top, and max values for ... + 1 0.01( 0.01) -3 -2 9 10 + 1:powerlaw:PhoIndex>2.5 + 1 0.01( 0.01) 0 0 1e+20 1e+24 + 2:powerlaw:norm>500.0 + + ======================================================================== + Model powerlaw<1> Source No.: 1 Active/On + Model Model Component Parameter Unit Value + par comp + 1 1 powerlaw PhoIndex 2.50000 +/- 0.0 + 2 1 powerlaw norm 500.000 +/- 0.0 + ________________________________________________________________________ + + + Fit statistic : C-Statistic = 712.82 using 20 PHA bins and 18 degrees of freedom. + + Test statistic : Chi-Squared = 605.26 using 20 PHA bins. + Reduced chi-squared = 33.625 for 18 degrees of freedom + Null hypothesis probability = 6.662146e-117 + + ***Warning: Chi-square may not be valid due to bins with zero variance + in spectrum number(s): 1 + + Current data and model not fit yet. + +Note that +`Xspec `_ +energies are in keV, and the power law normalization is taken at a reference +energy of 1 keV. + +Finally, you are ready to do the spectral fitting using the ``fit`` +command: + +.. code-block:: bash + + XSPEC12>fit + Parameters + C-Statistic |beta|/N Lvl 1:PhoIndex 2:norm + 85.7044 11239.1 -3 2.49284 654.424 + 25.1265 3951.81 -4 2.48007 578.148 + 23.0302 769.247 -5 2.47567 541.654 + 22.9323 165.457 -6 2.47509 538.451 + 22.9315 13.9094 -7 2.47503 538.072 + ============================== + Variances and Principal Axes + 1 2 + 8.8843E-07| -1.0000 0.0001 + 2.8990E+04| -0.0001 -1.0000 + ------------------------------ + + ======================== + Covariance Matrix + 1 2 + 2.326e-04 2.592e+00 + 2.592e+00 2.899e+04 + ------------------------ + + ======================================================================== + Model powerlaw<1> Source No.: 1 Active/On + Model Model Component Parameter Unit Value + par comp + 1 1 powerlaw PhoIndex 2.47503 +/- 1.52504E-02 + 2 1 powerlaw norm 538.072 +/- 170.265 + ________________________________________________________________________ + + + Fit statistic : C-Statistic = 22.93 using 20 PHA bins and 18 degrees of freedom. + + Test statistic : Chi-Squared = 437.54 using 20 PHA bins. + Reduced chi-squared = 24.308 for 18 degrees of freedom + Null hypothesis probability = 1.317221e-81 + + ***Warning: Chi-square may not be valid due to bins with zero variance + in spectrum number(s): 1 + + XSPEC12> + +The fit results are compared in the table below to the simulated true values +and the values obtained using :ref:`ctlike` using the same spectral model. +The +`Xspec `_ +results are compatible with the simulated values and very close to the values +obtained using :ref:`ctlike`. + + +-----------+-------+-----------------+-----------------+ + | Parameter | Truth | Xspec | ctlike | + +===========+=======+=================+=================+ + | Prefactor | 601.4 | 538.0 +/- 170.3 | 536.4 +/- 184.2 | + +-----------+-------+-----------------+-----------------+ + | Index | 2.48 | 2.475 +/- 0.015 | 2.475 +/- 0.017 | + +-----------+-------+-----------------+-----------------+ diff --git a/doc/source/users/tutorials/howto/index.rst b/doc/source/users/tutorials/howto/index.rst index fa0780d..a8b22ab 100755 --- a/doc/source/users/tutorials/howto/index.rst +++ b/doc/source/users/tutorials/howto/index.rst @@ -17,6 +17,7 @@ CTA data analysis howto_advanced_python_fitting howto_connect_irfs howto_connect_models + howto_xspec CTA first Data Challenge data analysis -- libgit2 0.21.4