Veronika Vodeb / ctools

Commit 3836359d46e2f3e52fd84a13615c6ec4e3ad71c6

Authored by L. Tibaldo
1 parent bc4de7c6

example script

Inline Side-by-side
Showing 1 changed file with 288 additions and 0 deletions
examples/show_residuals.py 0 → 100755
  View file @3836359
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Display residuals generated by csresspec
  4 +#
  5 +# Copyright (C) 2017- Luigi Tibaldo
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import sys
  22 +import math
  23 +import gammalib
  24 +import cscripts
  25 +
  26 +try:
  27 + import matplotlib.pyplot as plt
  28 + from matplotlib import gridspec
  29 +
  30 + plt.figure()
  31 + plt.close()
  32 +except (ImportError, RuntimeError):
  33 + print('This script needs the "matplotlib" module')
  34 + sys.exit()
  35 +
  36 +
  37 +# ================================= #
  38 +# Fill counts/model/residual arrays #
  39 +# ================================= #
  40 +def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
  41 + c_resid, algorithm):
  42 + """
  43 + Helper function that fills in the counts, model, residuals lists and
  44 + calculate relative errors for given table row
  45 + :param row: int row number
  46 + :param counts: list
  47 + :param model: list
  48 + :param resid: list
  49 + :param e_counts: list
  50 + :param e_resid: list
  51 + :param c_counts: `~gammalib.GFits.table.column'
  52 + :param c_model: `~gammalib.GFits.table.column'
  53 + :param c_resid: `~gammalib.GFits.table.column'
  54 + :param algorithm: string
  55 + :return: counts, model, resid, e_counts, e_resid
  56 + """
  57 + counts.append(c_counts.real(row))
  58 + model.append(c_model.real(row))
  59 + resid.append(c_resid.real(row))
  60 +
  61 + # calculate count error
  62 + err = math.sqrt(c_counts.real(row))
  63 + if err == 1:
  64 + # this prevents visualization problem in matplotlib
  65 + err = 0.99
  66 + else:
  67 + pass
  68 + e_counts.append(err)
  69 +
  70 + # calculate residual error
  71 + if algorithm == 'SUB':
  72 + e_resid.append(err)
  73 + elif algorithm == 'SUBDIV':
  74 + e_resid.append(err / c_model.real(row))
  75 + elif algorithm == 'SUBDIVSQRT':
  76 + e_resid.append(err / math.sqrt(c_model.real(row)))
  77 + elif algorithm == 'SIGNIFICANCE':
  78 + e_resid.append(1.)
  79 +
  80 + return counts, model, resid, e_counts, e_resid
  81 +
  82 +
  83 +# ============= #
  84 +# Plot spectrum #
  85 +# ============= #
  86 +def plot_residuals(filename, plotfile, hdu):
  87 + """
  88 + Plot spectrum
  89 +
  90 + Parameters
  91 + ----------
  92 + filename : str
  93 + Name of spectrum FITS file
  94 + plotfile : str
  95 + Plot file name
  96 + hdu : int
  97 + Number of observation to plot
  98 + """
  99 + # Read spectrum file
  100 + fits = gammalib.GFits(filename)
  101 + table = fits.table(1 + hdu)
  102 + c_emin = table['Emin']
  103 + c_emax = table['Emax']
  104 + c_counts = table['Counts']
  105 + c_model = table['Model']
  106 + c_resid = table['Residuals']
  107 + try:
  108 + c_counts_off = table['Counts_Off']
  109 + c_model_off = table['Model_Off']
  110 + c_resid_off = table['Residuals_Off']
  111 + is_onoff = True
  112 + except:
  113 + is_onoff = False
  114 +
  115 + # Initialise arrays to be filled
  116 + ebounds = []
  117 + ed_engs = []
  118 + eu_engs = []
  119 + em_engs = []
  120 + counts = []
  121 + e_counts = []
  122 + model = []
  123 + resid = []
  124 + e_resid = []
  125 + counts_off = []
  126 + e_counts_off = []
  127 + model_off = []
  128 + resid_off = []
  129 + e_resid_off = []
  130 +
  131 + # Residual algorithm
  132 + algorithm = table.header()['ALGORITHM'].string()
  133 +
  134 + # Loop over rows of the file
  135 + nrows = table.nrows()
  136 + # add first energy boundary
  137 + ebounds.append(c_emin.real(0))
  138 + for row in range(nrows):
  139 + # boundaries
  140 + ebounds.append(c_emax.real(row))
  141 + # geometrical mean energy and errors
  142 + emean = math.sqrt(c_emin.real(row) * c_emax.real(row))
  143 + em_engs.append(emean)
  144 + ed_engs.append(emean - c_emin.real(row))
  145 + eu_engs.append(c_emax.real(row) - emean)
  146 + # counts, model, residuals
  147 + counts, model, resid, e_counts, e_resid = fill_cmr(row, counts, model,
  148 + resid, e_counts,
  149 + e_resid, c_counts,
  150 + c_model, c_resid,
  151 + algorithm)
  152 + # onoff
  153 + if is_onoff:
  154 + counts_off, model_off, resid_off, e_counts_off, e_resid_off = fill_cmr(
  155 + row, counts_off,
  156 + model_off, resid_off,
  157 + e_counts_off,
  158 + e_resid_off,
  159 + c_counts_off,
  160 + c_model_off, c_resid_off,
  161 + algorithm)
  162 +
  163 + # Add model value to be compatible with plt.step
  164 + model = [model[0]] + model
  165 + if is_onoff:
  166 + model_off = [model_off[0]] + model_off
  167 +
  168 + # Plot
  169 + axarr = []
  170 + if is_onoff:
  171 + f = plt.figure(figsize=(8, 4))
  172 + gs = gridspec.GridSpec(2, 2, height_ratios=[2, 1])
  173 + gs.update(hspace=0, left=0.09, right=0.97)
  174 + axarr.append(f.add_subplot(gs[0]))
  175 + axarr.append(f.add_subplot(gs[2]))
  176 + axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
  177 + axarr.append(f.add_subplot(gs[3], sharex=axarr[1]))
  178 + else:
  179 + f = plt.figure()
  180 + gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
  181 + gs.update(hspace=0)
  182 + axarr.append(f.add_subplot(gs[0]))
  183 + axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
  184 + plt.setp(axarr[1].get_xticklabels(), visible=False)
  185 + axarr[0].set_xscale('log')
  186 + axarr[1].set_xscale('log')
  187 + axarr[0].set_yscale('log')
  188 + # Counts and model
  189 + axarr[0].errorbar(em_engs, counts, yerr=e_counts, xerr=[ed_engs, eu_engs],
  190 + fmt='ko', capsize=0, linewidth=2, zorder=2, label='Data')
  191 + axarr[0].step(ebounds, model, color='0.5', linewidth=2, zorder=1,
  192 + label='Model')
  193 + axarr[0].set_ylabel('Counts')
  194 + # Residuals
  195 + axarr[1].errorbar(em_engs, resid, yerr=e_resid, xerr=[ed_engs, eu_engs],
  196 + fmt='ko', capsize=0, linewidth=2, zorder=2)
  197 + axarr[1].axhline(0, color='0.5', linestyle='--')
  198 + axarr[1].set_xlabel('Energy (TeV)')
  199 + if algorithm == 'SUB':
  200 + axarr[1].set_ylabel('Residuals (counts)')
  201 + elif algorithm == 'SUBDIV':
  202 + axarr[1].set_ylabel('Residuals (fraction)')
  203 + elif algorithm == 'SUBDIVSQRT' or algorithm == 'SIGNIFICANCE':
  204 + axarr[1].set_ylabel(r'Residuals ($\sigma$)')
  205 + # On/Off
  206 + if is_onoff:
  207 + axarr[0].set_title('ON')
  208 + axarr[2].set_title('OFF')
  209 + plt.setp(axarr[3].get_xticklabels(), visible=False)
  210 + axarr[2].set_xscale('log')
  211 + axarr[3].set_xscale('log')
  212 + axarr[2].set_yscale('log')
  213 + # Counts and model
  214 + axarr[2].errorbar(em_engs, counts_off, yerr=e_counts_off,
  215 + xerr=[ed_engs, eu_engs],
  216 + fmt='ko', capsize=0, linewidth=2, zorder=2,
  217 + label='Data')
  218 + axarr[2].step(ebounds, model_off, color='0.5', linewidth=2, zorder=1,
  219 + label='Model')
  220 + # Residuals
  221 + axarr[3].errorbar(em_engs, resid_off, yerr=e_resid_off,
  222 + xerr=[ed_engs, eu_engs],
  223 + fmt='ko', capsize=0, linewidth=2, zorder=2)
  224 + axarr[3].axhline(0, color='0.5', linestyle='--')
  225 + axarr[3].set_xlabel('Energy (TeV)')
  226 +
  227 + # Add spectra of individual components
  228 + skiplist = ['Counts', 'Model', 'Residuals',
  229 + 'Counts_Off', 'Model_Off', 'Residuals_Off',
  230 + 'Emin', 'Emax']
  231 + for s in range(table.ncols()):
  232 + if table[s].name() in skiplist:
  233 + pass
  234 + else:
  235 + component = []
  236 + for row in range(nrows):
  237 + component.append(table[s].real(row))
  238 + component = [component[0]] + component
  239 + axarr[0].step(ebounds, component, zorder=0,
  240 + label=table[s].name())
  241 +
  242 + # Add legend
  243 + axarr[0].legend(loc='best')
  244 +
  245 + # Show figure
  246 + if len(plotfile) > 0:
  247 + plt.savefig(plotfile)
  248 + else:
  249 + plt.show()
  250 +
  251 + # Return
  252 + return
  253 +
  254 +
  255 +# =============== #
  256 +# Show residuals #
  257 +# =============== #
  258 +def show_residuals():
  259 + """
  260 + Show residuals
  261 + """
  262 + # Set usage string
  263 + usage = 'show_spectrum.py [-p plotfile] [-h hdu] [file]'
  264 +
  265 + # Set default options
  266 + options = [{'option': '-p', 'value': ''},
  267 + {'option': '-h', 'value': 0}]
  268 +
  269 + # Get arguments and options from command line arguments
  270 + args, options = cscripts.ioutils.get_args_options(options, usage)
  271 +
  272 + # Extract script parameters from options
  273 + plotfile = options[0]['value']
  274 + hdu = options[1]['value']
  275 +
  276 + # Show residuals
  277 + plot_residuals(args[0], plotfile, hdu)
  278 +
  279 + # Return
  280 + return
  281 +
  282 +
  283 +# ======================== #
  284 +# Main routine entry point #
  285 +# ======================== #
  286 +if __name__ == '__main__':
  287 + # Show spectrum
  288 + show_residuals()
... ...