ctools

#! /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")

# ============================ #

# Overplot south sensitivities #

# ============================ #

def south_sensitivities():

    """

    """

    # Set sensitivities

    energies_50h = [0.02522, 0.03965, 0.06359, 0.10000, 0.15725, 0.24726,

                    0.39654, 0.63594, 1.00000, 1.60372, 2.47265, 3.96543,

                    6.23551, 10.0000, 16.0372, 25.2180, 39.6543, 62.3551,

                    98.0512, 157.247]

    sens_50h =     [-10.84000, -11.41818, -11.73455, -12.07273, -12.29091, -12.47636,

                    -12.67273, -12.82545, -12.97818, -13.09818, -13.24000, -13.29455,

                    -13.29455, -13.21818, -13.07636, -12.91273, -12.70545, -12.52000,

                    -12.32364, -12.07273]

    energies_5h  = [0.02522, 0.04044, 0.06359, 0.10000, 0.15725, 0.25218,

                    0.39654, 0.63594, 1.00000, 1.60372, 2.47265, 3.96543,

                    6.35944, 10.1988, 16.0372, 24.7265, 39.6543, 63.5944,

                    100.000, 160.372]

    sens_5h  =     [-10.54545, -11.04727, -11.23273, -11.56000, -11.76727, -11.95273,

                    -12.10545, -12.23636, -12.36727, -12.45455, -12.49818, -12.47636,

                    -12.43273, -12.24727, -12.10545, -11.90909, -11.70182, -11.52727,

                    -11.30909, -11.08000]

    energies_05h = [0.02522, 0.04044, 0.06486, 0.10199, 0.15725, 0.25218,

                    0.40442, 0.63594, 1.00000, 1.60372, 2.57191, 4.04425,

                    6.35944, 10.0000, 16.0372, 25.2179, 39.6543, 63.5944,

                    100.000, 160.372]

    sens_05h =     [-10.02182, -10.53455, -10.69818, -11.01455, -11.21091, -11.35273,

                    -11.49454, -11.57091, -11.64727, -11.69091, -11.65818, -11.57091,

                    -11.44000, -11.28727, -11.10182, -10.91636, -10.70909, -10.51273,

                    -10.32727, -10.08727]

    # Convert sensitivity

    flux_50h = [math.pow(10.0,x) for x in sens_50h]

    flux_5h  = [math.pow(10.0,x) for x in sens_5h]

    flux_05h = [math.pow(10.0,x) for x in sens_05h]

    # Overplot sensitivities

    plt.loglog(energies_50h, flux_50h, 'b-.')

    plt.loglog(energies_5h,  flux_5h,  'g-.')

    plt.loglog(energies_05h, flux_05h, 'r-.')

    # Return

    return

# ============================ #

# Overplot north sensitivities #

# ============================ #

def north_sensitivities():

    """

    """

    # Set sensitivities

    energies_50h = [0.02517, 0.04032, 0.06460, 0.10148, 0.15943, 0.25542,

                    0.40125, 0.64284, 1.02990, 1.68272, 2.49223, 3.91517,

                    6.15053, 10.0492, 15.7868, 25.2921, 39.7326, 63.6558,

                    100.000, 160.210]

    sens_50h =     [-10.64364, -10.99273, -11.40727, -11.80000, -12.00727, -12.29091,

                    -12.47636, -12.62909, -12.78182, -12.91273, -12.92364, -12.92364,

                    -12.81455, -12.67273, -12.49818, -12.33455, -12.12727, -11.88727,

                    -11.59273, -11.32000]

    energies_5h =  [0.02567, 0.03954, 0.06334, 0.09951, 0.15943, 0.25045,

                    0.39344, 0.64284, 1.00987, 1.58645, 2.54166, 3.91517,

                    6.15053, 9.85379, 16.0998, 24.8002, 39.7326, 62.4179,

                    100.000, 157.094]

    sens_5h =      [-10.14182, -10.66545, -10.86182, -11.25455, -11.46182, -11.72364,

                    -11.89818, -11.98545, -12.09455, -12.16000, -12.12727, -12.06182,

                    -11.89818, -11.70182, -11.49454, -11.30909, -11.09091, -10.84000,

                    -10.54545, -10.28364]

    energies_05h = [0.02776, 0.04032, 0.06334, 0.10148, 0.15943, 0.25045,

                    0.40125, 0.64284, 0.99023, 1.58645, 2.54166, 3.99282,

                    6.27251, 10.0492, 15.7868, 25.2921, 38.9599, 62.4179,

                    77.4686]

    sens_05h =     [-9.71636, -10.15273, -10.33818, -10.68727, -10.88364, -11.13455,

                    -11.25455, -11.30909, -11.33091, -11.29818, -11.22182, -11.11273,

                    -10.90545, -10.70909, -10.48000, -10.29455, -10.08727, -9.84727,

                    -9.72727]

    # Convert sensitivity

    flux_50h = [math.pow(10.0,x) for x in sens_50h]

    flux_5h  = [math.pow(10.0,x) for x in sens_5h]

    flux_05h = [math.pow(10.0,x) for x in sens_05h]

    # Overplot sensitivities

    plt.loglog(energies_50h, flux_50h, 'b-.')

    plt.loglog(energies_5h,  flux_5h,  'g-.')

    plt.loglog(energies_05h, flux_05h, 'r-.')

    # Return

    return

# ===================== #

# 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

    # Check where file contains differential or integral sensitivity

    mode = "Integral"

    if sensitivity.has_key("emax"):

        emax_ref = -1.0

        for value in sensitivity["emax"]:

            if emax_ref < 0.0:

                emax_ref = value

            elif emax_ref != value:

                mode = "Differential"

                break

    # Add mode

    sensitivity["mode"] = mode

    # Add linear energy values

    if mode == "Differential":

        if sensitivity.has_key("loge"):

            name   = 'energy'

            values = []

            for value in sensitivity["loge"]:

                values.append(math.pow(10.0, value))

            sensitivity[name] = values

    else:

        if sensitivity.has_key("emin"):

            name   = 'energy'

            values = []

            for value in sensitivity["emin"]:

                values.append(value)

            sensitivity[name] = values

    # Return

    return sensitivity

# ===================== #

# Show sensitivity data #

# ===================== #

def show_sensitivity(sensitivity, filename):

    """

    Plot sensitivity.

    """

    # Build title

    title = sensitivity['mode'] + ' sensitivity ('+filename+')'

    # Create figure

    plt.figure()

    # Set plot attributes

    plt.title(title)

    plt.loglog()

    plt.grid()

    # Show differential sensitivity

    plt.plot(sensitivity['energy'], sensitivity['sensitivity\r'], 'ro-')

    # Set labels

    plt.xlabel("Energy (TeV)")

    if sensitivity['mode'] == "Differential":

        plt.ylabel(r"E$\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)")

    else:

        plt.ylabel(r"E$\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)")

    # Overplot sensitivities

    south_sensitivities()

    # 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]

    # Read sensitivity data

    sensitivity = read_sensitivity(filename)

    # Show sensitivity data

    show_sensitivity(sensitivity, filename)