ctools

#! /usr/bin/env python

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

# Show spill over for a given observation and model

#

# Copyright (C) 2018 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 os

import math

import gammalib

import ctools

import cscripts

import matplotlib.pyplot as plt

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

# Plot residual spectrum #

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

def plot_resspec(frame, energies, resspec, errspec):

    """

    Plot residual spectrum

    Parameters

    ----------

    frame : pyplot

        Frame for spectrum

    """

    # Set axes scales and limit

    frame.set_xscale('log')

    frame.set_xlim([0.5, 40])

    # Counts and model

    frame.errorbar(energies, resspec, yerr=errspec,

                   fmt='ko', capsize=0, linewidth=2, zorder=2, label='Data')

    frame.plot(frame.get_xlim(),[0,0], 'r-')

    frame.set_xlabel('Energy (TeV)')

    frame.set_ylabel('Stacked - Joint Binned')

    # Return

    return

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

# Plot residual map #

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

def plot_resmap(plt1, plt2, map, error, steps=1):

    """

    Plot pull histogram

    Parameters

    ----------

    plt1 : pyplot

        Frame for first map

    plt2 : pyplot

        Frame for second map

    map : `~gammalib.GSkyMap()`

        Sky map

    error : `~gammalib.GSkyMap()`

        Sky map

    steps : int

        Number of steps for rebinning

    """

    # Create longitude array from skymap

    x_lon = []

    y_lon = []

    e_lon = []

    for ix in range(0, map.nx(), steps):

        value = 0.0

        err   = 0.0

        for istep in range(steps):

            if istep+ix < map.nx():

                for iy in range(map.ny()):

                    index  = istep+ix+iy*map.nx()

                    value += map[index]

                    err   += error[index]*error[index]

        x_lon.append(ix)

        y_lon.append(value)

        e_lon.append(math.sqrt(err))

    # Create latitude array from skymap

    x_lat = []

    y_lat = []

    e_lat = []

    for iy in range(0, map.ny(), steps):

        value = 0.0

        err   = 0.0

        for istep in range(steps):

            if istep+iy < map.ny():

                for ix in range(map.nx()):

                    index  = ix+(istep+iy)*map.nx()

                    value += map[index]

                    err   += error[index]*error[index]

        x_lat.append(iy)

        y_lat.append(value)

        e_lat.append(math.sqrt(err))

    # Plot longitude distribution

    plt1.errorbar(x_lon, y_lon, yerr=e_lon,

                  fmt='ko', capsize=0, linewidth=2, zorder=2)

    plt1.axhline(0, color='r', linestyle='--')

    plt1.set_xlabel('Right Ascension (pixels)')

    plt1.set_ylabel('Stacked - Joint Binned')

    plt1.grid(True)

    # Plot latitude distribution

    plt2.errorbar(x_lat, y_lat, yerr=e_lat,

                  fmt='ko', capsize=0, linewidth=2, zorder=2)

    plt2.axhline(0, color='r', linestyle='--')

    plt2.set_xlabel('Declination (pixels)')

    plt2.set_ylabel('Stacked - Joint Binned')

    plt2.grid(True)

    # Return

    return

# ======== #

# Plot map #

# ======== #

def plot_map(frame, map, smooth=0.0):

    """

    Plot Aitoff map

    Parameters

    ----------

    frame : pyplot

        Frame for map

    map : `~gammalib.GSkyMap()`

        Sky map

    smooth : float, optional

        Map smoothing parameter (degrees)

    """

    # Optionally smooth map

    if smooth > 0.0:

        _map = map.copy()

        _map.smooth('DISK', smooth)

    else:

        _map = map

    # Create array from skymap

    array = []

    v_max = 0.0

    for iy in range(_map.ny()):

        row = []

        for ix in range(_map.nx()):

            index = ix+iy*_map.nx()

            value = _map[index]

            row.append(value)

        array.append(row)

    # Show Aitoff projection

    c = frame.imshow(array, aspect=1.0)

    cbar = plt.colorbar(c, orientation='vertical', shrink=0.7)

    frame.grid(True)

    frame.set_xlabel('Right Ascension (pixels)')

    frame.set_ylabel('Declination (pixels)')

    # Return

    return

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

# Determine spatial residuals #

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

def residual_spatial(modcube1, modcube2):

    """

    Determine spatial residuals

    Parameters

    ----------

    modcube1 : `~gammalib.GCTAEventCube`

        First model cube

    modcube2 : `~gammalib.GCTAEventCube`

        Second model cube

    Returns

    -------

    map, error : `~gammalib.GSkyMap()`

        Residual sky map and error

    """

    # Get model cube maps

    map1 = modcube1.counts()

    map2 = modcube2.counts()

    # Generate residual map

    map = map1 - map2

    map.stack_maps()

    # Generate residual error bars

    error = map2.copy()

    error.stack_maps()

    error = error.sqrt()

    # Return residual map and error map

    return map, error

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

# Determine spectral residuals #

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

def residual_spectral(modcube1, modcube2):

    """

    Determine spectral residuals

    Parameters

    ----------

    modcube1 : `~gammalib.GCTAEventCube`

        First model cube

    modcube2 : `~gammalib.GCTAEventCube`

        Second model cube

    Returns

    -------

    map, error : list

        Residual vector

    """

    # Get model cube maps

    map1 = modcube1.counts()

    map2 = modcube2.counts()

    # Compute spectra

    spec1    = []

    spec2    = []

    residual = []

    error    = []

    energies = []

    for ieng in range(map1.nmaps()):

        sum1 = 0.0

        sum2 = 0.0

        for ipix in range(map1.npix()):

            sum1 += map1[ipix, ieng]

            sum2 += map2[ipix, ieng]

        spec1.append(sum1)

        spec2.append(sum2)

        residual.append(sum1 - sum2)

        error.append(math.sqrt(sum2))

        energies.append(modcube2.energy(ieng).TeV())

    # Return residual and error

    return energies, residual, error

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

# Create model cube #

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

def create_modcube(obs, ebins=20, edisp=False):

    """

    Create model cube

    Parameters

    ----------

    obs : `~gammalib.GObservations`

        Observation container

    ebins : int, optional

        Number of energy bins for binned or stacked analysis

    edisp : boolean, optional

        Use energy dispersion

    """

    # Create model cube

    ctmodel = ctools.ctmodel(obs)

    ctmodel['edisp']    = edisp

    ctmodel['incube']   = 'NONE'

    ctmodel['ebinalg']  = 'LOG'

    ctmodel['emin']     = 0.67

    ctmodel['emax']     = 30.0

    ctmodel['enumbins'] = ebins

    ctmodel['coordsys'] = 'CEL'

    ctmodel['proj']     = 'CAR'

    ctmodel['xref']     = 83.63

    ctmodel['yref']     = 22.01

    ctmodel['nxpix']    = 50

    ctmodel['nypix']    = 50

    ctmodel['binsz']    = 0.02

    #ctmodel['debug']    = True

    ctmodel.run()

    #print(ctmodel.cube())

    # Return model cube

    return ctmodel.cube().copy()

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

# Create stacked observation #

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

def create_stacked(obs, ebins=20, edisp=False, enumbins=300):

    """

    Create stacked observation

    Parameters

    ----------

    obs : `~gammalib.GObservations`

        Observation container

    ebins : int, optional

        Number of energy bins for binned or stacked analysis

    edisp : boolean, optional

        Use energy dispersion

    """

    # Create counts cube

    ctbin = ctools.ctbin(obs)

    ctbin['stack']    = True

    ctbin['ebinalg']  = 'LOG'

    ctbin['emin']     = 0.67

    ctbin['emax']     = 30.0

    ctbin['enumbins'] = ebins

    ctbin['coordsys'] = 'CEL'

    ctbin['proj']     = 'CAR'

    ctbin['xref']     = 83.63

    ctbin['yref']     = 22.01

    ctbin['nxpix']    = 50

    ctbin['nypix']    = 50

    ctbin['binsz']    = 0.02

    ctbin.run()

    # Create exposure cube

    ctexpcube = ctools.ctexpcube(obs)

    ctexpcube['incube']  = 'NONE'

    ctexpcube['ebinalg']  = 'LOG'

    ctexpcube['emin']     =   0.1

    ctexpcube['emax']     = 100.0

    ctexpcube['enumbins'] = enumbins

    ctexpcube['coordsys'] = 'CEL'

    ctexpcube['proj']     = 'CAR'

    ctexpcube['xref']     = 83.63

    ctexpcube['yref']     = 22.01

    ctexpcube['nxpix']    = 50

    ctexpcube['nypix']    = 50

    ctexpcube['binsz']    = 0.02

    ctexpcube.run()

    # Create PSF cube

    ctpsfcube = ctools.ctpsfcube(obs)

    ctpsfcube['incube']    = 'NONE'

    ctpsfcube['ebinalg']   = 'LOG'

    ctpsfcube['emin']      =   0.1

    ctpsfcube['emax']      = 100.0

    ctpsfcube['enumbins']  = enumbins

    ctpsfcube['coordsys']  = 'CEL'

    ctpsfcube['proj']      = 'CAR'

    ctpsfcube['xref']      = 83.63

    ctpsfcube['yref']      = 22.01

    ctpsfcube['nxpix']     = 35

    ctpsfcube['nypix']     = 25

    ctpsfcube['binsz']     = 0.2

    ctpsfcube['amax']      = 0.7

    ctpsfcube['anumbins']  = 300

    ctpsfcube.run()

    # Create energy dispersion cube

    if edisp:

        ctedispcube = ctools.ctedispcube(obs)

        ctedispcube['incube']    = 'NONE'

        ctedispcube['ebinalg']   = 'LOG'

        ctedispcube['emin']      =   0.1 # Full energy range

        ctedispcube['emax']      = 100.0 # Full energy range

        ctedispcube['enumbins']  = enumbins

        ctedispcube['coordsys']  = 'CEL'

        ctedispcube['proj']      = 'CAR'

        ctedispcube['xref']      = 83.63

        ctedispcube['yref']      = 22.01

        ctedispcube['nxpix']     = 35

        ctedispcube['nypix']     = 25

        ctedispcube['binsz']     = 0.2

        ctedispcube['migramax']  = 5.0

        ctedispcube['migrabins'] = 300

        ctedispcube.run()

    # Create dummy background cube

    bkgcube = gammalib.GCTACubeBackground()

    # Get stacked observation

    obs_stacked = ctbin.obs().copy()

    # Build stacked observation

    if edisp:

        obs_stacked[0].response(ctexpcube.expcube(),

                                ctpsfcube.psfcube(),

                                ctedispcube.edispcube(),

                                bkgcube)

    else:

        obs_stacked[0].response(ctexpcube.expcube(),

                                ctpsfcube.psfcube(),

                                bkgcube)

    # Return stacked observation

    return obs_stacked

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

# Show spill over #

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

def show_spill_over(obsname, modelname, srcname, ebins=20, edisp=False):

    """

    """

    # Load observations

    obs_unbinned = gammalib.GObservations(obsname)

    # Create stacked observation

    obs_stacked = create_stacked(obs_unbinned, ebins=ebins, edisp=edisp)

    #print(obs_stacked)

    #print(obs_stacked[0].response())

    # Setup model

    models      = gammalib.GModels()

    models_full = gammalib.GModels(modelname)

    models.append(models_full[srcname])

    # Append model to observations

    obs_unbinned.models(models)

    obs_stacked.models(models)

    # Compute model cubes

    modcube_unbinned = create_modcube(obs_unbinned, ebins=ebins, edisp=edisp)

    modcube_stacked  = create_modcube(obs_stacked, ebins=ebins, edisp=edisp)

    # Create residual map

    resmap, errmap             = residual_spatial(modcube_stacked, modcube_unbinned)

    energies, resspec, errspec = residual_spectral(modcube_stacked, modcube_unbinned)

    # Create figure

    fig = plt.figure(figsize=(14,5))

    # Create subplot for residual spectrum

    frame1 = fig.add_subplot(131)

    plot_resspec(frame1, energies, resspec, errspec)

    # Create subplots for residual profiles

    frame2 = fig.add_subplot(232)

    frame3 = fig.add_subplot(235)

    plot_resmap(frame2, frame3, resmap, errmap)

    # Create subplot for residual map

    frame4 = fig.add_subplot(133)

    plot_map(frame4, resmap, smooth=0.0)

    # Show figure

    plt.show()

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Check if $HESSDATA has been set

    if 'HESSDATA' not in os.environ:

        print('HESSDATA environment variable not set. Please set HESSDATA to '

              'the root directory of the H.E.S.S. data.')

        exit(1)

    # Set parameters

    obsname   = 'obs_crab_selected.xml'

    modelname = 'crab_results_ptsrc_plaw_gauss_grad_hess.xml'

    srcname   = 'Crab'

    # Show spill over

    show_spill_over(obsname, modelname, srcname, ebins=40, edisp=False)