ctools

#! /usr/bin/env python

import gammalib

import ctools as ct

import matplotlib as mp

import matplotlib.pyplot as plt

def _scatter(array1, array2, xtitle="", ytitle="", maintitle=""):

    plt.scatter(array1, array2, marker=".")

    plt.ylabel(ytitle)

    plt.xlabel(xtitle)

    plt.show()

def _plot_energy(filename):

    events = gammalib.GCTAEventList(filename)

    event_num = [events[i].event_id() for i in range(events.size())]

    energy_TeV = [events[i].energy().TeV() for i in range(events.size())]

    _scatter(event_num, energy_TeV, "Event Number", "Energy [TeV]")

def _plot_phase(filename):

    events    = gammalib.GCTAEventList(filename)

    fits      = gammalib.GFits("phased_events.fits")

    event_num = [events[i].event_id() for i in range(events.size())]

    times     = [events[i].time().mjd() for i in range(events.size())]

    phases    = [fits[1]["PHASE"][i] for i in range(fits[1]["PHASE"].length())]

    _scatter(times, phases, "MJD", "Phase")

def _plot_energy_time(filename):

    events = gammalib.GCTAEventList(filename)

    energy_TeV = [events[i].energy().TeV() for i in range(events.size())]

    times = [events[i].time().mjd() for i in range(events.size())]

    _scatter(times, energy_TeV, "MJD", "Energy [TeV]")

def _plot_phase_distribution(filename):

    fits      = gammalib.GFits(filename)

    phases    = [fits[1]["PHASE"][i] for i in range(fits[1]["PHASE"].length())]

    plt.hist(phases, 10)

    plt.xlim(0.,1.)

    plt.ylabel("Counts")

    plt.xlabel("Phase")

    plt.show()

if __name__ == '__main__':

    # Simulate some data with no phase information

    obs = ct.ctobssim()

    obs["inobs"] = ""

    obs["inmodel"] = "test_model.xml"

    obs["outevents"] = "events.fits"

    obs["caldb"] = "prod2"

    obs["irf"] = "South_0.5h"

    obs["ra"] = 83.63

    obs["dec"] = 22.01

    obs["rad"] = 2.5

    obs["tmin"] = 0

    obs["tmax"] = 1800

    obs["emin"] = 0.1

    obs["emax"] = 100.0

    obs.execute()

    # Append phase information to the events data

    phase = ct.ctphase()

    phase["inobs"] = "events.fits"

    phase["outobs"] = "phased_events.fits"

    phase["inmodel"] = "test_model.xml"

    phase["srcname"] = "Crab"

    phase["mjd"] = 51544.5

    phase["phase"] = 0.0

    phase["f0"] = 0.2345

    phase["f1"] = -1.0e-10

    phase["f2"] = -1.0e-13

    phase.logFileOpen()

    phase.execute()

    # Plot the phase values

    #_plot_energy("events.fits")

    #_plot_energy("phased_events.fits")

    _plot_phase("phased_events.fits")

    _plot_phase_distribution("phased_events.fits")

    #_plot_energy_time("events.fits")