GammaLib

#! /usr/bin/env python

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

# This script simulates CTA observations of Arp220

# 

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

# Date: 20.3.2014

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

import gammalib as gl

import obsutils as obsu

def simarp220(obstime=100.*3600.,write=False):

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

    @param [float] obstime  Observation duration in seconds

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

    """

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

    #Create one observation (run)

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

    ra_degree=233.7382         #Arp220 position from NED

    dec_degree=23.503 

    offset_degree=-0.5            #Observing offset in dec from Arp220 in deg

    pntdir = gl.GSkyDir()

    pntdir.radec_deg(ra_degree, dec_degree+offset_degree)

    # Calibration database and response

    caldb = "aar"

    irf   = "DESY20140105_50h"

    #caldb = "i"

    #irf   = "IFAE20120510_50h"

    ob = obsu.set(pntdir,tstart=0.0,duration=obstime,deadc=0.95, \

                  emin=0.1, emax=10.0, rad=3.0, irf=irf, caldb=caldb)

    ob.id("run_1")

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

    # Create background model

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

    fixallpars=False

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

    bck  = gl.GCTAModelInstBackground(spec)

    bck.name("Background")

    bck.instruments("CTA")

    if fixallpars :

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

            bck[par].fix()

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

    # Create Arp220 model

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

    arp220dir = gl.GSkyDir()

    arp220dir.radec_deg(ra_degree, dec_degree)

    arp220_spatial = gl.GModelSpatialPointSource(arp220dir)

    arp220_spectrum = gl.GModelSpectralPlaw(1.0e-19, -2.2,gl.GEnergy(1,"TeV") )

    arp220 = gl.GModelSky(arp220_spatial, arp220_spectrum)

    arp220.name("Arp220")

    arp220.ids("")

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

    # Create GObservations and add models

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

    obs = gl.GObservations()

    obs.append(ob)

    obs.models().append(bck)

    obs.models().append(arp220)

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

    # Simulate, fit and write out

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

    obssim = obsu.sim(obs, log=write, debug=False, edisp=False, seed=4, nbins=40, binsz=0.02)

    if write:

        obssim.save("arp220_obs.xml")

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

        print "Creating counts map.."

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

                    binsz=0.02, nxpix=300, nypix=300, outname="arp220_cntmap.fits")

        print "Creating model map.."

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

                    binsz=0.02, nxpix=300, nypix=300, outname="arp220_modmap.fits")

        for ob in obssim:

            events = ob.events()

            events.save("arp220_events_run_1.fits")

    print "Fitting.."

    obsu.fit(obssim,log=write, debug=True, edisp=False)

if __name__ == '__main__':

    simarp220()