GammaLib

#! /usr/bin/env python

import gammalib

# Allocate response

rsp_dol = '/project-data/integral/ic/spi/rsp/spi_irf_grp_0021.fits'

rsp     = gammalib.GSPIResponse(rsp_dol)

print(rsp)

# Set response

og_dol  = '/Users/jurgen/git/gammalib/gammalib/inst/spi/test/data/obs/og_spi.fits'

obs     = gammalib.GSPIObservation(og_dol)

rsp.set(obs)

print(rsp)

# Save response

rsp.save('irf.fits', True)

# Loop over event bins

for bin in obs.events():

    # Print event bin

    #print(bin)

    # Get SPI_X direction for bin

    spi_x = obs.events().spi_x(bin.ipt())

    spi_z = obs.events().spi_z(bin.ipt())

    #print(spi_x)

    # Allocate sky map

    map = gammalib.GSkyMap('ARC', 'CEL', spi_x.ra_deg(), spi_x.dec_deg(), 0.1, 0.1, 2000, 1500)

    #print(map)

    # Compute response for all sky map pixels

    for i in range(map.npix()):

        # Get Sky direction

        dir = map.inx2dir(i)

        #print(dir)

        # Get IRF

        irf = rsp.irf(dir, bin, 0)

        #break

        # Set sky map bin

        map[i] = irf

    # Set x direction

    pixel = map.dir2pix(spi_x)

    print(pixel)

    for dx in range(10):

        p1 = map.pix2inx(gammalib.GSkyPixel(pixel.x()+dx, pixel.y()))

        p2 = map.pix2inx(gammalib.GSkyPixel(pixel.x()-dx, pixel.y()))

        p3 = map.pix2inx(gammalib.GSkyPixel(pixel.x(), pixel.y()+dx))

        p4 = map.pix2inx(gammalib.GSkyPixel(pixel.x(), pixel.y()-dx))

        map[p1] = 32000.0

        map[p2] = 32000.0

        map[p3] = 32000.0

        map[p4] = 32000.0

    # Set z direction

    pixel = map.dir2pix(spi_z)

    for dx in range(5):

        for dy in range(5):

            p1 = map.pix2inx(gammalib.GSkyPixel(pixel.x()+dx, pixel.y()+dy))

            p2 = map.pix2inx(gammalib.GSkyPixel(pixel.x()+dx, pixel.y()-dy))

            p3 = map.pix2inx(gammalib.GSkyPixel(pixel.x()-dx, pixel.y()+dy))

            p4 = map.pix2inx(gammalib.GSkyPixel(pixel.x()-dx, pixel.y()-dy))

            map[p1] = 32000.0

            map[p2] = 32000.0

            map[p3] = 32000.0

            map[p4] = 32000.0

    # Save sky map

    map.save('irf_map.fits', True)

    # Break

    break