GammaLib

#! /usr/bin/env python 

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

# This script simulates CTA observations of the Crab nebula and assesses

# to which extend a nebula extension can be detected. Several spatial and

# spectral models are tested

# 

# Author: Rolf Buehler (rolf.buehler@desy.de)

# Date: 20.3.2014

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

import gammalib as gl

from ctools import obsutils

import time

import os

import shutil

import ctools as ct

#import spectrum

import numpy as np

import json

import dicdata

def printandlog(text,tlastlog,logfile,spatial_type,toscreen=True):

    """Helper function to print to screen and log

    @param  [char]     Text text to print

    @param  [time]     Time last log entry was written

    @param  [file]     Log file to write to

    @param  [char]     Spatial model

    @param  [toscreen] Print also to screen

    @return [time]     Time log entry was written out

    """

    currenttime= time.time()

    fulltext=spatial_type+": "+text+"\n"+spatial_type+": (elapsed: "+str(currenttime-tlastlog)+" s)\n"

    if toscreen:

        print fulltext

    logfile.write(fulltext)

    return currenttime

def simcrab(obstime=0.5*3600,write=True,spatial_type="ellipdisk",dospec=False):

    """"Simulates and fits crab & background for one run

    @param [float] obstime  Observation duration in seconds

    @param [bool]  write    Write out simulated observations adn fit results

    """

    #-----------------------------------------------------

    #create sub-folder and setup log to write profiling

    #-----------------------------------------------------

    if os.path.isdir(spatial_type):

        print "Deleting directory:",spatial_type

        shutil.rmtree(spatial_type)

    os.mkdir(spatial_type)

    os.chdir(spatial_type)

    logfile = open("log_"+spatial_type+".log","w")

    lastlogtime=printandlog("Starting. Obstime is "+str(obstime)+" s",time.time(),logfile,spatial_type)

    #-----------------------------------------------------

    #Create one observation (run)

    #-----------------------------------------------------

    ra_degree=83.633212         #Crab position from NED

    dec_degree=22.014460 

    radius=8

    offset_degree=-1            #Observing offset in dec from Crab in deg

    pntdir = gl.GSkyDir()

    pntdir.radec_deg(ra_degree, dec_degree+offset_degree) 

    ob = obsutils.set_obs(pntdir,tstart=0.0,duration=obstime,deadc=0.95, \

                  emin=0.1, emax=100.0, rad=radius, irf="irf_file.fits",\

                  caldb="$PWD/caldb/data/cta/aar/bcf/DESY20140105_50h/")

    ob.id("run_1")

    #----------------------------------------------------

    # Create background model

    #----------------------------------------------------

    spec = gl.GModelSpectralPlaw(1,0,gl.GEnergy(1.0,"TeV"))

    filename="$PWD/caldb/data/cta/aar/bcf/DESY20140105_50h/irf_file.fits"

    bck = gl.GCTAModelIrfBackground(spec)

    bck.name("Background_run_1")

    bck.instruments("CTA")

    #----------------------------------------------------

    # Create Crab nebula model

    # Extension of the Crab in optical is approximatelly an

    # ellipse of: 0.12deg*0.08deg

    #----------------------------------------------------

    fixallpars=False

    crabdir = gl.GSkyDir()

    crabdir.radec_deg(ra_degree, dec_degree)

    if spatial_type=="point":

        crab_spatial = gl.GModelSpatialPointSource(crabdir)

    elif spatial_type=="ellipdisk":

        crab_spatial = gl.GModelSpatialEllipticalDisk(crabdir,0.025,0.015,90.)

    elif spatial_type=="raddisk":

        crab_spatial = gl.GModelSpatialRadialDisk(crabdir,0.02)

    elif spatial_type=="gauss":

        crab_spatial = gl.GModelSpatialRadialGauss(crabdir, 0.02)

    if fixallpars:

        for par in range(0,crab_spatial.size()):

            crab_spatial[par].fix()

    crab_spectrum = gl.GModelSpectralPlaw(3.45e-17, -2.63,gl.GEnergy(1,"TeV") ) #from HESS 2006

    #crab_spectrum = gl.GModelSpectralPlaw2(2.26e-11,-2.63,gl.GEnergy(1,"TeV"),gl.GEnergy(100,"TeV"))

    #crab_spectrum = gl.GModelSpectralLogParabola(3.2e-17,-2.47,gl.GEnergy(1,"TeV"),-0.24) # From MAGIC 2014c

    crab = gl.GModelSky(crab_spatial, crab_spectrum)

    crab.name("CrabNebula")

    crab.ids("")

    #----------------------------------------------------

    # Create GObservations and add models

    #----------------------------------------------------

    obs = gl.GObservations()

    obs.append(ob)

    obs.models().append(bck)

    obs.models().append(crab)

    #----------------------------------------------------

    # Simulate, fit and write out

    #----------------------------------------------------

    obssim = obsutils.sim(obs, log=write, debug=False, edisp=False, seed=0)

    lastlogtime=printandlog("Simulated..",lastlogtime,logfile,spatial_type)

    if write:

        degperpix=0.05

        npix=200

        obssim.save("crabextsim.xml")

        obssim.models().save("crabextsim_models.xml")

        obsutils.cntmap(obssim, proj="TAN", coord="EQU", xval=ra_degree, yval=dec_degree+offset_degree, \

                    binsz=degperpix, nxpix=npix, nypix=npix, outname="crabextsim_cntmap.fits")

        lastlogtime=printandlog("Created counts map..",lastlogtime,logfile,spatial_type)

        obsutils.modmap(obssim, eref=0.1, proj="TAN", coord="EQU", xval=ra_degree, yval=dec_degree+offset_degree, \

                    binsz=degperpix, nxpix=npix, nypix=npix, outname="crabextsim_modmap.fits")

        lastlogtime=printandlog("Created model map..",lastlogtime,logfile,spatial_type)

        for ob in obssim:

            events = ob.events()

            events.save("crabextsim_events_run_1.fits")

    like=obsutils.fit(obssim,log=True, debug=False, edisp=False,tscalc=True)

    lastlogtime=printandlog("Fitted..",lastlogtime,logfile,spatial_type)

    print "TS is",like.obs().models()["CrabNebula"].ts()

    if dospec:

        spec = obsutils.specpoints(obssim,"CrabNebula",obsutils.getlogbins())

        json.dump(spec,open("spectrum.json","w"),indent=4)

        lastlogtime=printandlog("Done Spectrum..",lastlogtime,logfile,spatial_type)

    #----------------------------------------------------

    # Close file and change back directory

    #----------------------------------------------------

    logfile.close()

    os.chdir("..")

if __name__ == '__main__':

    #----------------------------------------------------

    # Run either in multiprocesing or single core

    #----------------------------------------------------

    simcrab()