Merge branch 'devel' of https://cta-gitlab.irap.omp.eu/ctools/ctools into 2207-c…

…tskymap_RING_bkgsubtract

Merge branch 'devel' of https://cta-gitlab.irap.omp.eu/ctools/ctools into 2207-c…
…tskymap_RING_bkgsubtract
Josh Cardenzana
2 parents 9a32c2d4 5a2e99f2
Showing 20 changed files with 1071 additions and 89 deletions
ChangeLog
NEWS
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csmodelsois.par
cscripts/csmodelsois.py
doc/source/admin/release_history/1.5.rst
doc/source/users/reference_manual/csmodelsois.rst
doc/source/users/reference_manual/ctmapcube.rst
examples/show_spectrum.py
pyext/ctobservation.i
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
src/ctskymap/ctskymap.cpp
src/ctskymap/ctskymap.hpp
src/support/ctobservation.cpp
src/support/ctobservation.hpp
test/test_csmodelsois.py
test/test_ctskymap.py
test/test_python_cscripts.py
-2017-08-08
+2017-10-11
         * Version 1.5.0 released
         ========================
+        Add keywords to sky maps produced by ctskymap (#2194)
+        Enhance efficiency of computations in ctmodel (#2189)
+        Add csmodelsois script (#2180)
         Use t-type parameters for all "tmin" and "tmax" parameters (#1864)
 New Features and Important Changes in ctools 1.5.0
-8 August 2017
+11 October 2017
 Introduction
@@ -53,7 +53,9 @@ None
 ctskymap - CTA sky mapping
 --------------------------
-None
+Header keywords informing about the dates and times covered by the observations
+used to generate the skymap, as well as the event energy range, are now appended
+as header keywords to the FITS file (#2194).
 ctmodel - CTA counts cube modelling
@@ -61,6 +63,11 @@ ctmodel - CTA counts cube modelling
 Use t-type parameters for "tmin" and "tmax" to all specifying times in various
 formats (#1864).
+The tool now also checks which models actually overlap with the RoI of a given
+observation and then only use these models to compute the model cube. The
+parallelization of the tool has also been optimised so that diffuse models are
+not loaded over and over again (#2189).
+
 ctexpcube - CTA exposure cube generation
 ----------------------------------------
@@ -119,6 +126,9 @@ None
 cscripts
 --------
+Add csmodelsois script that replaces selected models in a model container by
+a diffuse map cube (#2180).
+
 Use t-type parameters for "tmin" and "tmax" to all specifying times in various
 formats (#1864).
@@ -28,6 +28,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
                 $(srcdir)/csmodelinfo.py \
                 $(srcdir)/csmodelmerge.py \
                 $(srcdir)/csmodelselect.py \
+                $(srcdir)/csmodelsois.py \
                 $(srcdir)/csobs2caldb.py \
                 $(srcdir)/csobsdef.py \
                 $(srcdir)/csobsinfo.py \
@@ -60,6 +61,7 @@ parfiles = $(srcdir)/cscaldb.par \
            $(srcdir)/csmodelinfo.par \
            $(srcdir)/csmodelmerge.par \
            $(srcdir)/csmodelselect.par \
+           $(srcdir)/csmodelsois.par \
            $(srcdir)/csobs2caldb.par \
            $(srcdir)/csobsdef.par \
            $(srcdir)/csobsinfo.par \
@@ -89,6 +91,7 @@ linkscripts = cscaldb \
               csmodelinfo \
               csmodelmerge \
               csmodelselect \
+              csmodelsois \
               csobs2caldb \
               csobsdef \
               csobsinfo \
@@ -123,6 +126,7 @@ testscripts = $(srcdir)/tests/__init__.py \
               $(top_srcdir)/test/test_csmodelinfo.py \
               $(top_srcdir)/test/test_csmodelmerge.py \
               $(top_srcdir)/test/test_csmodelselect.py \
+              $(top_srcdir)/test/test_csmodelsois.py \
               $(top_srcdir)/test/test_csobs2caldb.py \
               $(top_srcdir)/test/test_csobsdef.py \
               $(top_srcdir)/test/test_csobsinfo.py \
 # ==========================================================================
 # cscripts Python module
 #
-# Copyright (C) 2015-2016 Juergen Knoedlseder
+# Copyright (C) 2015-2017 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
@@ -31,6 +31,7 @@ __all__ = [
     "csmodelinfo",
     "csmodelmerge",
     "csmodelselect",
+    "csmodelsois",
     "csobs2caldb",
     "csobsdef",
     "csobsinfo",
@@ -66,6 +67,7 @@ from cscripts.cslightcrv    import cslightcrv
 from cscripts.csmodelinfo   import csmodelinfo
 from cscripts.csmodelmerge  import csmodelmerge
 from cscripts.csmodelselect import csmodelselect
+from cscripts.csmodelsois   import csmodelsois
 from cscripts.csobs2caldb   import csobs2caldb
 from cscripts.csobsdef      import csobsdef
 from cscripts.csobsinfo     import csobsinfo
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      csmodelsois.par                                               #
+#                                                                           #
+#  Author:    Josh Cardenzana                                               #
+#             josh.cardenzana@irap.omp.eu                                   #
+#             IRAP                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the csmodelsois application                #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+inmodel,  f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
+outcube,  f, a, mapcube.fits,,, "Output map cube file"
+ptsrcsig, r, h, 1.0,0,60, "Point source Gaussian sigma in arcmin"
+soilist,  s, a, ,,, "CSV list of sources to exclude from cube model"
+outmodel, f, a, model.xml,,, "Output model definition XML file"
+
+#
+# Energy binning parameters
+#==========================
+ebinalg,  s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins"
+emin,     r, a, 0.1,,, "Start value for first energy bin in TeV"
+emax,     r, a, 100.0,,, "Stop value for last energy bin in TeV"
+enumbins, i, a, 20,1,200, "Number of energy bins"
+ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
+
+#
+# Spatial binning parameters
+#===========================
+nxpix,    i, a, 200,,, "Size of the X axis in pixels"
+nypix,    i, a, 200,,, "Size of the Y axis in pixels"
+binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
+coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
+ra,       r, a, 83.63,0,360, "Right ascension of image center (in degrees)"
+dec,      r, a, 22.01,-90,90, "Declination of image center (in degrees)"
+glon,     r, a, 184.55,0,360, "Galactic longitude of image center (in degrees)"
+glat,     r, a, -5.80,-90,90, "Galactic latitude of image center (in degrees)"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish map cube on VO Hub?"
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,   b, h, no,,, "Debugging mode activated"
+mode,    s, h, "ql",,, "Mode of automatic parameters"
+logfile, f, h, "csmodelsois.log",,, "Log filename"
+#! /usr/bin/env python
+# ==========================================================================
+# Puts subset of sources in diffuse model cube
+#
+# Copyright (C) 2017 Josh Cardenzana
+#
+# 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/>.
+#
+# ==========================================================================
+import sys
+import gammalib
+import ctools
+
+
+# ================= #
+# csmodelsois class #
+# ================= #
+class csmodelsois(ctools.cscript):
+    """
+    Puts subset of sources in diffuse model cube
+
+    The csmodelsois class puts a subset of sources in a model definition
+    XML file into a model cube and generates a new model definition XML file
+    in which the subset of sources is replaced by a model of type
+    DiffuseMapCube.
+    """
+
+    # Constructor
+    def __init__(self, *argv):
+        """
+        Constructor
+        
+        Parameters
+        ----------
+        argv : list of parameters
+        
+        Raises
+        ------
+        TypeError
+            An invalid number of command line arguments was provided.
+        """
+        # Set name and version
+        self._name    = 'csmodelsois'
+        self._version = ctools.__version__
+        
+        # Initialize parameters
+        self._init_members()
+
+        # Initialise application by calling the appropriate class constructor
+        self._init_cscript(argv)
+        
+        # Return
+        return
+
+
+    # Private methods
+    def _init_members(self):
+        """
+        Initialize the parameters for this object
+        """
+        # Initialise parameters
+        self.cubegen    = ctools.ctmapcube() # Implements ctmapcube for generate the cube
+        self.m_models   = gammalib.GModels() # Stores the full list of input models
+        self.cubemodels = gammalib.GModels() # Stores the list of models to be written to file
+
+        # Return
+        return
+
+    def _get_parameters(self):
+        """
+        Get parameters from parfile
+        """
+        # Query input parameters
+        self['inmodel'].filename()
+        self['soilist'].string()
+
+        # Read input parameters
+        ptsrcsig = self['ptsrcsig'].real()
+        
+        # Get the cube centre coordinates
+        self.cubegen['coordsys'].string(self['coordsys'].string())
+        if self['coordsys'].string() == 'CEL':
+            self.cubegen['xref'].real(self['ra'].real())
+            self.cubegen['yref'].real(self['dec'].real())
+        else:
+            self.cubegen['xref'].real(self['glon'].real())
+            self.cubegen['yref'].real(self['glat'].real())
+
+        # Get spatial binning parameters
+        self.cubegen['binsz'].real(self['binsz'].real())
+        self.cubegen['nxpix'].integer(self['nxpix'].integer())
+        self.cubegen['nypix'].integer(self['nypix'].integer())
+        self.cubegen['proj'].string(self['proj'].string())
+        
+        # Get the energy binning parameters
+        self.cubegen['ebinalg'].string(self['ebinalg'].string())
+        if self['ebinalg'].string() == 'FILE':
+            self.cubegen['ebinfile'].filename(self['ebinfile'].filename())
+        else:
+            self.cubegen['emin'].real(self['emin'].real())
+            self.cubegen['emax'].real(self['emax'].real())
+            self.cubegen['enumbins'].integer(self['enumbins'].integer())
+
+        # Get remaining hidden parameters
+        self.cubegen['ptsrcsig'].real(ptsrcsig)
+        
+        # Read optionally output cube filenames
+        if self._read_ahead():
+            self.cubegen['outcube'].filename(self['outcube'].filename())
+            self['outmodel'].filename()
+
+        # Write parameters into logger
+        self._log_parameters(gammalib.TERSE);
+
+        return
+
+    def _gen_cubemodel(self):
+        """
+        Generates a binned model from the model file using ctmapcube
+        """
+        # Get a list of sources we DONT want to put into the cube
+        sources = self['soilist'].string().split(',')
+
+        # Load the model into a model container
+        if (self.m_models.size() == 0):
+            self.m_models = gammalib.GModels(self['inmodel'].filename())
+
+        # Store a copy of the models
+        self.cubemodels = self.m_models.copy()
+
+        # Loop through models and pull out all models not in the list
+        models2remove = []
+        for model in self.m_models:
+
+            # If model is not a sky model then continue
+            if model.classname() != 'GModelSky':
+                self.cubemodels.remove(model.name())
+
+            # ... otherwise, if model is in list of sources the remove it
+            # from the container
+            elif model.name() in sources:
+                self.cubemodels.remove(model.name())
+
+        # Log the number of models to be put into the cube
+        self._log_value(self['chatter'].integer(), 'Numner of cube models',
+                        self.cubemodels.size())
+
+        # Return
+        return
+
+
+    # Public methods
+    def run(self):
+        """
+        Implements the actual bulk of the script's tasks
+        """
+        # Switch screen logging on in debug mode
+        if self._logDebug():
+            self._log.cout(True)
+
+        # Get parameters
+        self._get_parameters()
+
+        # Set the models to put into the cube
+        self._gen_cubemodel()
+
+        # Only run if there still models to be filled into the cube
+        if not self.cubemodels.is_empty():
+
+            # Log action
+            self._log_string(gammalib.NORMAL, 'Generating model cube')
+
+            # Fill the cube models into the ctmapcube object
+            self.cubegen.models(self.cubemodels)
+
+            # Run the map generating method
+            self.cubegen.run()
+
+        else:
+            self._log_string(gammalib.NORMAL,
+                             'List of models is empty, nothing will be done')
+
+        # Return
+        return
+
+    def execute(self):
+        """
+        Execute the script
+        """
+        # Open logfile
+        self.logFileOpen()
+
+        # Read ahead output parameters
+        self._read_ahead(True)
+
+        # Run the script
+        self.run()
+
+        # Save model
+        self.save()
+
+        # Return
+        return
+
+    def save(self):
+        """
+        Save the cube and optionally the updated XML model file
+        """
+        # Save the generated cube if the cube and filename are not empty
+        if ((not self.cubegen.mapcube().cube().is_empty()) and
+            self._is_valid_filename(self['outcube'].filename())):
+            self.cubegen.mapcube().save(self['outcube'].filename(),
+                                        self['clobber'].boolean())
+        
+        # If requested, save the updated list of models
+        if self._is_valid_filename(self['outmodel'].filename()):
+
+            # Generate a list of models that will be output to file
+            outmodels = gammalib.GModels(self.m_models)
+            
+            # Remove all models used in the generated cube
+            for model in self.cubemodels:
+                outmodels.remove(model.name())
+
+            # Generate the actual cube model
+            spat     = gammalib.GModelSpatialDiffuseCube(self['outcube'].filename())
+            spec     = gammalib.GModelSpectralConst()
+            newmodel = gammalib.GModelSky(spat, spec)
+
+            # Give the model a default name
+            newmodel.name(self.cubemodelname())
+
+            # Now append the model to the list of output models
+            outmodels.append(newmodel)
+
+            # Save the model list
+            outmodels.save(self['outmodel'].filename())
+
+        # Return
+        return
+
+    def mapcube(self):
+        """
+        Return the mapcube generated by the underlying 'ctmapcube' object
+        """
+        # Return
+        return self.cubegen.mapcube()
+
+    def cubemodelname(self):
+        """
+        Return the name to be given to the cube model in the output XML file
+        """
+        # Return
+        return 'FixedSourcesCube'
+
+    def models(self, models):
+        """
+        Set the models
+
+        Parameters
+        -------
+        models : `~gammalib.GModels`
+            Model container
+        """
+        # Set models
+        self.m_models = models
+
+        # Return
+        return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    
+    # Create instance of application
+    app = csmodelsois(sys.argv)
+    
+    # Execute application
+    app.execute()
@@ -22,11 +22,16 @@ Bug fixes
 Improvements
 ------------
-* TBW
+* [`2194 <https://cta-redmine.irap.omp.eu/issues/2194>`_] -
+  Add keywords to sky maps produced by :ref:`ctskymap`
+* [`2189 <https://cta-redmine.irap.omp.eu/issues/2189>`_] -
+  Enhance efficiency of computations in :ref:`ctmodel`
 New features
 ------------
+* [`2180 <https://cta-redmine.irap.omp.eu/issues/2180>`_] -
+  Add :ref:`csmodelsois` script
 * [`1864 <https://cta-redmine.irap.omp.eu/issues/1864>`_] -
   Use t-type parameters for all ``tmin`` and ``tmax`` parameters
+.. _csmodelsois:
+
+csmodelsois
+===========
+
+Generate a map cube from a subset of models in an input model.
+
+
+Synopsis
+----------
+
+This tool generates a map cube from a subset of sky model components contained
+in an input model. A map cube is a series of sky maps, each spanned by Right
+Ascension or Galactic longitude and Declination or Galactic latitude. Each
+map corresponds to a specific photon energy. The energy binning of the map
+cube may be either linear, logarithmic, or custom defined using an input
+file.
+
+csmodelsois generates a map cube FITS file comprising two extensions. The
+primary extension contains a 3-dimensional image containing the map cube
+values. A second extension named ``ENERGIES`` contains a binary table that
+defines the energy values of each map in the cube.
+
+Optionally, csmodelsois can also output an updated model XML file in which all
+sources from the input model XML file used to generate the map cube are
+replaced with an updated model that uses the generated map cube. This is done
+by specifying the 'outmodel' parameter. The 'soilist' parameter specifies a 
+list of source names (separated by ',') to be excluded from the generated
+cube.
+
+General parameters
+------------------
+
+``inmodel [file]``
+    Input model XML file.
+
+``outmodel [file]``
+    Output model XML file (not generated if set to ``NONE``).
+
+``outcube [file]``
+    Output map cube file.
+
+``soilist [string]``
+    CSV list of sources to be excluded from the generated map cube.
+
+``(ptsrcsig = 1.0) [real]``
+    Sigma of Gaussian to be used for point sources (in arcmin).
+
+``ebinalg <FILE|LIN|LOG> [string]``
+    Algorithm for defining energy bins.
+
+``emin [real]``
+    Lower energy value for first energy bin (in TeV).
+
+``emax [real]``
+    Upper energy value for last energy bin (in TeV).
+
+``enumbins [integer]``
+    Number of energy bins.
+
+``ebinfile [file]``
+    Name of the file containing the energy bin definition.
+
+``nxpix [integer]``
+    Number of cube bins in Right Ascension or Galactic longitude.
+
+``nypix [integer]``
+    Number of cube bins in Declination or Galactic latitude.
+
+``binsz [real]``
+    Image scale (in degrees/pixel).
+
+``coordsys <CEL|GAL> [string]``
+    Coordinate system (CEL - celestial, GAL - galactic).
+
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
+``ra [real]``
+    Right Ascension of cube centre (J2000, in degrees).
+
+``dec [real]``
+    Declination of cube centre (J2000, in degrees).
+
+``glon [real]``
+    Galactic longitude of cube centre (in degrees).
+
+``glat [real]``
+    Galactic latitude of cube centre (in degrees).
+
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+    Specifies whether the map cube should be published on VO Hub.
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     ``chatter = 0``: no information will be logged
+
+     ``chatter = 1``: only errors will be logged
+
+     ``chatter = 2``: errors and actions will be logged
+
+     ``chatter = 3``: report about the task execution
+
+     ``chatter = 4``: detailed report about the task execution
+
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube file should be overwritten.
+
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = csmodelsois.log) [string]``
+    Name of log file.
+
+
+Related tools or scripts
+------------------------
+
+:ref:`ctmapcube`
@@ -36,40 +36,40 @@ General parameters
 ``ebinalg <FILE|LIN|LOG> [string]``
     Algorithm for defining energy bins.
- 	 	 
+
 ``emin [real]``
     Lower energy value for first energy bin (in TeV).
- 	 	 
+
 ``emax [real]``
     Upper energy value for last energy bin (in TeV).
- 	 	 
+
 ``enumbins [integer]``
     Number of energy bins.
- 	 	 
+
 ``ebinfile [file]``
     Name of the file containing the energy bin definition.
- 	 	 
+
 ``nxpix [integer]``
     Number of cube bins in Right Ascension or Galactic longitude.
- 	 	 
+
 ``nypix [integer]``
     Number of cube bins in Declination or Galactic latitude.
- 	 	 
+
 ``binsz [real]``
     Image scale (in degrees/pixel).
- 	 	  	 	 
+
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
- 	 	 
+
 ``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
- 	 	 
+
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
- 	 	 
+
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
- 	 	 
+
 Standard parameters
 -------------------
@@ -80,21 +80,21 @@ Standard parameters
 ``(chatter = 2) [integer]``
     Verbosity of the executable:
      ``chatter = 0``: no information will be logged
-     
+
      ``chatter = 1``: only errors will be logged
-     
+
      ``chatter = 2``: errors and actions will be logged
-     
+
      ``chatter = 3``: report about the task execution
-     
+
      ``chatter = 4``: detailed report about the task execution
- 	 	 
+
 ``(clobber = yes) [boolean]``
     Specifies whether an existing output counts cube file should be overwritten.
- 	 	 
+
 ``(debug = no) [boolean]``
     Enables debug mode. In debug mode the executable will dump any log file output to the console.
- 	 	 
+
 ``(mode = ql) [string]``
     Mode of automatic parameters (default is "ql", i.e. "query and learn").
@@ -105,4 +105,4 @@ Standard parameters
 Related tools or scripts
 ------------------------
-None
+:ref:`csmodelsois`
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Display spectrum generated by csspec
 #
-# Copyright (C) 2015-2016 Juergen Knoedlseder
+# Copyright (C) 2015-2017 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
@@ -91,6 +91,9 @@ def plot_spectrum(filename, plotfile):
             ul_ed_engs.append(c_ed.real(row))
             ul_eu_engs.append(c_eu.real(row))
+    # Set upper limit errors
+    yerr = [0.6 * x for x in ul_flux]
+
     # Plot the spectrum 
     plt.figure()
     plt.loglog()
@@ -98,7 +101,7 @@ def plot_spectrum(filename, plotfile):
     plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs],
                  fmt='ro')
     plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs],
-                 yerr=1.0e-11, uplims=True, fmt='ro')
+                 yerr=yerr, uplims=True, fmt='ro')
     plt.xlabel('Energy (TeV)')
     plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)')
 /***************************************************************************
  *             ctobservation - Base class for observation tools            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2016 by Juergen Knoedlseder                              *
+ *  copyright (C) 2016-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -89,6 +89,7 @@ public:
     // Public methods
     const GObservations& obs(void) const;
+    void                 write_ogip_keywords(GFitsHDU* hdu) const;
     // Make methods private in Python by prepending an underscore
     %rename(_first_unbinned_observation) first_unbinned_observation;
@@ -235,48 +235,55 @@ void ctmodel::run(void)
     log_header1(TERSE, "Generate model cube");
     // Loop over all observations in the container
-    #pragma omp parallel for
-    for (int i = 0; i < m_obs.size(); ++i) {
+    #pragma omp parallel
+    {
+        // Create a models object for this core
+        GModels models(m_obs.models());
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Skip observation if it's not CTA
-        if (obs == NULL) {
-            log_header3(TERSE, get_obs_header(m_obs[i]));
-            std::string msg = " Skipping "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-        }
-
-        // Fill cube and leave loop if we are binned mode (meaning we
-        // only have one binned observation)
-        else if (m_binned) {
-            fill_cube(obs);
-            i = m_obs.size();
-        }
-
-        // Skip observation if we have a binned observation
-        else if (obs->eventtype() == "CountsCube") {
-            log_header3(TERSE, get_obs_header(m_obs[i]));
-            std::string msg = " Skipping binned "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-        }
-
-        // Otherwise, everything seems to be fine, so fill the cube from obs
-        else {
-            // Fill the cube
-            fill_cube(obs);
+        #pragma omp for schedule(static,1)
+        for (int i = 0; i < m_obs.size(); ++i) {
-            // Dispose events to free memory if event file exists on disk
-            if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
-                obs->dispose_events();
+            // Get CTA observation
+            GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+            
+            // Skip observation if it's not CTA
+            if (obs == NULL) {
+                log_header3(TERSE, get_obs_header(m_obs[i]));
+                std::string msg = " Skipping "+m_obs[i]->instrument()+
+                " observation";
+                log_string(NORMAL, msg);
             }
-        }
-
-    } // endfor: looped over observations
+            
+            // Fill cube and leave loop if we are binned mode (meaning we
+            // only have one binned observation)
+            else if (m_binned) {
+                fill_cube(obs, models);
+                i = m_obs.size();
+            }
+            
+            // Skip observation if we have a binned observation
+            else if (obs->eventtype() == "CountsCube") {
+                log_header3(TERSE, get_obs_header(m_obs[i]));
+                std::string msg = " Skipping binned "+m_obs[i]->instrument()+
+                " observation";
+                log_string(NORMAL, msg);
+            }
+            
+            // Otherwise, everything seems to be fine, so fill the cube from obs
+            else {
+                // Fill the cube
+                fill_cube(obs, models);
+                
+                // Dispose events to free memory if event file exists on disk
+                if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
+                    obs->dispose_events();
+                }
+            }
+            
+        } // endfor: looped over observations
+    } // end openmp parallel section
+    
     // Write model cube into header
     log_header1(NORMAL, "Model cube");
     log_string(NORMAL, m_cube.print());
@@ -789,7 +796,7 @@ void ctmodel::extract_cube_properties(void)
  * GTI of the model cube so that cube GTI is a list of the GTIs of all
  * observations that were used to generate the model cube.
  ***************************************************************************/
-void ctmodel::fill_cube(const GCTAObservation* obs)
+void ctmodel::fill_cube(const GCTAObservation* obs, GModels& models)
 {
     // Get number of spatial and spectral bins in counts cube
     int npix   = m_cube.npix();
@@ -819,8 +826,8 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     int    num_outside_roi  = 0;
     // Extract copy of models from observation container
-    GModels models(m_obs.models());
-
+    GModels ex_models = trim_models(models, roi);
+    
     // Get pointer to event cube pixels
     double* pixels = const_cast<double*>(m_cube.counts().pixels());
@@ -877,6 +884,9 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     } // endfor: looped over all cube bins
+    // Re-append the models to the main model
+    models.extend(ex_models);
+    
     // Update GTIs
     #pragma omp critical(ctmodel_fill_cube)
     {
@@ -905,3 +915,66 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     // Return
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Find the models falling inside a defined region of interest.
+ *
+ * @param[in,out] all_models Model container to be trimmed.
+ * @param[in] roi Observation region of interest.
+ * @return New model container containing only models inside RoI.
+ *
+ * Note that a buffer is added to the observation region to ensure that
+ * point sources (which have radius=0) are appropriately included if they
+ * fall near the edge of the observation.
+ ***************************************************************************/
+GModels ctmodel::trim_models(GModels& all_models, const GCTARoi& roi)
+{
+    // Create the model container containing only models in RoI
+    GModels excluded_models;
+
+    // Do nothing if roi is not valid. This will be the case for binned data
+    // sets as they have the default roi radius of 0. In this case we would
+    // remove ALL models, which is obviously incorrect. So we keep all models
+    // when filling based on a counts cube.
+    if (!roi.is_valid()) {
+        return excluded_models;
+    }
+
+    // Remove all models that dont overlap with the region of interest. Note
+    // that an extra factor is used since point sources have regions of
+    // radius 0, which is a problem when they fall just barely outside the
+    // RoI. This ensures point sources are appropriately included.
+    GSkyRegionCircle obsreg(roi.centre().dir(), roi.radius()+0.5);
+
+    // Loop over the models in the passed container
+    for (int i = 0; i < all_models.size(); ++i) {
+
+        // Cast the model to a GModelSky object
+        GModelSky* model = dynamic_cast<GModelSky*>(all_models[i]);
+
+        // If model is not a GModelSky model than skip
+        if (model == NULL) {
+            continue;
+        }
+
+        // Otherwise, exclude the model if it isn't a diffuse cube and it
+        // doesn't overlap with the observation.
+        if (!(model->spatial()->code() == GMODEL_SPATIAL_DIFFUSE) &&
+            (!model->spatial()->region()->overlaps(obsreg))) {
+
+            // Append model to the excluded models
+            excluded_models.append(*model);
+
+            // Remove model
+            all_models.remove(i);
+
+            // Decrement to prevent skipping models
+            i--;
+        }
+
+    } // endfor: looped over all models
+
+    // Return excluded models
+    return excluded_models;
+}
@@ -74,14 +74,15 @@ public:
 protected:
     // Protected methods
-    void init_members(void);
-    void copy_members(const ctmodel& app);
-    void free_members(void);
-    void get_parameters(void);
-    void get_obs(void);
-    void extract_cube_properties(void);
-    void fill_cube(const GCTAObservation* obs);
-    bool has_cube(void) const;
+    void    init_members(void);
+    void    copy_members(const ctmodel& app);
+    void    free_members(void);
+    void    get_parameters(void);
+    void    get_obs(void);
+    void    extract_cube_properties(void);
+    void    fill_cube(const GCTAObservation* obs, GModels& models);
+    bool    has_cube(void) const;
+    GModels trim_models(GModels& all_models, const GCTARoi& roi);
     // User parameters
     GFilename m_outcube;      //!< Output model cube
@@ -281,23 +281,39 @@ void ctskymap::save(void)
         GFits fits;
         // Write sky map into FITS file
-        m_skymap.write(fits);
+        GFitsHDU* hdu = m_skymap.write(fits);
+
+        // Write keywords into sky map extension
+        write_ogip_keywords(hdu);
+        write_hdu_keywords(hdu);
         // If background subtraction is requested then write background map
         // and significance map to FITS file
         if (m_bkgsubtract != "NONE") {
             // Write background map into FITS file
-            m_bkgmap.write(fits);
+            hdu = m_bkgmap.write(fits);
             // Set background map extension name
-            fits[fits.size()-1]->extname("BACKGROUND");
+            if (hdu != NULL) {
+                hdu->extname("BACKGROUND");
+            }
+
+            // Write keywords into background extension
+            write_ogip_keywords(hdu);
+            write_hdu_keywords(hdu);
             // Write significance map into FITS file
-            m_sigmap.write(fits);
+            hdu = m_sigmap.write(fits);
             // Set significance map extension name
-            fits[fits.size()-1]->extname("SIGNIFICANCE");
+            if (hdu != NULL) {
+                hdu->extname("SIGNIFICANCE");
+            }
+
+            // Write keywords into significance extension
+            write_ogip_keywords(hdu);
+            write_hdu_keywords(hdu);
         } // endif: background subtraction was requested
@@ -525,7 +541,7 @@ void ctskymap::map_events(GCTAObservation* obs)
     // Make sure that the observation holds a CTA event list. If this
     // is not the case then throw an exception.
     if (events == NULL) {
-        throw GException::no_list(G_BIN_EVENTS);
+        throw GException::no_list(G_MAP_EVENTS);
     }
     // Setup energy range covered by data
@@ -647,7 +663,7 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
                           get_obs_header(obs)+" to compute IRF background. "
                           "Please specify response information or use "
                           "another background subtraction method.";
-        throw GException::invalid_value(G_BKG_SUBTRACT_IRF, msg);
+        throw GException::invalid_value(G_MAP_BACKGROUND_IRF, msg);
     }
     // Get IRF background template
@@ -659,7 +675,7 @@ void ctskymap::map_background_irf(GCTAObservation* obs)
                           "response function for "+
                           get_obs_header(obs)+". Please specify an instrument "
                           "response function containing a background template.";
-        throw GException::invalid_value(G_BKG_SUBTRACT_IRF, msg);
+        throw GException::invalid_value(G_MAP_BACKGROUND_IRF, msg);
     }
     // Compute natural logarithm of energy range in MeV
@@ -860,6 +876,30 @@ void ctskymap::map_background_ring(GCTAObservation* obs)
     return;
 }
+/***********************************************************************//**
+ * @brief Write keywords in FITS HDU
+ *
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Writes keywords in FITS HDU.
+ ***************************************************************************/
+void ctskymap::write_hdu_keywords(GFitsHDU* hdu) const
+{
+    // Continue only if pointer is valid
+    if (hdu != NULL) {
+
+        // Set keywords
+        hdu->card("BKGSUB", m_bkgsubtract, "Background substraction method");
+        hdu->card("E_MIN",  m_emin, "[TeV] Lower energy boundary");
+        hdu->card("E_MAX",  m_emax, "[TeV] Upper energy boundary");
+        hdu->card("EUNIT",  "TeV",  "Units for E_MIN and E_MAX");
+
+    } // endif: pointer was valid
+
+    // Return
+    return;
+}
+
 /***********************************************************************//**
  * @brief Computes the bin-by-bin significance if background subtraction requested
@@ -80,6 +80,7 @@ protected:
                              const GSkyMap& sensitivity,
                              const GSkyDir& position,
                              double& non, double& noff, double& alpha);
+    void write_hdu_keywords(GFitsHDU* hdu) const;
     // Background integration kernel
     class irf_kern : public GFunction {
@@ -192,6 +192,53 @@ ctobservation&amp; ctobservation::operator=(const ctobservation&amp; app)
  =                                                                         =
  ==========================================================================*/
+/***********************************************************************//**
+ * @brief Write OGIP keywords in FITS HDU
+ *
+ * @param[in,out] hdu Pointer to FITS HDU.
+ *
+ * Writes OGIP keywords in FITS HDU.
+ ***************************************************************************/
+void ctobservation::write_ogip_keywords(GFitsHDU* hdu) const
+{
+    // Continue only if pointer is valid
+    if (hdu != NULL) {
+
+        // Set creator
+        std::string creator = this->name() + " v" + this->version();
+
+        // Compute times
+        std::string utc_obs  = m_ogip_tstart.utc();
+        std::string utc_end  = m_ogip_tstop.utc();
+        std::string date_obs = utc_obs.substr(0, 10);
+        std::string time_obs = utc_obs.substr(11, 8);
+        std::string date_end = utc_end.substr(0, 10);
+        std::string time_end = utc_end.substr(11, 8);
+
+        // Compute deadtime correction
+        double deadc = (m_ogip_ontime > 0.0) ? m_ogip_livetime / m_ogip_ontime : 1.0;
+
+        // Set observation information
+        hdu->card("CREATOR",  creator,          "Program which created the file");
+        hdu->card("TELESCOP", m_ogip_telescope, "Telescope");
+
+        // Set observation time information
+        hdu->card("DATE_OBS", date_obs,        "Observation start date");
+        hdu->card("TIME_OBS", time_obs,        "Observation start time");
+        hdu->card("DATE_END", date_end,        "Observation end date");
+        hdu->card("TIME_END", time_end,        "Observation end time");
+        hdu->card("TELAPSE",  m_ogip_telapse,  "[s] Elapsed time");
+        hdu->card("ONTIME",   m_ogip_ontime,   "[s] Total good time including deadtime");
+        hdu->card("LIVETIME", m_ogip_livetime, "[s] Total livetime");
+        hdu->card("EXPOSURE", m_ogip_exposure, "[s] Exposure time");
+        hdu->card("DEADC",    deadc,           "Deadtime correction factor");
+
+    } // endif: pointer was valid
+
+    // Return
+    return;
+}
+
 /*==========================================================================
  =                                                                         =
@@ -208,6 +255,13 @@ void ctobservation::init_members(void)
     m_obs.clear();
     // Initialise private members
+    m_ogip_telescope.clear();
+    m_ogip_tstart.clear();
+    m_ogip_tstop.clear();
+    m_ogip_telapse   = 0.0;
+    m_ogip_exposure  = 0.0;
+    m_ogip_ontime    = 0.0;
+    m_ogip_livetime  = 0.0;
     m_index_unbinned = 0;
     // Return
@@ -226,6 +280,13 @@ void ctobservation::copy_members(const ctobservation&amp; app)
     m_obs = app.m_obs;
     // Copy private members
+    m_ogip_telescope = app.m_ogip_telescope;
+    m_ogip_tstart    = app.m_ogip_tstart;
+    m_ogip_tstop     = app.m_ogip_tstop;
+    m_ogip_telapse   = app.m_ogip_telapse;
+    m_ogip_exposure  = app.m_ogip_exposure;
+    m_ogip_ontime    = app.m_ogip_ontime;
+    m_ogip_livetime  = app.m_ogip_livetime;
     m_index_unbinned = app.m_index_unbinned;
     // Return
@@ -258,6 +319,15 @@ const GCTAObservation* ctobservation::first_unbinned_observation(void) const
     // Initialise index
     m_index_unbinned = 0;
+    // Initialise OGIP members
+    m_ogip_telescope.clear();
+    m_ogip_tstart.clear();
+    m_ogip_tstop.clear();
+    m_ogip_telapse  = 0.0;
+    m_ogip_exposure = 0.0;
+    m_ogip_ontime   = 0.0;
+    m_ogip_livetime = 0.0;
+
     // Get next unbinned CTA observation
     const GCTAObservation* obs = next_unbinned_observation();
@@ -323,6 +393,31 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const
     } // endfor: looped over all observation
+    // Update OGIP members
+    if (obs != NULL) {
+        if (m_ogip_telapse == 0.0) {
+            m_ogip_telescope = obs->instrument();
+            m_ogip_tstart    = obs->gti().tstart();
+            m_ogip_tstop     = obs->gti().tstop();
+            m_ogip_telapse   = obs->gti().telapse();
+            m_ogip_exposure  = obs->ontime();
+            m_ogip_ontime    = obs->ontime();
+            m_ogip_livetime  = obs->livetime();
+        }
+        else {
+            if (obs->gti().tstart() < m_ogip_tstart) {
+                m_ogip_tstart = obs->gti().tstart();
+            }
+            if (obs->gti().tstop() > m_ogip_tstop) {
+                m_ogip_tstop = obs->gti().tstop();
+            }
+            m_ogip_telapse  += obs->gti().telapse();
+            m_ogip_exposure += obs->ontime();
+            m_ogip_ontime   += obs->ontime();
+            m_ogip_livetime += obs->livetime();
+        }
+    }
+
     // Return next CTA observation
     return obs;
 }
@@ -63,6 +63,7 @@ public:
     // Methods
     const GObservations& obs(void) const;
+    void                 write_ogip_keywords(GFitsHDU* hdu) const;
 #ifndef SWIG
 protected:
@@ -82,11 +83,18 @@ protected:
     void save_events_xml(void);
     // Protected members
-    GObservations m_obs;            //!< Observation container
+    GObservations m_obs;             //!< Observation container
 private:
     // Private members
-    mutable int m_index_unbinned;   //!< Current index of unbinned observation
+    mutable std::string m_ogip_telescope; //!< Name of telescope
+    mutable GTime       m_ogip_tstart;    //!< Start time for OGIP keywords
+    mutable GTime       m_ogip_tstop;     //!< Stop time for OGIP keywords
+    mutable double      m_ogip_telapse;   //!< Elapsed time
+    mutable double      m_ogip_exposure;  //!< Exposure time
+    mutable double      m_ogip_ontime;    //!< Ontime for OGIP keywords
+    mutable double      m_ogip_livetime;  //!< Livetime for OGIP keywords
+    mutable int         m_index_unbinned; //!< Current index of unbinned observation
 };
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the csmodelsois tool.
+#
+# Copyright (C) 2017 Josh Cardenzana
+#
+# 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/>.
+#
+# ==========================================================================
+import os
+import gammalib
+import cscripts
+from testing import test
+
+
+# =============================== #
+# Test class for csmodelsois tool #
+# =============================== #
+class Test(test):
+    """
+    Test class for csmodelsois tool
+
+    This test class makes unit tests for the csmodelsois tool by using it from
+    the command line and from Python.
+    """
+
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions
+        """
+        # Set test name
+        self.name('csmodelsois')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test csmodelsois on command line')
+        self.append(self._test_python, 'Test csmodelsois from Python')
+
+        # Return
+        return
+
+    # Test csmodelsois on command line
+    def _test_cmd(self):
+        """
+        Test csmodelsois on the command line
+        """
+        # Set tool name
+        csmodelsois = self._script('csmodelsois')
+
+        # Setup csmodelsois command
+        cmd = csmodelsois+' inmodel="'+self._model+'"'+ \
+                          ' outcube="csmodelsois_cmd1.fits"'+\
+                          ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \
+                          ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \
+                          ' soilist="" outmodel=NONE'+ \
+                          ' ra=83.63 dec=22.01 proj=CAR'+ \
+                          ' logfile="csmodelsois_cmd1.log" chatter=1'
+
+        # Check if execution of wrong command fails
+        self.test_assert(self._execute('command_that_does_not_exist') != 0,
+             'Self test of test script')
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check map cube
+        self._check_result_file('csmodelsois_cmd1.fits')
+
+        # Setup csmodelsois command
+        cmd = csmodelsois+' inmodel="events_that_do_not_exist.fits"'+ \
+                          ' outcube="ccsmodelsois_cmd2.fits"'+\
+                          ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \
+                          ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \
+                          ' soilist="" outmodel=NONE'+ \
+                          ' ra=83.63 dec=22.01 proj=CAR'+ \
+                          ' logfile="csmodelsois_cmd2.log"'
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check invalid input file when executed from command line')
+
+        # Check csmodelsois --help
+        self._check_help(csmodelsois)
+
+        # Return
+        return
+
+    # Test csmodelsois from Python
+    def _test_python(self):
+        """
+        Test csmodelsois from Python
+        """
+        # Allocate csmodelsois
+        modelsois = cscripts.csmodelsois()
+
+        # Check that empty csmodelsois tool holds an map cube that has no
+        # energy bins
+        self._check_cube(modelsois.mapcube(), nmaps=0, npixels=0)
+
+        # Check that saving does nothing
+        modelsois['outcube']  = 'csmodelsois_py0.fits'
+        modelsois['outmodel'] = 'NONE'
+        modelsois['logfile']  = 'csmodelsois_py0.log'
+        modelsois.logFileOpen()
+        modelsois.save()
+        self.test_assert(not os.path.isfile('csmodelsois_py0.fits'),
+             'Check that no map cube has been created')
+
+        # Check that clearing does not lead to an exception or segfault
+        modelsois.clear()
+
+        # Now set csmodelsois parameters
+        modelsois['inmodel']  = self._model
+        modelsois['ebinalg']  = 'LOG'
+        modelsois['emin']     = 0.1
+        modelsois['emax']     = 100.0
+        modelsois['enumbins'] = 10
+        modelsois['nxpix']    = 100
+        modelsois['nypix']    = 100
+        modelsois['binsz']    = 0.04
+        modelsois['coordsys'] = 'CEL'
+        modelsois['proj']     = 'CAR'
+        modelsois['ra']       = 83.63
+        modelsois['dec']      = 22.01
+        modelsois['ptsrcsig'] = 1.0
+        modelsois['outcube']  = 'csmodelsois_py1.fits'
+        modelsois['logfile']  = 'csmodelsois_py1.log'
+        modelsois['chatter']  = 2
+        modelsois['soilist']  = ''
+        modelsois['outmodel'] = 'NONE'
+
+        # Run csmodelsois tool
+        modelsois.logFileOpen()   # Make sure we get a log file
+        modelsois.run()
+
+        # Save map cube
+        modelsois.save()
+
+        # Check map cube
+        self._check_result_file('csmodelsois_py1.fits')
+
+        # Allocate csmodelsois scripts, set models and do now a linear
+        # binning; also do not use Gaussian point sources
+        modelsois = cscripts.csmodelsois()
+        modelsois['inmodel']  = self._model
+        modelsois['ebinalg']  = 'LIN'
+        modelsois['emin']     = 0.1
+        modelsois['emax']     = 100.0
+        modelsois['enumbins'] = 10
+        modelsois['nxpix']    = 100
+        modelsois['nypix']    = 100
+        modelsois['binsz']    = 0.04
+        modelsois['coordsys'] = 'CEL'
+        modelsois['proj']     = 'CAR'
+        modelsois['ra']       = 83.63
+        modelsois['dec']      = 22.01
+        modelsois['ptsrcsig'] = 0.0
+        modelsois['outcube']  = 'csmodelsois_py2.fits'
+        modelsois['logfile']  = 'csmodelsois_py2.log'
+        modelsois['chatter']  = 4
+        modelsois['soilist']  = ''
+        modelsois['outmodel'] = 'NONE' 
+
+        # Execute mapcube
+        modelsois.logFileOpen()  # Needed to get a new log file
+        modelsois.execute()
+
+        # Check result file
+        self._check_result_file('csmodelsois_py2.fits')
+
+        # Update output filenames and output a new model file
+        modelsois['outcube']  = 'csmodelsois_py3.fits'
+        modelsois['logfile']  = 'csmodelsois_py3.log'
+        modelsois['outmodel'] = 'csmodelsois_py3.xml'
+
+        # Execute the file
+        modelsois.logFileOpen()
+        modelsois.execute()
+
+        # Now check that the output model file doesnt contain a Crab Sources
+        models3 = gammalib.GModels('csmodelsois_py3.xml')
+        self.test_assert(not models3.contains('Crab'),
+                         'Check Crab model is not present')
+        self.test_assert(models3.contains(modelsois.cubemodelname()), 
+                         'Check cube model is present')
+        self._check_result_file('csmodelsois_py3.fits')
+
+        # Return
+        return
+
+    # Check result file
+    def _check_result_file(self, filename):
+        """
+        Check content of map cube
+
+        Parameters
+        ----------
+        filename : str
+            Map cube file name
+        """
+        # Load map cube
+        cube = gammalib.GModelSpatialDiffuseCube(filename)
+
+        # Check map cube
+        self._check_cube(cube)
+
+        # Return
+        return
+
+    # Check map cube
+    def _check_cube(self, cube, nmaps=11, npixels=10000):
+        """
+        Check map cube
+
+        Parameters
+        ----------
+        cube : `~gammalib.GModelSpatialDiffuseCube`
+            Map cube
+        nmaps : int, optional
+            Number of maps
+        npixels : int, optional
+            Number of pixels
+        """
+        # Get energies (this forces loading in case the map cube is not
+        # loaded)
+        energies = cube.energies()
+        self.test_value(len(cube.energies()), nmaps, 'Check number of energies')
+        
+        # Check dimensions
+        self.test_value(cube.maps(), nmaps, 'Check number of maps')
+        self.test_value(cube.pixels(), npixels, 'Check number of pixels')
+
+        # Return
+        return
@@ -221,7 +221,7 @@ class Test(test):
         skymap['proj']        = 'CAR'
         skymap['xref']        = 184.5575
         skymap['yref']        =  -5.7844
-        skymap['bkgsubtract'] =  'IRF'
+        skymap['bkgsubtract'] = 'IRF'
         skymap['caldb']       = self._caldb
         skymap['irf']         = self._irf
         skymap['outmap']      = 'ctskymap_py4.fits'
@@ -234,6 +234,8 @@ class Test(test):
         # Check result file
         self._check_result_file('ctskymap_py4.fits')
+        self._check_result_file('ctskymap_py4.fits[BACKGROUND]')
+        self._check_result_file('ctskymap_py4.fits[SIGNIFICANCE]')
         # ================================
         # TEST RING BACKGROUND SUBTRACTION
@@ -268,6 +270,8 @@ class Test(test):
         # Check result file
         self._check_result_file('ctskymap_py5.fits', 100, 100)
+        self._check_result_file('ctskymap_py5.fits[BACKGROUND]', 100, 100)
+        self._check_result_file('ctskymap_py5.fits[SIGNIFICANCE]', 100, 100)
         # Return
         return
@@ -32,6 +32,7 @@ import test_cslightcrv
 import test_csmodelinfo
 import test_csmodelmerge
 import test_csmodelselect
+import test_csmodelsois
 import test_csobs2caldb
 import test_csobsdef
 import test_csobsinfo
@@ -121,6 +122,7 @@ def test(installed=False, debug=False):
              test_csmodelinfo.Test(),
              test_csmodelmerge.Test(),
              test_csmodelselect.Test(),
+             test_csmodelsois.Test(),
              test_csobs2caldb.Test(),
              test_csobsdef.Test(),
              test_csobsinfo.Test(),
@@ -145,7 +147,7 @@ def test(installed=False, debug=False):
         # Check for VHEFITS environment variable
         if 'VHEFITS' in os.environ:
-            
+
             # If the environment variable exists then unset it for the test
             # cases. Since Python is executed in a subprocess this will not
             # impact the environment variable in the parent shell.