ctools

#! /usr/bin/env python

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

#

# Copyright (C) 2015-2016 Juergen Knoedlseder

#

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

#

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

# GNU General Public License for more details.

#

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#

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

# This script simulates the CTA Galactic Plane Scan KSP. It requires on

# input the following files:

#

#   gps.xml                    Observation definition file for GPS KSP

#   model/

#         map_ics.fits         Inverse Compton scattering diffuse map

#         map_pi0.fits         Pion decay diffuse map

#         map_RXJ1713.fits     RXJ1713 supernova remnant (X-ray template)

#         map_VelaJunior.fits  Vela Junior supernova remnant

#         model_bkg.xml        Background model definition

#         model_cygnus.xml     Cygnus region model definition

#         model_ics.xml        Inverse Compton scattering model definition

#         model_pi0.xml        Pion decay model definition

#         model_pwne.xml       PWNe population model definition (-TevCAT)

#         model_snrs.xml       SNRs population model definition

#         model_tevcat.xml     TeVCat model definition (thanks, Deidre)

#

# It will produce a file named gps_residual_1-10TeV.fits that contains the

# residual counts above background.

#

# Usage:

# ./make_gps.py

#

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

import gammalib

import ctools

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

# Create response files #

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

def create_response(obsname, model, \

                    emin=0.1, emax=1.0, enumbins=1, \

                    debug=True, force=False, \

                    expcube="expcube.fits", psfcube="psfcube.fits", \

                    bkgcube="bkgcube.fits", bkgmod="model_bkg.xml"):

    """

    Create response files for stacked analysis.

    Parameters:

     obsname  - Filename of observation definition XML file

     model    - Filename of model XML file

    Keywords:

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

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

     enumbins - Number of energy bins in cube (default: 20)

     debug    - Enable debugging (default: True)

     force    - Force computation (default: False)

     expcube  - Exposure cube file name

     psfcube  - PSF cube file name

     bkgcube  - Background cube file name

     bkgmod   - Model file name

    """

    # Create exposure cube

    if force or not gammalib.file_exists(expcube):

        ctexpcube = ctools.ctexpcube()

        ctexpcube["inobs"]    = obsname

        ctexpcube["incube"]   = "NONE"

        ctexpcube["outcube"]  = expcube

        ctexpcube["debug"]    = debug

        ctexpcube["ebinalg"]  = "LOG"

        ctexpcube["emin"]     = emin

        ctexpcube["emax"]     = emax

        ctexpcube["enumbins"] = enumbins

        ctexpcube["nxpix"]    = 3600

        ctexpcube["nypix"]    = 200

        ctexpcube["binsz"]    = 0.1

        ctexpcube["coordsys"] = "GAL"

        ctexpcube["proj"]     = "CAR"

        ctexpcube["xref"]     = 0.0

        ctexpcube["yref"]     = 0.0

        ctexpcube.execute()

    # Create Psf cube

    if force or not gammalib.file_exists(psfcube):

        ctpsfcube = ctools.ctpsfcube()

        ctpsfcube["inobs"]    = obsname

        ctpsfcube["incube"]   = "NONE"

        ctpsfcube["outcube"]  = psfcube

        ctpsfcube["ebinalg"]  = "LOG"

        ctpsfcube["emin"]     = emin

        ctpsfcube["emax"]     = emax

        ctpsfcube["enumbins"] = enumbins

        ctpsfcube["nxpix"]    = 360

        ctpsfcube["nypix"]    = 20

        ctpsfcube["binsz"]    = 1.0

        ctpsfcube["coordsys"] = "GAL"

        ctpsfcube["proj"]     = "CAR"

        ctpsfcube["xref"]     = 0.0

        ctpsfcube["yref"]     = 0.0

        ctpsfcube["debug"]    = debug

        ctpsfcube.execute()

    # Create Background cube

    if force or not gammalib.file_exists(bkgcube):

        ctbkgcube = ctools.ctbkgcube()

        ctbkgcube["inobs"]    = obsname

        ctbkgcube["inmodel"]  = model

        ctbkgcube["incube"]   = "NONE"

        ctbkgcube["outcube"]  = bkgcube

        ctbkgcube["outmodel"] = bkgmod

        ctbkgcube["ebinalg"]  = "LOG"

        ctbkgcube["emin"]     = emin

        ctbkgcube["emax"]     = emax

        ctbkgcube["enumbins"] = enumbins

        ctbkgcube["nxpix"]    = 7200

        ctbkgcube["nypix"]    = 400

        ctbkgcube["binsz"]    = 0.05

        ctbkgcube["coordsys"] = "GAL"

        ctbkgcube["proj"]     = "CAR"

        ctbkgcube["xref"]     = 0.0

        ctbkgcube["yref"]     = 0.0

        ctbkgcube["debug"]    = debug

        ctbkgcube.execute()

    # Return

    return

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

# Create model #

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

def create_model(expcube, psfcube, bkgcube, model, \

                 emin=0.1, emax=1.0, enumbins=1, \

                 nxpix=7200, nypix=400, binsz=0.05, \

                 modcube="modcube.fits", debug=True, force=False):

    """

    Create one specific model for stacked analysis.

    Parameters:

     expcube  - Exposure cube file name

     psfcube  - PSF cube file name

     bkgcube  - Background cube file name

     model    - Filename of model XML file

    Keywords:

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

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

     enumbins - Number of energy bins in cube (default: 20)

     nxpix    - Number of pixels in Galactic longitude (default: 7200)

     nypix    - Number of pixels in Galactic latitude (default: 400)

     binsz    - Pixel size (default: 0.05 deg)

     modcube  - Model cube file name

     debug    - Enable debugging (default: True)

     force    - Force computation (default: False)

    """

    # Create model map

    if force or not gammalib.file_exists(modcube):

        ctmodel = ctools.ctmodel()

        ctmodel["incube"]   = "NONE"

        ctmodel["inmodel"]  = model

        ctmodel["outcube"]  = modcube

        ctmodel["inobs"]    = "NONE"

        ctmodel["expcube"]  = expcube

        ctmodel["psfcube"]  = psfcube

        ctmodel["bkgcube"]  = bkgcube

        ctmodel["ebinalg"]  = "LOG"

        ctmodel["emin"]     = emin

        ctmodel["emax"]     = emax

        ctmodel["enumbins"] = enumbins

        ctmodel["nxpix"]    = nxpix

        ctmodel["nypix"]    = nypix

        ctmodel["binsz"]    = binsz

        ctmodel["coordsys"] = "GAL"

        ctmodel["proj"]     = "CAR"

        ctmodel["xref"]     = 0.0

        ctmodel["yref"]     = 0.0

        ctmodel["debug"]    = debug

        ctmodel.execute()

    # Return

    return

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

# Add and flatten models #

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

def add_and_flatten(models, modmap, name, estr):

    """

    Add and flatten models.

    Parameters:

     models   - List of models

     name     - Name part of file name

     estr     - Energy string part of file name

     modmap   - Flattened model file name

    """

    # Initialise first map flag

    first = True

    # Loop over all models

    for model in models:

        # Build model filename

        filename = "%s_%s_%s.fits" % (name, model, estr)

        # Load model

        model = gammalib.GCTAEventCube(filename)

        # Get sky map and flatten

        map = model.map()

        map.stack_maps()

        # If this is the first map then store it

        if first:

            outmap = map.copy()

            first  = False

        else:

            outmap += map

    # Save map

    outmap.save(modmap, True)

    # Return

    return

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

# Simulate events #

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

def simulate_events(modmap, evtmap):

    """

    Simulate events from model map.

    Parameters:

     modmap   - Flattened model file name

     evtmap   - Event map file name

    """

    # Allocate random number generator

    ran = gammalib.GRan()

    # Load model map

    map = gammalib.GSkyMap(modmap)

    # Randomize

    for i in range(map.npix()):

        map[i] = ran.poisson(map[i])

    # Save event map

    map.save(evtmap, True)

    # Return

    return

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

# Create residual map #

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

def create_residual(evtmap, bkgcomp, resmap):

    """

    Create residual map.

    Parameters:

     evtmap   - Event map file name

     bkgcomp  - Background model file name

     resmap   - Residual map file name

    """

    # Create flat skymap from background component

    bkg    = gammalib.GCTAEventCube(bkgcomp)

    bkgmap = bkg.map()

    bkgmap.stack_maps()

    # Load event map

    events = gammalib.GSkyMap(evtmap)

    # Subtract background model

    for i in range(events.npix()):

        events[i] -= bkgmap[i]

    # Save residual map

    events.save(resmap, True)

    # Return

    return

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

# Main routine entry point #

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

if __name__ == '__main__':

    """

    """

    # Dump header

    print("************************************************")

    print("*    CTA Galactic Plane Scan KSP simulation    *")

    print("************************************************")

    # Set simulation parameters

    name     = "gps"

    emin     =  1.0

    emax     = 10.0

    enumbins = 1

    # Set model components that should be processed

    components = ["model_pi0", \

                  "model_ics", \

                  "model_pwne", \

                  "model_snrs", \

                  "model_tevcat", \

                  "model_cygnus"]

    # Set model components that should be co-added

    models = []

    models.append("model_pi0")

    models.append("model_ics")

    models.append("model_pwne")

    models.append("model_snrs")

    models.append("model_tevcat")

    models.append("model_cygnus")

    models.append("model_bkg")

    # Build names

    if emin < 1.0:

        estr = "%dGeV-%dTeV" % (emin*1.0e3, emax)

    else:

        estr = "%d-%dTeV" % (emin, emax)

    obsname = "%s.xml" % (name)

    expcube = "%s_expcube_%s.fits" % (name, estr)

    psfcube = "%s_psfcube_%s.fits" % (name, estr)

    bkgcube = "%s_bkgcube_%s.fits" % (name, estr)

    bkgmod  = "%s_bkgmodel_%s.xml" % (name, estr)

    bkgcomp = "%s_model_bkg_%s.fits" % (name, estr)

    evtmap  = "%s_events_%s.fits" % (name, estr)

    modmap  = "%s_model_%s.fits" % (name, estr)

    resmap  = "%s_residual_%s.fits" % (name, estr)

    dmodmap = "%s_diff_model_%s.fits" % (name, estr)

    devtmap = "%s_diff_events_%s.fits" % (name, estr)

    dresmap = "%s_diff_residual_%s.fits" % (name, estr)

    # Create response components (exposure, psf, background)

    create_response(obsname, "models/model_bkg.xml", \

                    emin=emin, emax=emax, enumbins=enumbins, \

                    expcube=expcube, psfcube=psfcube, \

                    bkgcube=bkgcube, bkgmod=bkgmod)

    # Create background model component

    create_model(expcube, psfcube, bkgcube, bkgmod, modcube=bkgcomp, \

                 emin=emin, emax=emax, enumbins=enumbins)

    # Create source model components

    for component in components:

        modname = "models/%s.xml" % (component)

        modcomp = "%s_%s_%s.fits" % (name, component, estr)

        create_model(expcube, psfcube, bkgcube, modname, modcube=modcomp, \

                     emin=emin, emax=emax, enumbins=enumbins)

    # Add components and flatten

    add_and_flatten(models, modmap, name, estr)

    # Simulate events

    simulate_events(modmap, evtmap)

    # Create residual map

    create_residual(evtmap, bkgcomp, resmap)