GammaLib

#! /Applications/Anaconda/bin/python

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

# This script tests the ON-OFF analysis functionality of the gammalib

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

from ctools import *

from gammalib import *

from math import *

from cscripts import obsutils

from copy import deepcopy

import os

import glob

import sys

import numpy

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

# Setup single observation #

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

def single_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \

            emin=0.1, emax=100.0, rad=5.0, \

            irf="South_50h", caldb="prod2", id="000000", name="Obs", instrument="CTA"):

    """

    Returns a single CTA observation.

    Parameters:

     pntdir   - Pointing direction [GSkyDir]

    Keywords:

     tstart     - Start time [seconds] (default: 0.0)

     duration   - Duration of observation [seconds] (default: 1800.0)

     deadc      - Deadtime correction factor (default: 0.95)

     emin       - Minimum event energy [TeV] (default: 0.1)

     emax       - Maximum event energy [TeV] (default: 100.0)

     rad        - ROI radius used for analysis [deg] (default: 5.0)

     irf        - Instrument response function (default: cta_dummy_irf)

     caldb      - Calibration database path (default: "dummy")

     id         - Run identifier (default: "000000")

     instrument - Intrument (default: "CTA")

    """

    # Allocate CTA observation

    obs_cta = GCTAObservation()

    # Set calibration database

    db = GCaldb()

    if (os.path.isdir(caldb)):

        db.rootdir(caldb)

    else:

        db.open("cta", caldb)

    # Set pointing direction

    pnt = GCTAPointing()

    pnt.dir(pntdir)

    obs_cta.pointing(pnt)

    # Set ROI

    roi     = GCTARoi()

    instdir = GCTAInstDir()

    instdir.dir(pntdir)

    roi.centre(instdir)

    roi.radius(rad)

    # Set GTI

    gti = GGti()

    gti.append(GTime(tstart), GTime(tstart+duration))

    # Set energy boundaries

    ebounds = GEbounds(GEnergy(emin, "TeV"),GEnergy(emax, "TeV"))

    # Allocate event list

    events = GCTAEventList()

    events.roi(roi)

    events.gti(gti)

    events.ebounds(ebounds)

    obs_cta.events(events)

    # Set instrument response

    obs_cta.response(irf, db)

    # Set ontime, livetime, and deadtime correction factor

    obs_cta.ontime(duration)

    obs_cta.livetime(duration*deadc)

    obs_cta.deadc(deadc)

    obs_cta.id(id)

    obs_cta.name(name)

    # Return CTA observation

    return obs_cta

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

# Add CTA background model to observations #

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

def add_background_model(models,name="Background",instru="CTA",id="",bgd_file="",bgd_sigma=0.0):

    """

    Add CTA background model to a model container.

    Parameters:

     models - a model container

     bgd_file - file function for background spectral dependence

     bgd_sigma - sigma parameter for the gaussian offset angle dependence of the background rate

    """

    # Gaussian radial dependence model

    if bgd_sigma > 0.0:

        # Spatial part

        bgd_radial = GCTAModelRadialGauss(bgd_sigma)

        # Spectral part

        if len(bgd_file) > 0:

            bgd_spectrum = GModelSpectralFunc(bgd_file,1.0)

        else:

            pivot_nrj=GEnergy(1.0e6, "MeV")

            bgd_spectrum = GModelSpectralPlaw(1.0e-6, -2.0, pivot_nrj)

            bgd_spectrum["Prefactor"].value(61.8e-6)

            bgd_spectrum["Prefactor"].scale(1.0e-6)

            bgd_spectrum["PivotEnergy"].value(1.0)

            bgd_spectrum["PivotEnergy"].scale(1.0e6)

            bgd_spectrum["Index"].value(-1.85)

            bgd_spectrum["Index"].scale(1.0)

            # Full model

            bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)

    # IRF-defined model

    else:

            # Spectral model

            pivot_nrj=GEnergy(1.0e6, "MeV")

            bgd_spectrum = GModelSpectralPlaw(1.0e-6, -2.0, pivot_nrj)

            bgd_spectrum["Prefactor"].value(1.0)

            bgd_spectrum["Prefactor"].scale(1.0)

            bgd_spectrum["PivotEnergy"].value(1.0)

            bgd_spectrum["PivotEnergy"].scale(1.0e6)

            bgd_spectrum["Index"].value(0.0)

            bgd_spectrum["Index"].scale(1.0)

            # Full model

            bgd_model = GCTAModelIrfBackground(bgd_spectrum)

    # Set other background properties

    if (len(name) > 0):

        bgd_model.name(name)

    else:

        bgd_model.name("Background")

    if (len(instru) > 0):

        bgd_model.instruments(instru)

    else:

        bgd_model.instruments("CTA")

    if (len(id) > 0):

        bgd_model.ids(id)

    # By default, set model to be fitted

    if (len(bgd_file) > 0):

        bgd_model['Normalization'].free()

    else:

        bgd_model['Prefactor'].free()

        bgd_model['Index'].free()

    # Add background model to container

    models.append(bgd_model)

    # Return model container

    return models

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

# Set any general type source #

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

def add_source_model(models, name, ra, dec, coord='CEL', \

                     skymap='', specfile='', \

                     flux=5.7e-16, index=-2.5, pivot=3.0e5, \

                     type="point", sigma=1.0, radius=1.0, width=0.1, \

                     nnode=0, emin=1.0, emax=10.0, \

                     pivot_scale=1e6, flux_scale=1.0e-16):

    """

    Adds a single point-source with power-law spectrum to a model container.

    The default is crab-like.

    Parameters:

     models - a model container

     name   - Unique source name

     ra     - RA of source location [deg]

     dec    - Declination of source location [deg]

    Keywords:

     flux   - Source flux density in 1.0e-16 ph/cm2/s/MeV

     index  - Source spectral index

     pivot  - Pivot energy in TeV

     type   - Source spatial model

     sigma/radius/width  - Spatial model parameters

     nnode  - Number of nodes (if > 0, spectral model is a node function)

    """

    # Set source spectral model

    if (nnode > 0) and (emin != emax):

        spectrum = GModelSpectralNodes()

        dloge=log10(emax/emin)/nnode

        logemin=log10(emin)

        spectrum.reserve(nnode)

        for i in range(nnode):

            enode = GEnergy()

            enode.TeV(10.0**(logemin+(i+0.5)*dloge))

            inode=flux*(enode.MeV()/pivot)**(index)

            spectrum.append(enode, inode)

            spectrum[i*2].scale(pivot_scale)

            spectrum[i*2].fix()

            spectrum[i*2+1].scale(flux_scale)

            spectrum[i*2+1].free()    

    # Set source spectral model

    elif (len(specfile) > 0):

        spectrum = GModelSpectralFunc(specfile)

        spectrum["Normalization"].value(1.0)

        spectrum["Normalization"].free()

    else:

        pivot_nrj=GEnergy()

        pivot_nrj.TeV(pivot*pivot_scale/1e6)

        spectrum = GModelSpectralPlaw(flux, index,pivot_nrj)

        spectrum["Prefactor"].scale(flux_scale)

        spectrum["Prefactor"].value(flux)

        spectrum["Prefactor"].free()

        spectrum["PivotEnergy"].scale(pivot_scale)

        spectrum["PivotEnergy"].value(pivot)

        spectrum["Index"].value(index)

        spectrum["Index"].scale(1.0)

        spectrum["Index"].free()        

    # Set source spatial model

    location = GSkyDir()

    if (coord == 'CEL'):

        location.radec_deg(ra, dec)

    else:

        location.lb_deg(ra, dec)    

    if (len(skymap) > 0):

        spatial = GModelSpatialDiffuseMap(skymap,1.0)

        source = GModelSky(spatial, spectrum)

        spatial[0].free()

    elif type == "point":

        spatial = GModelSpatialPointSource(location)

        spatial[0].free()

        spatial[1].free()

        source  = GModelSky(spatial, spectrum)

    elif type == "gauss":

        radial = GModelSpatialRadialGauss(location, sigma)

        radial[0].free()

        radial[1].free()

        radial[2].free()

        source = GModelSky(radial, spectrum)

    elif type == "disk":

        radial = GModelSpatialRadialDisk(location, radius)

        radial[0].free()

        radial[1].free()

        radial[2].free()

        source = GModelSky(radial, spectrum)

    elif type == "shell":

        radial = GModelSpatialRadialShell(location, radius, width)

        radial[0].free()

        radial[1].free()

        radial[2].free()

        radial[3].free()

        source = GModelSky(radial, spectrum)

    else:

        print "ERROR: Unknown source type '"+type+"'."

        return None

    # Set source name

    source.name(name)

    # Append source model to model container

    models.append(source)

    # Return source

    return models

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

# Set the fit parameters for a given model #

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

def set_fit_param(model, parname=[], parfit=[], parval=[]):

    """

    Set the fit parameters for a given model

    Arguments

    - model: a model container

    - parname: array of parameters names

    - parfit: array of fit/fix flags (True=fit,False=fix)

    - parval: array of initial/fix parameters values

    Notes

    - Arrays should have same numbers of elements

    - A way to not set a value but just the fit/fix flag is to pass/leave parval empty

    - If among several parameters some values must be set and others not, separate these and call the function two times

    """

    # Count/check number of parameters to be set

    if (len(parname) != len(parfit)):

        print "Incorrect input. Fit parameters not set or modified."

        return 0

    else:

        num_toset=len(parname)

    # Loop over parameters

    for i in range(model.size()):

        name=model[i].name()

        if (name in parname):

            # Get index of matching name

            i_par=parname.index(name)

            # Set the parameter to be fix/free

            if (parfit[i_par]):

                model[name].free()

            else:

                model[name].fix()

            # Optionally set the value

            if (len(parval) > 0): 

                model[name].value(parval[i_par])

    # Exit point

    return model

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

# Simulate observations #

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

def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,

        emin=None, emax=None, nbins=0, addbounds=False,

        binsz=0.05, npix=200, proj='TAN', coord='GAL'):

    """

    Simulate events for all observations in the container

    Simulate events for all observations using ctobssim. If the number of

    energy bins is positive, the events are binned into a counts cube using

    ctbin. If multiple observations are simulated, the counts cube is a

    stacked cube and the corresponding response cubes are computed using

    ctexpcube, ctpsfcube, ctbkgcube and optionally ctedispcube. The response

    cubes are attached to the first observation in the container, which

    normally is the observation with the counts cube.

    Parameters

    ----------

    obs : `~gammalib.GObservations`

        Observation container without events

    log : bool, optional

        Create log file(s)

    debug : bool, optional

        Create console dump?

    chatter : int, optional

        Chatter level

    edisp : bool, optional

        Apply energy dispersion?

    seed : int, optional

        Seed value for simulations

    emin : float, optional

        Minimum energy of counts cube for binned (TeV)

    emax : float, optional

        Maximum energy of counts cube for binned (TeV)

    nbins : int, optional

        Number of energy bins (0=unbinned)

    addbounds : bool, optional

        Add boundaries at observation energies

    binsz : float, optional

        Pixel size for binned simulation (deg/pixel)

    npix : int, optional

        Number of pixels in X and Y for binned simulation

    proj : str, optional

        Projection for binned simulation

    coord : str, optional

        Coordinate system for binned simulation

    Returns

    -------

    obs : `~gammalib.GObservations`

        Observation container filled with simulated events

    """

    # Allocate ctobssim application and set parameters

    sim = ctobssim(obs)

    sim['seed']    = seed

    sim['edisp']   = edisp

    sim['chatter'] = chatter

    sim['debug']   = debug

    # Optionally open the log file

    if log:

        sim.logFileOpen()

    # Run ctobssim application. This will loop over all observations in the

    # container and simulation the events for each observation. Note that

    # events are not added together, they still apply to each observation

    # separately.

    sim.run()

    # Binned option?

    if nbins > 0:

        # If energy boundaries are not given then determine the minimum and

        # the maximum energies from all observations and use these values

        # as energy boundaries. The energy boundaries are given in TeV.

        if emin == None or emax == None:

            emin = 1.0e30

            emax = 0.0

            for run in sim.obs():

                emin = min(run.events().ebounds().emin().TeV(), emin)

                emax = max(run.events().ebounds().emax().TeV(), emax)

        # Allocate ctbin application and set parameters

        bin = ctbin(sim.obs())

        bin['ebinalg']  = 'LOG'

        bin['emin']     = emin

        bin['emax']     = emax

        bin['enumbins'] = nbins

        bin['usepnt']   = True # Use pointing for map centre

        bin['nxpix']    = npix

        bin['nypix']    = npix

        bin['binsz']    = binsz

        bin['coordsys'] = coord

        bin['proj']     = proj

        bin['chatter']  = chatter

        bin['debug']    = debug

        # Optionally open the log file

        if log:

            bin.logFileOpen()

        # Run ctbin application. This will loop over all observations in

        # the container and bin the events in counts maps

        bin.run()

        # If we have multiple input observations then create stacked response

        # cubes and append them to the observation

        if len(sim.obs()) > 1:

            # Get stacked response (use pointing for map centre)

            response = get_stacked_response(sim.obs(), None, None,

                                            binsz=binsz, nxpix=npix, nypix=npix,

                                            emin=emin, emax=emax, enumbins=nbins,

                                            edisp=edisp,

                                            coordsys=coord, proj=proj,

                                            addbounds=addbounds,

                                            log=log, debug=debug,

                                            chatter=chatter)

            # Set stacked response

            if edisp:

                bin.obs()[0].response(response['expcube'],

                                      response['psfcube'],

                                      response['edispcube'],

                                      response['bkgcube'])

            else:

                bin.obs()[0].response(response['expcube'],

                                      response['psfcube'],

                                      response['bkgcube'])

            # Set new models

            bin.obs().models(response['models'])

        # Make a deep copy of the observation that will be returned

        # (the ctbin object will go out of scope one the function is

        # left)

        obs = bin.obs().copy()

    else:

        # Make a deep copy of the observation that will be returned

        # (the ctobssim object will go out of scope one the function is

        # left)

        obs = sim.obs().copy()

    # Delete the simulation

    del sim

    # Return observation container

    return obs

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

# Get stacked response #

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

def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,

                         emin=0.1, emax=100.0, enumbins=20, edisp=False,

                         coordsys='GAL', proj='TAN', addbounds=False,

                         log=False, debug=False, chatter=2):

    """

    Get stacked response cubes

    The number of energies bins are set to at least 30 bins per decade, and

    the "enumbins" parameter is only used if the number of bins is larger

    than 30 bins per decade.

    If the xref or yref arguments are "None" the response cube centre will be

    determined from the pointing information in the observation container.

    Parameters

    ----------

    obs : `~gammalib.GObservations`

        Observation container

    xref : float

        Right Ascension or Galactic longitude of response centre (deg)

    yref : float

        Declination or Galactic latitude of response centre (deg)

    binsz : float, optional

        Pixel size (deg/pixel)

    nxpix : int, optional

        Number of pixels in X direction

    nypix : int, optional

        Number of pixels in Y direction

    emin : float, optional

        Minimum energy (TeV)

    emax : float, optional

        Maximum energy (TeV)

    enumbins : int, optional

        Number of energy bins

    edisp : bool, optional

        Apply energy dispersion?

    coordsys : str, optional

        Coordinate system

    proj : str, optional

        Projection

    addbounds : bool, optional

        Add boundaries at observation energies

    log : bool, optional

        Create log file(s)

    debug : bool, optional

        Create console dump?

    chatter : int, optional

        Chatter level

    Returns

    -------

    result : dict

        Dictionary of response cubes

    """

    # If no xref and yref arguments have been specified then use the pointing

    # information

    if xref == None or yref == None:

        usepnt = True

    else:

        usepnt = False

    # Set number of energy bins to at least 30 per energy decade

    _enumbins = int((log10(emax) - log10(emin)) * 30.0)

    if enumbins > _enumbins:

        _enumbins = enumbins

    # Compute spatial binning for point spread function and energy dispersion

    # cubes. The spatial binning is 10 times coarser than the spatial binning

    # of the exposure and background cubes. At least 2 spatial are required.

    psf_binsz = 10.0 * binsz

    psf_nxpix = max(nxpix // 10, 2)  # Make sure result is int

    psf_nypix = max(nypix // 10, 2)  # Make sure result is int

    # Create exposure cube

    expcube = ctexpcube(obs)

    expcube['incube']    = 'NONE'

    expcube['usepnt']    = usepnt

    expcube['ebinalg']   = 'LOG'

    expcube['binsz']     = binsz

    expcube['nxpix']     = nxpix

    expcube['nypix']     = nypix

    expcube['enumbins']  = _enumbins

    expcube['emin']      = emin

    expcube['emax']      = emax

    expcube['coordsys']  = coordsys

    expcube['proj']      = proj

    expcube['addbounds'] = addbounds

    expcube['debug']     = debug

    expcube['chatter']   = chatter

    if not usepnt:

        expcube['xref'] = xref

        expcube['yref'] = yref

    if log:

        expcube.logFileOpen()

    expcube.run()

    # Create point spread function cube

    psfcube = ctpsfcube(obs)

    psfcube['incube']    = 'NONE'

    psfcube['usepnt']    = usepnt

    psfcube['ebinalg']   = 'LOG'

    psfcube['binsz']     = psf_binsz

    psfcube['nxpix']     = psf_nxpix

    psfcube['nypix']     = psf_nypix

    psfcube['enumbins']  = _enumbins

    psfcube['emin']      = emin

    psfcube['emax']      = emax

    psfcube['coordsys']  = coordsys

    psfcube['proj']      = proj

    psfcube['addbounds'] = addbounds

    psfcube['debug']     = debug

    psfcube['chatter']   = chatter

    if not usepnt:

        psfcube['xref'] = xref

        psfcube['yref'] = yref

    if log:

        psfcube.logFileOpen()

    psfcube.run()

    # Create background cube

    bkgcube = ctbkgcube(obs)

    bkgcube['incube']    = 'NONE'

    bkgcube['usepnt']    = usepnt

    bkgcube['ebinalg']   = 'LOG'

    bkgcube['binsz']     = binsz

    bkgcube['nxpix']     = nxpix

    bkgcube['nypix']     = nypix

    bkgcube['enumbins']  = _enumbins

    bkgcube['emin']      = emin

    bkgcube['emax']      = emax

    bkgcube['coordsys']  = coordsys

    bkgcube['proj']      = proj

    bkgcube['addbounds'] = addbounds

    bkgcube['debug']     = debug

    bkgcube['chatter']   = chatter

    if not usepnt:

        bkgcube['xref'] = xref

        bkgcube['yref'] = yref

    if log:

        bkgcube.logFileOpen()

    bkgcube.run()

    # If energy dispersion is requested then create energy dispersion cube

    if edisp:

        edispcube = ctedispcube(obs)

        edispcube['incube']    = 'NONE'

        edispcube['usepnt']    = usepnt

        edispcube['ebinalg']   = 'LOG'

        edispcube['binsz']     = psf_binsz

        edispcube['nxpix']     = psf_nxpix

        edispcube['nypix']     = psf_nypix

        edispcube['enumbins']  = _enumbins

        edispcube['emin']      = emin

        edispcube['emax']      = emax

        edispcube['coordsys']  = coordsys

        edispcube['proj']      = proj

        edispcube['addbounds'] = addbounds

        edispcube['debug']     = debug

        edispcube['chatter']   = chatter

        if not usepnt:

            edispcube['xref'] = xref

            edispcube['yref'] = yref

        if log:

            edispcube.logFileOpen()

        edispcube.run()

    # Build response dictionary

    response = {}

    response['expcube'] = expcube.expcube().copy()

    response['psfcube'] = psfcube.psfcube().copy()

    response['bkgcube'] = bkgcube.bkgcube().copy()

    response['models']  = bkgcube.models().copy()

    if edisp:

        response['edispcube'] = edispcube.expcube().copy()

    # Return response cubes

    return response

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

# Get stacked observation container #

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

def get_stacked_obs(cls, obs):

    """

    Bin an observation and return an observation container with a single

    binned observation

    Parameters

    ----------

    cls : `~cscript`

        cscript class

    obs : `~gammalib.GObservations`

        Observation container

    Returns

    -------

    obs : `~gammalib.GObservations`

        Observation container where the first observation is a binned observation

    """

    # Write header

    if cls._logExplicit():

        cls._log.header3('Binning events')

    # Bin events

    cntcube = ctbin(obs)

    cntcube['usepnt']   = False

    cntcube['ebinalg']  = 'LOG'

    cntcube['xref']     = cls['xref'].real()

    cntcube['yref']     = cls['yref'].real()

    cntcube['binsz']    = cls['binsz'].real()

    cntcube['nxpix']    = cls['nxpix'].integer()

    cntcube['nypix']    = cls['nypix'].integer()

    cntcube['enumbins'] = cls['enumbins'].integer()

    cntcube['emin']     = cls['emin'].real()

    cntcube['emax']     = cls['emax'].real()

    cntcube['coordsys'] = cls['coordsys'].string()

    cntcube['proj']     = cls['proj'].string()

    cntcube.run()

    # Write header

    if cls._logExplicit():

        cls._log.header3('Creating stacked response')

    # Get stacked response

    response = get_stacked_response(obs,

                                    cls['xref'].real(),

                                    cls['yref'].real(),

                                    binsz=cls['binsz'].real(),

                                    nxpix=cls['nxpix'].integer(),

                                    nypix=cls['nypix'].integer(),

                                    emin=cls['emin'].real(),

                                    emax=cls['emax'].real(),

                                    enumbins=cls['enumbins'].integer(),

                                    edisp=cls['edisp'].boolean(),

                                    coordsys=cls['coordsys'].string(),

                                    proj=cls['proj'].string())

    # Retrieve a new oberservation container

    new_obs = cntcube.obs().copy()

    # Set stacked response

    if cls['edisp'].boolean():

        new_obs[0].response(response['expcube'], response['psfcube'],

                            response['edispcube'], response['bkgcube'])

    else:

        new_obs[0].response(response['expcube'], response['psfcube'],

                            response['bkgcube'])

    # Get new models

    models = response['models']

    # Set models for new oberservation container

    new_obs.models(models)

    # Return new oberservation container

    return new_obs

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

# Make a simple model #

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

def make_simple_model(models,name,ra,dec,flux,idx,num_obs,instru='CTA'):

    # Add one background model for each observation

    for i_obs in range(num_obs):

      model_name="Background %d" % i_obs

      models = add_background_model(models,bgd_sigma=0.0,instru=instru,id="%d" % i_obs,name=model_name)

      set_fit_param(models[model_name],parname=['Sigma','Normalization','Prefactor','Index'],parfit=[False,True,True,False])

    # Add source model

    model_name="Test source"

    models = add_source_model(models, name, ra, dec, flux=flux, index=idx, pivot=1.0e6, type="point", pivot_scale=1e6, flux_scale=1.0e-16)

    set_fit_param(models[model_name],parname=['Prefactor','Index','RA','DEC'],parfit=[True,True,False,False])

    # Exit

    return models

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

# Make a simple model #

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

def reset_simple_model(models,ra,dec,flux,idx,num_obs):

    # Add background model (one for each pointing)

    for i_obs in range(num_obs):

      model_name="Background %d" % i_obs

      set_fit_param(models[model_name],parname=['Sigma','Normalization','Prefactor','Index'],parfit=[False,True,True,False],parval=[4.0,1.0,1.0,0.0])

    set_fit_param(models["Test source"],parname=['Prefactor','Index','RA','DEC'],parfit=[True,True,False,False],parval=[flux,idx,ra,dec])

    # Exit

    return models

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

# Print out fit results #

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

def print_fit_results(obslist):

    # Print

    num_obs=len(obslist)

    for i_obs in range(num_obs):

      model_name="Background %d" % i_obs

      print "%s prefactor: %.3e +/- %.3e " % (model_name,obslist.models()[model_name]['Prefactor'].value(),obslist.models()[model_name]['Prefactor'].error())

      print "%s index: %.3e +/- %.3e " % (model_name,obslist.models()[model_name]['Index'].value(),obslist.models()[model_name]['Index'].error())

    print "Source prefactor: %.3e +/- %.3e " % (obslist.models()["Test source"]['Prefactor'].value(),obslist.models()["Test source"]['Prefactor'].error())

    print "Source index: %.3e +/- %.3e " % (obslist.models()["Test source"]['Index'].value(),obslist.models()["Test source"]['Index'].error())

    print "Fit ended with value %.3e and status %d after %d iterations." % (opt.value(),opt.status(),opt.iter())

    # Exit

    return

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

# Routine entry point #

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

if __name__ == '__main__':

    """

    Aims:

    This script tests the ON-OFF analysis functionality of the gammalib

    Arguments:

    - None

    Notes:

    - Hard-coded source, observations, and regions (more flexible script later...)

    - Approximate method to set OFF regions (works only on equator)

    - Have to set up path to data base (db and irf)

    - Source flux set as flux density at 1TeV

    """

    # Job parameters

    # CTA

    db="prod2"

    irf="South_50h"

    bgd=""

    duration=1800.0

    rad=5.0

    # Source

    name="Test source"

    ra=0.0

    dec=0.0

    flux=5.0e-17

    idx=-2.5

    dflux=1.5  # Offset factor by which true value is rescaled before fit

    didx=0.3   # Offset factor by which true value is rescaled before fit

    # Observations

    num_obs=2

    emin=0.1

    emax=100.0

    nbins=9

    onshift=1.0

    onsize=0.3

    noff=3

    nx=500

    ny=500

    dx=0.02

    dy=0.02

    npix=nx*ny

    # Others

    chatter=4

    seed=10

    # Debug

    print 'Got params...'

    # Create observations container

    obslist = GObservations()

    dirlist = []

    for i_obs in range(num_obs):

      obsdir = GSkyDir()

      obsdir.radec_deg(ra+((-1)**i_obs)*onshift, dec)

      obs=single_obs(obsdir, tstart=0.0, duration=duration, deadc=0.95, \

                   emin=emin, emax=emax, rad=rad, \

                   irf=irf, caldb=db, id="%d" % i_obs, name="Obs%d" % i_obs, instrument="CTA")

      obslist.append(obs)

      dirlist.append(obsdir)

    print 'Created %d observations...' % (num_obs)

    # Create model container

    models = GModels()   

    models_onoff_fit = GModels()   

    # Make a simple model for two observations

    models = make_simple_model(models,name,ra,dec,flux,idx,num_obs,instru='CTA')

    models_onoff_fit = make_simple_model(models_onoff_fit,name,ra,dec,flux,idx,num_obs,instru='CTAOnOff')

    # Add model to the container

    obslist.models(models)

    obslist.models().save("sim_model.xml")

    print 'Created source and background model...'    

    # Make events list

    obslist = sim(obslist, log=False, debug=False, chatter=chatter, edisp=False, seed=seed, nbins=0)

    # Save events and observations

    for obs in obslist:

        obs.save("%s_events.fits" % obs.name())

    obslist.save("UNBINNED-observations.xml")

    print 'Created events list...' 

    # Make counts cube

    bin = ctbin(obslist)

    bin['ebinalg']  = 'LOG'

    bin['emin']     = emin

    bin['emax']     = emax

    bin['enumbins'] = nbins

    bin['usepnt']   = True

    bin['nxpix']    = 200

    bin['nypix']    = 200

    bin['binsz']    = 0.05

    bin['coordsys'] = 'CEL'

    bin['proj']     = 'TAN'

    bin['outcube']  = "cntcube.fits"

    bin.run()

    bin.save()

    print 'Created events cube...'    

    # Make ON map

    srcdir = GSkyDir()

    srcdir.radec_deg(ra,dec)

    onmap= GSkyMap('TAN','CEL',ra,dec,dx,dy,nx,ny,1)

    # Loop over all pixels in map and flag those falling into regions

    for i in range(npix):

        pixdir=onmap.inx2dir(i)

        if pixdir.dist_deg(srcdir) < onsize:

            onmap[i,0]= 1.0

    # Make ON region

    onreg=GSkyRegionMap(onmap)

    print 'Created ON map'

    # Make OFF maps

    offregs=[]

    for i_obs in range(num_obs):

      offmap= GSkyMap('TAN','CEL',ra,dec,dx,dy,nx,ny,1)

      offdirs=[]

      # Prepare directions to circular OFF regions

      for i in range(noff):

          phi=(i+1)*360./(noff+1.0)

          offdir = GSkyDir(dirlist[i_obs])

          offdir.rotate_deg(phi-((-1)**i_obs)*90., onshift)

          offdirs.append(offdir)

      # Loop over all pixels in map and flag those falling into OFF regions    

      for i in range(npix):

          # Direction to pixel in map

          pixdir=offmap.inx2dir(i)

          # Loop over centers of OFF regions and compute distance to this pixel

          for cendir in offdirs:

              if pixdir.dist_deg(cendir) < onsize:

                  offmap[i,0]= 1.0

                  break

      # Make a sky region out of the map

      offregs.append(GSkyRegionMap(offmap))

    # Define energy binning for ON-OFF observations

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

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

    # Make and save ON-OFF observations

    onofflist = GObservations()

    for obs,offreg,i_obs in zip(obslist,offregs,range(num_obs)):

      onoffobs=GCTAOnOffObservation(obs, etrue, ereco, onreg, offreg, srcdir)

      onoffobs.name("Obs%d" % i_obs)

      onoffobs.id("%d" % i_obs)

      onoffobs.arf().save("ARF_%d.fits" % i_obs)

      onoffobs.rmf().save("RMF_%d.fits" % i_obs)

      onoffobs.on_spec().save("ON-spectrum_%d.fits" % i_obs)

      onoffobs.off_spec().save("OFF-spectrum_%d.fits" % i_obs)

      onoffobs.on_regions().map().save("ON-regions_%d.fits" % i_obs)

      onoffobs.off_regions().map().save("OFF-regions_%d.fits" % i_obs)

      onofflist.append(onoffobs)

    # Save ON-OFF observations

    onofflist.save("ONOFF-observations.xml")

    print "Created and saved ON-OFF observations"

    # Print ON-OFF data

    for oneobs in onofflist:

      print "\nObservation %s with %.1fs ON time" % (oneobs.name(),oneobs.ontime())

      for i in range(nbins):

        print "ONobs %d - OFFobs %d - Aeff %.2e - Bgd %.2e" % (oneobs.on_spec()[i],oneobs.off_spec()[i],oneobs.arf()[i],oneobs.arf()['BACKGROUND'][i])

      print '\n'   

    # Reset model container and append it to observation list

    models_onoff_fit=reset_simple_model(models_onoff_fit,ra,dec,dflux*flux,idx-didx,num_obs)

    onofflist.models(models_onoff_fit)

    # Set optimizer and logger

    log=GLog('fit_onoff.log',True)

    log.chatter(chatter)

    opt=GOptimizerLM(log)

    # Optimize for ON-OFF analysis

    print '\nON-OFF fitting...'

    onofflist.optimize(opt)

    onofflist.errors(opt)

    print_fit_results(onofflist)

    # Free memory

    del opt

    del log

    # Reset model container and append it to observation list

    models=reset_simple_model(models,ra,dec,dflux*flux,idx-didx,num_obs)

    obslist.models(models)

    # Set optimizer and logger

    log=GLog('fit_unbinned.log',True)

    log.chatter(chatter)

    opt=GOptimizerLM(log)

    # Optimize for Fermi-style analysis

    print '\nUnbinned fitting...'

    obslist.optimize(opt)

    obslist.errors(opt)

    print_fit_results(obslist)

    # Free memory

    del opt

    del log