Merge branch '2180-create_csmodelsois' into integration

# Conflicts: # ChangeLog

Merge branch '2180-create_csmodelsois' into integration
# Conflicts: # ChangeLog
Jürgen Knödlseder
2 parents c5d45282 80f82c6d
Showing 11 changed files with 766 additions and 22 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
test/test_csmodelsois.py
test/test_python_cscripts.py
@@ -5,6 +5,7 @@
  
         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)
  
  
@@ -126,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()
@@ -31,5 +31,7 @@ Improvements
 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`
+#! /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
@@ -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.