ctools

import gammalib

import ctools

import cscripts

import matplotlib.pyplot as plt

import numpy as np

from astropy.wcs import WCS

# Init observation

runlist = [20136]

cobs = cscripts.csiactobs()

cobs["prodname"] = "hess_dl3_dr1"

cobs["bkgpars"] = 2

cobs["psf_hiera"] = "psf_table"

cobs["debug"] = True

cobs["bkg_mod_hiera"] = "aeff"

cobs.runlist(runlist)

cobs.run()

observation = cobs.obs().copy()

# Select events

select = ctools.ctselect(observation)

select["usepnt"] = True

select["usethres"] = "DEFAULT"

select["rad"] = 2.5

select["emin"] = 0.2

select["emax"] = 92

select.run()

sel_observation = select.obs().copy()

# Build radial disk model

spec = gammalib.GModelSpectralPlaw(1e-14, -2.0, gammalib.GEnergy(1.0, 'TeV'))

spec['Prefactor'].fix()

spec['Index'].fix()

ra = 228.53

dec = -59.16

direction = gammalib.GSkyDir()

direction.radec_deg(ra,dec)

spat = gammalib.GModelSpatialRadialDisk(direction, 1.0)

spat["RA"].fix()

spat["DEC"].fix()

spat["Radius"].fix()

sourceModel = gammalib.GModelSky(spat, spec)

models = gammalib.GModels()

models.append(sourceModel)

sel_observation.models(models)

# Simulate observation

sim = ctools.ctobssim(sel_observation)

sim['debug'] = True

sim['nthreads'] = 40

sim['emin'] = 0.2

sim['emax'] = 91

sim['seed'] = 0

sim['tmin'] = "UNDEF"

sim.run()

sim_obs = sim.obs().copy()

# Bin observation

binpars = {'ebinalg': 'LOG', 'emin': 0.1, 'emax': 100,

           'enumbins': 24, 'nxpix': 150, 'nypix': 150,

           'binsz': 0.02, 'coordsys': "CEL", 'xref': 228.53,

           'yref': -59.16, 'proj': "CAR"}

cbin = ctools.ctbin(sim_obs.copy())

cbin.pardict(binpars)

cbin['debug'] = True

cbin.run()

cntcube = cbin.cube().copy()

# Model observation

mcube = ctools.ctmodel(sim_obs)

mcube.pardict(binpars)

mcube['incube'] = "NONE"

mcube['debug'] = True

mcube['nthreads'] = 30

mcube.run()

modcube = mcube.cube().copy()

# Compare counts and model cubes (note the difference in number of events)

print "Counts cube:"

print cntcube

print "\n Model cube:"

print modcube

# Generate resmap

cres = cscripts.csresmap(sim_obs)

cres.pardict(binpars)

cres['algorithm'] = "SUB"

cres.run()

# Define some helper functions for plotting

def get_wcs_from_skymap(skymap):

    from astropy.wcs import WCS

    import numpy as np

    pro = skymap.projection()

    wcs = WCS(naxis=2)

    wcs.wcs.crpix = [pro.crpix(0), pro.crpix(1)]

    wcs.wcs.cdelt = np.array([pro.cdelt(0), pro.cdelt(1)])

    wcs.wcs.crval = [pro.crval(0), pro.crval(1)]

    ctype_1 = "RA--" if pro.coordsys() == "EQU" else "GLON"

    ctype_1 += "-" + pro.code()

    ctype_2 = "DEC-" if pro.coordsys() == "EQU" else "GLAT"

    ctype_2 += "-" + pro.code()

    wcs.wcs.ctype = [ctype_1, ctype_2]

    return wcs

def plot_skymap(skymap, ax=None, label="Counts", cmap='afmhot', title="", correlate=False):

    skymap = skymap.copy()

    skymap.stack_maps()

    wcs = get_wcs_from_skymap(skymap)

    array = skymap.array()

    if correlate:

        array = tophat_correlate(array, pixelsize=wcs.wcs.cdelt[1], normalize_kernel=True)

    plot_array(array, wcs, ax, cmap=cmap, title=title, label=label)

def plot_array(array, wcs, ax=None, cmap='afmhot', title="", vmin=None, vmax=None, label=""):

    import matplotlib.pyplot as plt

    if ax== None:

        fig = plt.figure()

        ax = fig.add_subplot(1, 1, 1, projection=wcs)

    if wcs.wcs.ctype[0][:3] == "RA-":

        ax.set_xlabel('Right Ascension [deg]')

        ax.set_ylabel('Declination [deg]')

    else:

        ax.set_xlabel('Galactic logitude [deg]')

        ax.set_ylabel('Galactic latitude [deg]')

    ax.set_title(title)

    image = ax.imshow(array, origin='lower', interpolation='None', cmap=cmap, vmin=vmin, vmax=vmax)

    bar = ax.figure.colorbar(image, label=label)

def tophat_correlate( inarray, size=0.07, pixelsize=0.02, normalize_kernel=False ):

    """                                                                                         

    Correlates numpy array with tophat kernel                                                   

    :param inarray: Input array                                                                 

    :param size: size for kernel in degrees                                                     

    :param pixelsize: size of the pixels in input array in degrees                              

    """

    from astropy.convolution import Tophat2DKernel

    from astropy.convolution import convolve

    kernel = Tophat2DKernel(size/pixelsize)

    if not normalize_kernel:

        kernel = kernel.array/kernel.array[kernel.array.shape[0]/2,kernel.array.shape[0]/2]

    return convolve(inarray, kernel, normalize_kernel=normalize_kernel)

# Plot resmap

plot_skymap(cres._resmap, correlate=True)

plt.show()