GammaLib

#! /usr/bin/env python

import gammalib

import ctools

import cscripts

import numpy as np

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

# Main routine entry point #

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

if __name__ == '__main__':

    ra         = 0.0

    dec        = 0.0

    obs_offset = 1.0

    noff       = 4

    size       = 0.3 #degrees

    off_offset = 1.0

    nbins = 31

    duration = 1800

    emin     = 0.03

    emax     = 100

    rad      = 5.0

    caldb    = 'prod3b'

    irf      = 'South_z20_30m'

    pattern = cscripts.obsutils.set_obs_patterns("single", ra=ra, dec=dec, offset=obs_offset)

    obs = cscripts.obsutils.set_obs_list(pattern, duration=duration, emin=emin, emax=emax, rad=rad, caldb=caldb, irf=irf)

    xmlname  = './crab.xml'

    obs.models(gammalib.GModels(xmlname))

    obs = cscripts.obsutils.sim(obs)

    # Define ON region

    ondir = gammalib.GSkyDir()

    ondir.radec_deg(ra,dec)

    on    = gammalib.GSkyRegions()

    onreg = gammalib.GSkyRegionCircle(ondir,size)

    on.append(onreg)

    # Define OFF regions

    off = gammalib.GSkyRegions()

    ra_off  = [0,1,0,-1]

    dec_off = [1,0,-1,-1]

    for i in range(noff):

        offdir = gammalib.GSkyDir()

        offdir.radec_deg(ra_off[i], dec_off[i])

        offreg = gammalib.GSkyRegionCircle(offdir, size)

        off.append(offreg)

    # Define energy binning

    etrue = gammalib.GEbounds(nbins, gammalib.GEnergy(emin, "TeV"), gammalib.GEnergy(emax, "TeV"))

    ereco = gammalib.GEbounds(nbins, gammalib.GEnergy(emin, "TeV"), gammalib.GEnergy(emax, "TeV"))

    onoff_container = gammalib.GObservations()

    obs_onoff = gammalib.GCTAOnOffObservation(obs[0], etrue, ereco, on, off)

    obs_onoff.name("TestCrab")

    obs_onoff.id("0001")

    onoff_container.append(obs_onoff)

    onoff_container.models(gammalib.GModels('crab_fitmodel.xml'))

    # Set optimizer and logger

    opt = gammalib.GOptimizerLM()

    onoff_container.optimize(opt)

    onoff_container.errors(opt)

    print(onoff_container.models())

    fit_models = onoff_container.models()

    gradient  = gammalib.GVector()

    curvature = gammalib.GMatrixSparse()

    npred = []

    value = onoff_container[0].likelihood(fit_models, gradient, curvature, npred)