Olga Sergijenko / ctools

Commit 1dbfd762872213946660deda9e39be24a3426389

Authored by Michael Mayer
1 parent 3cdb23e6

Added example script to visualise observation container

Inline Side-by-side
Showing 1 changed file with 137 additions and 0 deletions
examples/show_obs.py 0 → 100644
  View file @1dbfd76
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script displays a summary of an observation container stored in an XML file
  4 +#
  5 +# Required 3rd party modules:
  6 +# - matplotlib
  7 +#
  8 +# Copyright (C) 2016 Michael Mayer
  9 +#
  10 +# This program is free software: you can redistribute it and/or modify
  11 +# it under the terms of the GNU General Public License as published by
  12 +# the Free Software Foundation, either version 3 of the License, or
  13 +# (at your option) any later version.
  14 +#
  15 +# This program is distributed in the hope that it will be useful,
  16 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  18 +# GNU General Public License for more details.
  19 +#
  20 +# You should have received a copy of the GNU General Public License
  21 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  22 +#
  23 +# ==========================================================================
  24 +
  25 +import gammalib
  26 +import ctools
  27 +import cscripts
  28 +import sys
  29 +try:
  30 + import matplotlib.pyplot as plt
  31 +except:
  32 + sys.exit("This script requires matplotlib")
  33 +
  34 +
  35 +# ======================== #
  36 +# Main routine entry point #
  37 +# ======================== #
  38 +if __name__ == "__main__":
  39 + """
  40 + Print and display details of the observation container using matplotlib.
  41 + """
  42 +
  43 + # Print help if wrong number of arguments provided
  44 + if not len(sys.argv) > 1:
  45 + msg = "usage: inspect_obs.py <inobs.xml> <ra> <dec>\n"
  46 + msg += "or\n"
  47 + msg += "inspect_obs.py <inobs.xml>"
  48 + sys.exit(msg)
  49 +
  50 + # Get input observation XML file
  51 + inobs = gammalib.GFilename(sys.argv[1])
  52 +
  53 + # Assert more than one observation (circumvents having a binned observation without e.g.
  54 + # zenith and azimuth
  55 + if inobs.is_fits():
  56 + raise RuntimeError("Input observation container should be an XML file")
  57 +
  58 + # Run csobsinfo
  59 + info = cscripts.csobsinfo()
  60 + info["offset"] = False
  61 +
  62 + # Add offset if coordinates are provided
  63 + if len(sys.argv) == 4:
  64 + ra = float(sys.argv[2])
  65 + dec = float(sys.argv[3])
  66 + info["ra"] = ra
  67 + info["dec"] = dec
  68 + info["offset"] = True
  69 + info["inobs"] = inobs.url()
  70 + info["debug"] = True
  71 + info.run()
  72 +
  73 + # Retrieve observation info
  74 + zeniths = info.zeniths()
  75 + azimuths = info.azimuths()
  76 + offsets = info.offsets()
  77 + ebounds = info.ebounds()
  78 + gti = info.gti()
  79 +
  80 + # Plot zenith angle distribution
  81 + plt.figure()
  82 + zmin = min(zeniths)
  83 + zmax = max(zeniths)
  84 + plt.hist(zeniths, bins=30, range=(zmin, zmax), fc='blue')
  85 + plt.xlabel("Zenith Angle [deg]")
  86 + plt.ylabel("Abundance")
  87 + plt.title("Zenith angle distribution")
  88 +
  89 + # Plot azimuth angle distribution
  90 + plt.figure()
  91 + amin = min(azimuths)
  92 + amax = max(azimuths)
  93 + plt.hist(azimuths, bins=30, range=(amin, amax),fc='blue')
  94 + plt.xlabel("Azimuth Angle [deg]")
  95 + plt.ylabel("Abundance")
  96 + plt.title("Azimuth distribution")
  97 +
  98 + # Plot offset if possible
  99 + if info["offset"].boolean():
  100 + plt.figure()
  101 + omin = min(offsets)
  102 + omax = max(offsets)
  103 + plt.hist(offsets, bins=30, range=(omin, omax),fc='blue')
  104 + plt.xlabel("Offset from (RA, DEC)=("+str(ra)+","+str(dec)+") [deg]")
  105 + plt.ylabel("Abundance")
  106 + plt.title("Offset distribution")
  107 +
  108 + # Plot energy thresholds if possible
  109 + if ebounds.size():
  110 +
  111 + emin = []
  112 + emax = []
  113 + for i in range(ebounds.size()):
  114 + emin.append(ebounds.emin(i).log10TeV())
  115 + emax.append(ebounds.emax(i).log10TeV())
  116 + plt.figure()
  117 + plt.hist(emin, bins=80, range=(-1.0, 2.0), fc='red',label="emin")
  118 + plt.hist(emax, bins=80, range=(-1.0, 2.0), fc='blue', label ="emax")
  119 + plt.xlabel("Energy threshold [log10 (E/TeV)]")
  120 + plt.ylabel("Abundance")
  121 + plt.legend(loc="upper left")
  122 + plt.title("Energy threshold")
  123 +
  124 + # Plot observation point in time wrt zenith angle
  125 + plt.figure()
  126 + times = []
  127 + for i in range(gti.size()):
  128 + tmean = gti.tstart(i) + 0.5*(gti.tstart(i)-gti.tstart(i))
  129 + times.append(tmean.mjd())
  130 + plt.plot(times, zeniths, 'o',lw=2.0, color='black')
  131 + plt.xlabel("Time [MJD]")
  132 + plt.ylabel("Zenith Angle [deg]")
  133 + plt.title("Observation time")
  134 +
  135 + # Show plots
  136 + plt.show()
  137 +
... ...