GammaLib

#! /usr/bin/env python

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

# Test file function

#

# Copyright (C) 2017 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 gammalib

import math

import matplotlib.pyplot as plt

# ========= #

# Make plot #

# ========= #

def make_plot(tmin, tmax, rate=100.0, bins=100):

    """

    Make plot

    """

    # Load temporal file function

    filefct = gammalib.GModelTemporalFunc("temp_filefunction.fits")

    # Set reference time (should be able to extract that from file function?)

    tref = gammalib.GTimeReference(51544.5,"s")

    # Set random number generator

    ran = gammalib.GRan()

    # Create times

    times = filefct.mc(rate, gammalib.GTime(tmin, tref),

                             gammalib.GTime(tmax, tref), ran)

    # Convert times to seconds

    seconds = [time.convert(tref) for time in times]

    # Create time vector for eval

    norm       = rate/float(bins)

    eval_time  = [float(i) for i in range(1800)]

    eval_value = [filefct.eval(gammalib.GTime(t, tref))*norm for t in eval_time]

    # Create figure

    plt.figure(1)

    # Plot file function

    plt.plot(eval_time, eval_value)

    # Plot histogram

    plt.hist(seconds, bins=bins)

    # Show plot

    plt.show()

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Make plot

    #make_plot(0.0, 1800.0, rate=1000000.0)

    #make_plot(100.0, 1450.0, rate=1000000.0)

    #make_plot(410.0, 1000.0, rate=1000000.0)

    #make_plot(1178.0, 3000.0, rate=1000000.0)

    #make_plot(2000.0, 3000.0, rate=1000000.0)

    make_plot(-1000.0, -100.0, rate=1000000.0)