Michael Mayer / ctools

Commit aabec45be400d8c0f9cfad1693b772779795c07b

Authored by Jürgen Knödlseder
1 parent 628c46f3

Add show_spectrum.py script

Inline Side-by-side
Showing 1 changed file with 101 additions and 0 deletions
examples/show_spectrum.py 0 → 100755
  View file @aabec45
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script shows a spectrum created with csspec
  4 +#
  5 +# Copyright (C) 2015 Juergen Knoedlseder
  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 +
  22 +import gammalib
  23 +import sys
  24 +try:
  25 + import matplotlib.pyplot as plt
  26 +except:
  27 + sys.exit("This script needs matplotlib")
  28 +
  29 +
  30 +# ============ #
  31 +# Script entry #
  32 +# ============ #
  33 +if __name__ == "__main__":
  34 + """
  35 + """
  36 +
  37 + # Check for given spectrum file
  38 + if not len(sys.argv) == 2:
  39 + sys.exit("Usage: show_spectrum.py spectrum.fits")
  40 +
  41 + # Read spectrum file
  42 + filename = sys.argv[1]
  43 + fits = gammalib.GFits(filename)
  44 + table = fits.table(1)
  45 + c_energy = table["Energy"]
  46 + c_ed = table["ed_Energy"]
  47 + c_eu = table["eu_Energy"]
  48 + c_flux = table["Flux"]
  49 + c_eflux = table["e_Flux"]
  50 + c_ts = table["TS"]
  51 + c_upper = table["UpperLimit"]
  52 +
  53 + # Initialise arrays to be filled
  54 + energies = []
  55 + flux = []
  56 + ed_engs = []
  57 + eu_engs = []
  58 + e_flux = []
  59 + ul_energies = []
  60 + ul_ed_engs = []
  61 + ul_eu_engs = []
  62 + ul_flux = []
  63 +
  64 + # Loop over rows of the file
  65 + nrows = table.nrows()
  66 + for row in range(nrows):
  67 +
  68 + # Get TS
  69 + ts = c_ts.real(row)
  70 + flx = c_flux.real(row)
  71 + e_flx = c_eflux.real(row)
  72 +
  73 + # Switch
  74 + if ts > 9.0 and e_flx < flx:
  75 +
  76 + # Add information
  77 + energies.append(c_energy.real(row))
  78 + flux.append(c_flux.real(row))
  79 + ed_engs.append(c_ed.real(row))
  80 + eu_engs.append(c_eu.real(row))
  81 + e_flux.append(c_eflux.real(row))
  82 +
  83 + #
  84 + else:
  85 +
  86 + # Add information
  87 + ul_energies.append(c_energy.real(row))
  88 + ul_flux.append(c_upper.real(row))
  89 + ul_ed_engs.append(c_ed.real(row))
  90 + ul_eu_engs.append(c_eu.real(row))
  91 +
  92 + # Plot the spectrum
  93 + plt.figure()
  94 + plt.loglog()
  95 + plt.grid()
  96 + plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], fmt='ro')
  97 + plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], yerr=1.0e-11, uplims=True, fmt='ro')
  98 + plt.xlabel("Energy [TeV]")
  99 + plt.ylabel(r"dN/dE [erg cm$^{-2}$ s$^{-1}$]")
  100 + plt.show()
  101 +
0 102 \ No newline at end of file
... ...