ctools

#! /usr/bin/env python

import ctools

import gammalib

from math import *

import os

import glob

import sys

import time

import numpy

# Energy range for the simulation

sim_emin = 0.1

sim_emax = 100.

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

# Setup Galactic plane survey #

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

def survey_gplane2(lrange=6, lstep=2, brange=2, bstep=2, duration=36000):

    """

    Creates a single observation survey for test purposes.

    Keywords:

     lrange   - Longitude range (deg)

     lstep    - Longitude step size (deg)

     brange   - Latitude range (deg)

     bstep    - Latitude step size (deg)

     duration - Observation duration (integer seconds)

    """

    # Allocate observation container

    obs = gammalib.GObservations()

    # Loop over longitudes

    run_counter=0

    for l in numpy.linspace(-lrange,lrange,2*lrange/lstep+1):

        for b in numpy.linspace(-brange,brange,2*brange/bstep+1):

            # Set pointing

            pntdir = gammalib.GSkyDir()

            pntdir.lb_deg(l, b)

            run = ctools.obsutils.set_obs(pntdir, duration=duration, irf="irf_file.fits", caldb="$CALDB/data/cta/e/bcf/IFAE20120510_50h", emin=sim_emin, emax=sim_emax)

            run.id("run_counter_%i" %run_counter)

            run_counter+=1

            obs.append(run) 

    # Return observation container

    return obs

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

# Main routine entry point #

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

if __name__ == '__main__':

    # Set parameters

    outpath         = "TEST_out/"

    extension       = "490h"

    xml_model_name  = "obs_models_%s.xml" %(extension)

    xml_obs_name    = "%sobs_stacked%s.xml" %(outpath,extension)

    xml_smodel_name = "%sobs_models_%s.xml" %(outpath,extension)

    events_name     = "%sevents%s.xml" %(outpath,extension)

    cntmap_name     = "%scntmap%s.fits" %(outpath,extension)

    expcube_name    = "%sexpcube%s.fits" %(outpath,extension)

    psfcube_name    = "%spsfcube%s.fits" %(outpath,extension)

    bkgcube_name    = "%sbkgcube%s.fits" %(outpath,extension)

    likeout_model   = "%scube_results%s.xml" %(outpath,extension)

    model_name      = "%smodmap%s.fits" %(outpath,extension)

    nxpix           = 100

    nypix           = 100

    binsz           = 0.06

    # Setup observation survey

    obs = survey_gplane2(lrange=2.5, lstep=1, brange=3.5, bstep=1, duration=100)

    obs.models(gammalib.GModels(xml_model_name))

    # Run ctobssim (simulates event lists)

    print("Simulate events")

    sim = ctools.ctobssim(obs)

    sim.logFileOpen()

    sim["outfile"].filename(events_name)

    sim["prefix"].string("%ssim_events_%s" %(outpath,extension))

    sim.execute()

    print("Done")

    # Run ctbin (builds stacked event cube)

    print("\nBin events")

    bin = ctools.ctbin()

    bin.logFileOpen()

    bin["evfile"].filename(events_name)

    bin["outfile"].filename(cntmap_name)

    bin["ebinalg"].string("LOG")

    bin["emin"].real(0.1)

    bin["emax"].real(100.0)

    bin["enumbins"].integer(20)

    bin["usepnt"].boolean(False)

    bin["nxpix"].integer(nxpix)

    bin["nypix"].integer(nypix)

    bin["binsz"].real(binsz)

    bin["coordsys"].string("GAL")

    bin["proj"].string("CAR")

    bin["xref"].real(0.0)

    bin["yref"].real(0.0)

    bin.execute()

    print("Done")

    # Run ctexpcube (builds exposure cube)

    print("\nBuild exposure cube")

    expcube = ctools.ctexpcube(sim.obs())

    expcube.logFileOpen()

    expcube["cntmap"].filename(cntmap_name)

    expcube["outfile"].filename(expcube_name)

    #expcube["caldb"].string(caldb)    # Not needed, response in XML file

    #expcube["irf"].string(irf)        # Not needed, response in XML file

    expcube.execute()

    print("Done")

    # Run ctpsfcube (builds PSF cube)

    print("\nBuild PSF cube")

    psfcube = ctools.ctpsfcube(sim.obs())

    psfcube.logFileOpen()

    psfcube["cntmap"].filename("NONE")

    psfcube["outfile"].filename(psfcube_name)

    #psfcube["caldb"].string(caldb)    # Not needed, response in XML file

    #psfcube["irf"].string(irf)        # Not needed, response in XML file

    psfcube["ebinalg"].string("LOG")

    psfcube["emin"].real(0.1)

    psfcube["emax"].real(100.0)

    psfcube["enumbins"].integer(20)

    psfcube["nxpix"].integer(nxpix)

    psfcube["nypix"].integer(nypix)

    psfcube["binsz"].real(binsz)

    psfcube["coordsys"].string("GAL")

    psfcube["proj"].string("CAR")

    psfcube["xref"].real(0.)

    psfcube["yref"].real(0.)

    psfcube["amax"].real(0.3)

    psfcube["anumbins"].integer(10)

    psfcube.execute()

    print("Done")

    # Run ctbkgcube (builds background cube)

    print("\nBuild background cube")

    bkgcube = ctools.ctbkgcube(sim.obs())

    bkgcube.logFileOpen()

    bkgcube["infile"].filename(events_name)

    bkgcube["cntmap"].filename(cntmap_name)

    bkgcube["bkgmdl"].filename(xml_model_name)

    bkgcube["outfile"].filename(bkgcube_name)

    bkgcube.execute()

    print("Done")

    # Create XML file for stacked analysis. We need this as so far there is

    # no way to directly feed the stacked response files to ctlike and

    # ctmodel. This should be changed in the future, but for the moment

    # that's how to do the job ...

    obs      = gammalib.GObservations()

    cta      = gammalib.GCTAObservation()

    exposure = gammalib.GCTAExposure(expcube_name)

    psf      = gammalib.GCTAMeanPsf(psfcube_name)

    cta.load(cntmap_name)

    cta.response(exposure, psf)

    obs.append(cta)

    obs.save(xml_obs_name)

    # Create model file for stacked analysis. To do this I remove the

    # original background model from the model container and add a diffuse

    # cube background with a constant fitting parameter.

    models = gammalib.GModels(xml_model_name)

    models.remove("Background")

    bgd_spatial  = gammalib.GModelSpatialDiffuseCube(bkgcube_name)

    bgd_spectral = gammalib.GModelSpectralConst(1.0)

    bgd          = gammalib.GCTAModelCubeBackground(bgd_spatial, bgd_spectral)

    models.append(bgd)

    models.save(xml_smodel_name)

    # Run ctlike

    print("\nMaximum likelihood fitting")

    like = ctools.ctlike()

    like.logFileOpen()

    like["infile"].filename(xml_obs_name)

    #like["irf"].string(irf)        # Not needed, response in XML file

    #like["caldb"].string(caldb)    # Not needed, response in XML file

    like["srcmdl"].filename(xml_smodel_name)

    #like["tscalc"].boolean(False)  # This parameter has been removed from ctlike

    like["debug"].boolean(True)

    like["outmdl"].filename(likeout_model)

    like.execute()

    print("Done")

    # Make model map

    print("\nMake model map")

    model = ctools.ctmodel()

    model.logFileOpen()

    model["infile"].filename(cntmap_name)

    #model["obsfile"].filename(events_name) # Use the unstacked response

    model["obsfile"].filename(xml_obs_name) # Use the stacked response (faster!)

    model["outfile"].filename(model_name)

    #model["caldb"].string(caldb)           # Not needed, response in XML file

    #model["irf"].string(irf)               # Not needed, response in XML file

    model["srcmdl"].filename(likeout_model)

    model.execute()

    print("Done")