Mathieu Servillat / ctools

Commit a6dd93fdac084eed879b9803b50552ff6ec09da1

Authored by Jürgen Knödlseder
1 parent 7a8e2181

Add PHA spectrum display script.

Inline Side-by-side
Showing 1 changed file with 110 additions and 0 deletions
examples/show_pha.py 0 → 100755
  View file @a6dd93f
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script displays a spectrum in PHA format.
  4 +#
  5 +# Required 3rd party modules:
  6 +# - matplotlib
  7 +# - numpy
  8 +#
  9 +# Copyright (C) 2013 Juergen Knoedlseder
  10 +#
  11 +# This program is free software: you can redistribute it and/or modify
  12 +# it under the terms of the GNU General Public License as published by
  13 +# the Free Software Foundation, either version 3 of the License, or
  14 +# (at your option) any later version.
  15 +#
  16 +# This program is distributed in the hope that it will be useful,
  17 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  19 +# GNU General Public License for more details.
  20 +#
  21 +# You should have received a copy of the GNU General Public License
  22 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  23 +#
  24 +# ==========================================================================
  25 +import sys
  26 +import math
  27 +from gammalib import *
  28 +
  29 +
  30 +# ================= #
  31 +# Show PHA spectrum #
  32 +# ================= #
  33 +def show_pha(pha):
  34 + """
  35 + """
  36 + # Create figure
  37 + plt.figure(1)
  38 + plt.title("PHA spectrum (" + pha.filename() + ")")
  39 +
  40 + # Generate energy vector
  41 + energy = []
  42 + channels = False
  43 + if pha.ebounds().size() > 0:
  44 + for i in range(pha.ebounds().size()):
  45 + energy.append(pha.ebounds().elogmean(i).TeV())
  46 + else:
  47 + for i in range(pha.size()):
  48 + energy.append(float(i+1))
  49 + channels = True
  50 +
  51 + # Generate data and error vector
  52 + counts = [0.0 for i in range(pha.size())]
  53 + for i in range(pha.size()):
  54 + counts[i] = pha[i]
  55 + error = [math.sqrt(c) for c in counts]
  56 +
  57 + # Plot data
  58 + if channels:
  59 + plt.semilogy(energy, counts, 'ro')
  60 + else:
  61 + plt.loglog(energy, counts, 'ro')
  62 + plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  63 +
  64 + # Set axes
  65 + if channels:
  66 + plt.xlabel("Channels")
  67 + else:
  68 + plt.xlabel("Energy (TeV)")
  69 + plt.ylabel("Counts")
  70 +
  71 + # Show plot
  72 + plt.show()
  73 +
  74 + # Return
  75 + #return
  76 +
  77 +
  78 +# ======================== #
  79 +# Main routine entry point #
  80 +# ======================== #
  81 +if __name__ == '__main__':
  82 + """
  83 + Display spectrum in PHA format.
  84 + """
  85 + # Print usage information
  86 + usage = "Usage: show_pha filename" \
  87 + " [-n bins] [-c column] [-t title] [-p plot]"
  88 + if len(sys.argv) < 2:
  89 + print(usage)
  90 + sys.exit()
  91 +
  92 + # Extract parameters
  93 + filename = sys.argv[1]
  94 +
  95 + # Load PHA spectrum
  96 + pha = GPha(filename)
  97 + #print(pha)
  98 +
  99 + # Try importing matplotlib
  100 + try:
  101 + import matplotlib.pyplot as plt
  102 + has_matplotlib = True
  103 + except ImportError:
  104 + print("Matplotlib is not (correctly) installed on your system.")
  105 + has_matplotlib = False
  106 +
  107 + # Show PHA spectrum
  108 + if has_matplotlib:
  109 + show_pha(pha)
  110 +
0 111 \ No newline at end of file
... ...