From 1dbfd762872213946660deda9e39be24a3426389 Mon Sep 17 00:00:00 2001
From: Michael Mayer <michael.mayer@desy.de>
Date: Thu, 23 Jun 2016 17:03:47 +0200
Subject: [PATCH] Added example script to visualise observation container

---
 examples/show_obs.py | 137 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 137 insertions(+), 0 deletions(-)
 create mode 100644 examples/show_obs.py

diff --git a/examples/show_obs.py b/examples/show_obs.py
new file mode 100644
index 0000000..6dc8ee0
--- /dev/null
+++ b/examples/show_obs.py
@@ -0,0 +1,137 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This script displays a summary of an observation container stored in an XML file
+#
+# Required 3rd party modules:
+# - matplotlib
+#
+# Copyright (C) 2016 Michael Mayer
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+
+import gammalib
+import ctools
+import cscripts
+import sys
+try:
+    import matplotlib.pyplot as plt
+except:
+    sys.exit("This script requires matplotlib")
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == "__main__":
+    """
+    Print and display details of the observation container using matplotlib.
+    """
+
+    # Print help if wrong number of arguments provided
+    if not len(sys.argv) > 1:
+        msg = "usage: inspect_obs.py <inobs.xml> <ra> <dec>\n"
+        msg += "or\n"
+        msg += "inspect_obs.py <inobs.xml>"
+        sys.exit(msg)
+    
+    # Get input observation XML file
+    inobs = gammalib.GFilename(sys.argv[1])
+    
+    # Assert more than one observation (circumvents having a binned observation without e.g.
+    # zenith and azimuth
+    if inobs.is_fits():
+        raise RuntimeError("Input observation container should be an XML file")    
+    
+    # Run csobsinfo
+    info = cscripts.csobsinfo()
+    info["offset"] = False
+    
+    # Add offset if coordinates are provided
+    if len(sys.argv) == 4:
+        ra = float(sys.argv[2])
+        dec = float(sys.argv[3])
+        info["ra"] = ra
+        info["dec"] = dec
+        info["offset"] = True 
+    info["inobs"] = inobs.url()
+    info["debug"] = True
+    info.run()
+    
+    # Retrieve observation info
+    zeniths  = info.zeniths()
+    azimuths = info.azimuths()
+    offsets  = info.offsets()
+    ebounds  = info.ebounds()
+    gti      = info.gti()
+    
+    # Plot zenith angle distribution
+    plt.figure()
+    zmin = min(zeniths)
+    zmax = max(zeniths)
+    plt.hist(zeniths, bins=30, range=(zmin, zmax), fc='blue')
+    plt.xlabel("Zenith Angle [deg]")
+    plt.ylabel("Abundance")
+    plt.title("Zenith angle distribution")
+    
+    # Plot azimuth angle distribution
+    plt.figure()
+    amin = min(azimuths)
+    amax = max(azimuths)
+    plt.hist(azimuths, bins=30, range=(amin, amax),fc='blue')
+    plt.xlabel("Azimuth Angle [deg]")
+    plt.ylabel("Abundance")
+    plt.title("Azimuth distribution")
+    
+    # Plot offset if possible
+    if info["offset"].boolean():
+        plt.figure()
+        omin = min(offsets)
+        omax = max(offsets)  
+        plt.hist(offsets, bins=30, range=(omin, omax),fc='blue')
+        plt.xlabel("Offset from (RA, DEC)=("+str(ra)+","+str(dec)+") [deg]")
+        plt.ylabel("Abundance")
+        plt.title("Offset distribution")
+    
+    # Plot energy thresholds if possible
+    if ebounds.size():
+        
+        emin = []
+        emax = []
+        for i in range(ebounds.size()):
+            emin.append(ebounds.emin(i).log10TeV())
+            emax.append(ebounds.emax(i).log10TeV())
+        plt.figure()
+        plt.hist(emin, bins=80, range=(-1.0, 2.0), fc='red',label="emin")
+        plt.hist(emax, bins=80, range=(-1.0, 2.0), fc='blue', label ="emax")  
+        plt.xlabel("Energy threshold [log10 (E/TeV)]")
+        plt.ylabel("Abundance")
+        plt.legend(loc="upper left")
+        plt.title("Energy threshold")
+    
+    # Plot observation point in time wrt zenith angle
+    plt.figure()
+    times = []
+    for i in range(gti.size()):
+        tmean = gti.tstart(i) + 0.5*(gti.tstart(i)-gti.tstart(i))
+        times.append(tmean.mjd())
+    plt.plot(times, zeniths, 'o',lw=2.0, color='black')
+    plt.xlabel("Time [MJD]")
+    plt.ylabel("Zenith Angle [deg]")
+    plt.title("Observation time")
+
+    # Show plots
+    plt.show()
+    
--
libgit2 0.21.4