diff --git a/doc/source/users/tutorials/hess_dr1/classical.rst b/doc/source/users/tutorials/hess_dr1/classical.rst new file mode 100644 index 0000000..cbdeff2 --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/classical.rst @@ -0,0 +1,130 @@ +.. _hess_dr1_classical: + +Performing a classical On/Off analysis +-------------------------------------- + + .. admonition:: What you will learn + + You will learn **how to perform a classical On/Off analysis of the + data**. + + +.. code-block:: bash + + $ csphagen use_irf_bkg=no + Input event list or observation definition XML file [obs.xml] obs_crab_selected.xml + Input model definition XML file (if NONE, use point source) [NONE] crab_models.xml + Source name [Crab] + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] + Start value for first energy bin in TeV [0.1] 0.66 + Stop value for last energy bin in TeV [100.0] + Number of energy bins [120] 40 + Stack multiple observations into single PHA, ARF and RMF files? [no] + Output observation definition XML file [onoff_obs.xml] onoff_obs_crab.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] + Radius of source region circle (deg) (0-180) [0.2] + +.. code-block:: xml + + + + + + + + + + + + + + + + +.. code-block:: bash + + $ ctlike statistic=wstat + Input event list, counts cube or observation definition XML file [events.fits] onoff_obs_crab.xml + Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_models_onoff.xml + Output model definition XML file [crab_results.xml] + +.. code-block:: none + + 2018-09-13T11:49:06: +=========================================+ + 2018-09-13T11:49:06: | Maximum likelihood optimisation results | + 2018-09-13T11:49:06: +=========================================+ + 2018-09-13T11:49:06: === GOptimizerLM === + 2018-09-13T11:49:06: Optimized function value ..: 60.048 + 2018-09-13T11:49:06: Absolute precision ........: 0.005 + 2018-09-13T11:49:06: Acceptable value decrease .: 2 + 2018-09-13T11:49:06: Optimization status .......: converged + 2018-09-13T11:49:06: Number of parameters ......: 6 + 2018-09-13T11:49:06: Number of free parameters .: 2 + 2018-09-13T11:49:06: Number of iterations ......: 5 + 2018-09-13T11:49:06: Lambda ....................: 1e-08 + 2018-09-13T11:49:06: Maximum log likelihood ....: -60.048 + 2018-09-13T11:49:06: Observed events (Nobs) ...: 605.000 + 2018-09-13T11:49:06: Predicted events (Npred) ..: 604.191 (Nobs - Npred = 0.808986318662733) + 2018-09-13T11:49:06: === GModels === + 2018-09-13T11:49:06: Number of models ..........: 1 + 2018-09-13T11:49:06: Number of parameters ......: 6 + 2018-09-13T11:49:06: === GModelSky === + 2018-09-13T11:49:06: Name ......................: Crab + 2018-09-13T11:49:06: Instruments ...............: all + 2018-09-13T11:49:06: Instrument scale factors ..: unity + 2018-09-13T11:49:06: Observation identifiers ...: all + 2018-09-13T11:49:06: Model type ................: PointSource + 2018-09-13T11:49:06: Model components ..........: "PointSource" * "PowerLaw" * "Constant" + 2018-09-13T11:49:06: Number of parameters ......: 6 + 2018-09-13T11:49:06: Number of spatial par's ...: 2 + 2018-09-13T11:49:06: RA .......................: 83.633 deg (fixed,scale=1) + 2018-09-13T11:49:06: DEC ......................: 22.015 deg (fixed,scale=1) + 2018-09-13T11:49:06: Number of spectral par's ..: 3 + 2018-09-13T11:49:06: Prefactor ................: 4.77687464316323e-17 +/- 3.12433509735449e-18 [0,infty[ ph/cm2/s/MeV (free,scale=5.7e-18,gradient) + 2018-09-13T11:49:06: Index ....................: -2.60709757251519 +/- 0.0767624316833204 [9.92,-9.92] (free,scale=-2.48,gradient) + 2018-09-13T11:49:06: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient) + 2018-09-13T11:49:06: Number of temporal par's ..: 1 + 2018-09-13T11:49:06: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) + +The intensity at 1 TeV of the Crab was fitted to +:math:`(4.78 \pm 0.31) \times 10^{-11}\,{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm TeV}^{-1}` +the spectral index of the power law is +:math:`-2.61 \pm 0.08`. +This can be compared to the values +:math:`(4.64 \pm 0.26) \times 10^{-11}\,{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm TeV}^{-1}` +and +:math:`-2.61 \pm 0.07` +that were obtained before with the unbinned maximum likelihood fit using the +energy dispersion. + +.. code-block:: bash + + $ csresspec components=yes + Input event list, counts cube, or observation definition XML file [events.fits] onoff_obs_crab.xml + Stack observations? [no] yes + Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_results_onoff.xml + Residuals computation algorithm (SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE) [SIGNIFICANCE] + Output residual spectrum file [resspec.fits] resspec_onoff.fits + +The script will produce the FITS file ``resspec.fits`` that contains the +spectral residuals for each of the model components. +To display the residuals you can use the ``show_residuals.py`` script in the +example folder as follows: + +.. code-block:: bash + + $CTOOLS/share/examples/python/show_residuals.py resspec.fits + +The figure below shows the resulting residual spectrum. +Although the model provides a reasonable first order approximation of the +spectral distribution of the data, there are still significant spectral +residuals that cannot be explained by the model. + +.. figure:: classical_resspec.png + :width: 600px + :align: center + + *Residual counts spectrum after model fitting* diff --git a/doc/source/users/tutorials/hess_dr1/classical_resspec.png b/doc/source/users/tutorials/hess_dr1/classical_resspec.png new file mode 100644 index 0000000..39a37eb --- /dev/null +++ b/doc/source/users/tutorials/hess_dr1/classical_resspec.png diff --git a/doc/source/users/tutorials/hess_dr1/index.rst b/doc/source/users/tutorials/hess_dr1/index.rst index 740515e..af0a9ea 100644 --- a/doc/source/users/tutorials/hess_dr1/index.rst +++ b/doc/source/users/tutorials/hess_dr1/index.rst @@ -10,3 +10,4 @@ Analysing public H.E.S.S. data select_events skymap fitting + classical diff --git a/examples/show_resmap.py b/examples/show_resmap.py new file mode 100755 index 0000000..fdd48e7 --- /dev/null +++ b/examples/show_resmap.py @@ -0,0 +1,119 @@ +#! /usr/bin/env python +# ========================================================================== +# Shows histogram of residual map +# +# Copyright (C) 2018 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import sys +import gammalib +import cscripts +try: + import matplotlib.pyplot as plt + import matplotlib.mlab as mlab + plt.figure() + plt.close() +except (ImportError, RuntimeError): + print('This script needs the "matplotlib" module') + sys.exit() +try: + import numpy as np +except ImportError: + print('This script needs the "numpy" module') + sys.exit() + + +# ================= # +# Plot residual map # +# ================= # +def plot_resmap(filename, nbins, plotfile): + """ + Plot pull histogram + + Parameters + ---------- + filename : str + Residual map filename + nbins : int + Number of hisogram bins + plotfile : str + Plot filename + """ + # Open residual map + map = gammalib.GSkyMap(filename) + + # Create array of values + values = map.array() + + # Create histogram + _, bins, _ = plt.hist(values, nbins, range=[-4.0,4.0], + normed=True, facecolor='green') + + # Create expected distribution + y = mlab.normpdf(bins, 0.0, 1.0) + plt.plot(bins, y, 'r-', linewidth=2) + + # Set plot + plt.xlabel('Pull') + plt.ylabel('Arbitrary units') + #plt.title(parname) + plt.grid(True) + + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# ================= # +# Show residual map # +# ================= # +def show_resmap(): + """ + Show residual map + """ + # Set usage string + usage = 'show_resmap.py [-n bins] [-p plotfile] file' + + # Set default options + options = [{'option': '-n', 'value': '50'}, + {'option': '-p', 'value': ''}] + + # Get arguments and options from command line arguments + args, options = cscripts.ioutils.get_args_options(options, usage) + + # Extract script parameters from options + nbins = int(options[0]['value']) + plotfile = options[1]['value'] + + # Plot residual map + plot_resmap(args[0], nbins, plotfile) + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Show residual map + show_resmap()