Add show_sensitivity.py script.

Jürgen Knödlseder
1 parent cbe8c747
Showing 2 changed files with 119 additions and 1 deletions
examples/README
examples/show_sensitivity.py
@@ -69,7 +69,10 @@ show_pull_histogram.py
   Displays pull histograms generated by the cspull script.
  
 show_response.py
-  Display one or several CTA response functions.
+  Displays one or several CTA response functions.
+
+show_sensitivity.py
+  Displays sensitivity information generated by the cssens script.
  
 show_ts_distribution.py
   Displays Test Statistic histograms generated using the cstsdist
+#! /usr/bin/env python
+# ==========================================================================
+# This script displays a differential sensitivity plot generated by cssens.
+#
+# Copyright (C) 2015 Juergen Knoedlseder
+#
+# 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 sys
+import math
+import gammalib
+try:
+    import matplotlib.pyplot as plt
+except:
+    sys.exit("This script needs matplotlib")
+
+
+# ===================== #
+# Read sensitivity data #
+# ===================== #
+def read_sensitivity(filename):
+    """
+    Read sensitivity information from CSV file.
+    """
+    # Read filename
+    csv = gammalib.GCsv(filename,',')
+
+    # Create dictionnary
+    sensitivity = {}
+    for column in range(csv.ncols()):
+        name   = csv[0,column]
+        values = []
+        for row in range(csv.nrows()-1):
+            values.append(float(csv[row+1,column]))
+        sensitivity[name] = values
+
+    # Add linear energy values
+    if sensitivity.has_key("loge"):
+        name   = 'energy'
+        values = []
+        for value in sensitivity["loge"]:
+            values.append(math.pow(10.0, value))
+        sensitivity[name] = values
+
+    # Return
+    return sensitivity
+
+
+# ===================== #
+# Show sensitivity data #
+# ===================== #
+def show_sensitivity(sensitivity, title="Sensitivity"):
+    """
+    Plot sensitivity.
+    """
+    # Create figure
+    plt.figure()
+
+    # Set plot attributes
+    plt.title(title)
+    plt.loglog()
+    plt.grid()
+
+    # Show differential sensitivity
+    plt.plot(sensitivity['energy'], sensitivity['eflux'], 'ro-')
+
+    # Set labels
+    plt.xlabel("Energy (TeV)")
+    plt.ylabel(r"Sensitivity (erg cm$^{-2}$ s$^{-1}$)")
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    """
+    # Print usage information
+    usage = "Usage: show_sensitivity [filename]"
+    if len(sys.argv) < 1:
+        sys.stdout.write(usage+"\n")
+        sys.exit()
+
+    # Extract parameters
+    if len(sys.argv) == 1:
+        filename = 'sensitivity.dat'
+    else:
+        filename = sys.argv[1]
+
+    # Set title
+    title = 'Differential sensitivity ('+filename+')'
+
+    # Read sensitivity data
+    sensitivity = read_sensitivity(filename)
+
+    # Show sensitivity data
+    show_sensitivity(sensitivity, title=title)