Integrate code

Jürgen Knödlseder
1 parent ae4014dd
Showing 10 changed files with 181 additions and 154 deletions
ChangeLog
NEWS
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
-2017-08-08
+2017-10-11
         * Version 1.5.0 released
         ========================
+        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
@@ -119,6 +119,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).
 # ==========================================================================
 # 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
@@ -16,11 +16,11 @@
 #
 # General parameters
 #===================
-inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
-outcube,  f, a, "mapcube.fits",,, "Output map cube file"
+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, "NONE",,, "Output model definition XML file"
+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
@@ -52,4 +52,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
 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, "ctmapcube.log",,, "Log filename"
+logfile, f, h, "csmodelsois.log",,, "Log filename"
 #! /usr/bin/env python
 # ==========================================================================
-# Allows generation of DiffuseMapCube model from subset of sources in a
-# model definition XML file
+# Puts subset of sources in diffuse model cube
 #
 # Copyright (C) 2017 Josh Cardenzana
 #
@@ -23,19 +22,24 @@ import sys
 import gammalib
 import ctools
-#=============#
-# csmodelsois #
-#=============#
+
+# ================= #
+# csmodelsois class #
+# ================= #
 class csmodelsois(ctools.cscript):
     """
-    Generates a DiffuseMapCube model from a subset of sources in a model
-    definition XML file. All background models will be ignored!
+    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.
+        Constructor
         Parameters
         ----------
@@ -44,7 +48,7 @@ class csmodelsois(ctools.cscript):
         Raises
         ------
         TypeError
-        An invalid number of command line arguments was provided.
+            An invalid number of command line arguments was provided.
         """
         # Set name and version
         self._name    = 'csmodelsois'
@@ -60,56 +64,60 @@ class csmodelsois(ctools.cscript):
         return
+    # Private methods
     def _init_members(self):
         """
         Initialize the parameters for this object
         """
-        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
+        # 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 _set_cube_parameters(self):
+    def _get_parameters(self):
         """
-        Fill in the parameters associated with this object
+        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() )
+        self.cubegen['coordsys'].string(self['coordsys'].string())
         if self['coordsys'].string() == 'CEL':
-            # Center in celestial coordinates
-            self.cubegen['xref'].real( self['ra'].real() )
-            self.cubegen['yref'].real( self['dec'].real() )
+            self.cubegen['xref'].real(self['ra'].real())
+            self.cubegen['yref'].real(self['dec'].real())
         else:
-            # Otherwise assume center in galactic coordinates
-            self.cubegen['xref'].real( self['glon'].real() )
-            self.cubegen['yref'].real( self['glat'].real() )
+            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() )
+        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":
-            # Read energy binning from a file
-            self.cubegen['ebinfile'].filename( self['ebinfile'].filename() )
+        self.cubegen['ebinalg'].string(self['ebinalg'].string())
+        if self['ebinalg'].string() == 'FILE':
+            self.cubegen['ebinfile'].filename(self['ebinfile'].filename())
         else:
-            # Create the energy binning information
-            self.cubegen['emin'].real( self['emin'].real() )
-            self.cubegen['emax'].real( self['emax'].real() )
-            self.cubegen['enumbins'].integer( self['enumbins'].integer() )
+            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( self['ptsrcsig'].real() )
+        self.cubegen['ptsrcsig'].real(ptsrcsig)
         # Read optionally output cube filenames
         if self._read_ahead():
-            self.cubegen['outcube'].filename( self['outcube'].filename() )
-            self['soilist'].string()
+            self.cubegen['outcube'].filename(self['outcube'].filename())
             self['outmodel'].filename()
         # Write parameters into logger
@@ -117,11 +125,9 @@ class csmodelsois(ctools.cscript):
         return
-
     def _gen_cubemodel(self):
         """
         Generates a binned model from the model file using ctmapcube
-        returns: updated model filename
         """
         # Get a list of sources we DONT want to put into the cube
         sources = self['soilist'].string().split(',')
@@ -129,52 +135,32 @@ class csmodelsois(ctools.cscript):
         # 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:
-            # Convert the model to a GModelSky
-            mod = model.classname()
-            # Check the model
-            if mod != "GModelSky":
-                # Model is most likely a background model, so continue
+            # 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:
-                # Model is a source of interest, so remove it from list of models
-                # to be used in generating the cube
                 self.cubemodels.remove(model.name())
         # Log the number of models to be put into the cube
-        self._log_value(self['chatter'].integer(), 'Cube models', self.cubemodels.size())
-
-        return
-
-
-    def mapcube(self):
-        """
-        Return the mapcube generated by the underlying 'ctmapcube' object
-        """
-        return self.cubegen.mapcube()
+        self._log_value(self['chatter'].integer(), 'Numner of cube models',
+                        self.cubemodels.size())
-
-    def cubemodelname(self):
-        """
-        Return the name to be given to the cube model in the output XML file
-        """
-        return "FixedSourcesCube"
-
-
-    def models(self, models:gammalib.GModels):
-        """
-        Set the models for analysis
-        """
-        self.m_models = models
+        # Return
         return
+    # Public methods
     def run(self):
         """
         Implements the actual bulk of the script's tasks
@@ -184,41 +170,62 @@ class csmodelsois(ctools.cscript):
             self._log.cout(True)
         # Get parameters
-        self._set_cube_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():
-            self._log_string(gammalib.NORMAL, 'Generting model cube')
+
+            # Log action
+            self._log_string(gammalib.NORMAL, 'Generating model cube')
             # Fill the cube models into the ctmapcube object
-            self.cubegen.models( self.cubemodels )
+            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')
+            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
         """
-        #cube = gammalib.GModelSpatialDiffuseCube(self.cubegen.mapcube())
-
         # 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() )
+            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() ):
+        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)
@@ -227,35 +234,47 @@ class csmodelsois(ctools.cscript):
                 outmodels.remove(model.name())
             # Generate the actual cube model
-            spat     = gammalib.GModelSpatialDiffuseCube( self['outcube'].filename() )
+            spat     = gammalib.GModelSpatialDiffuseCube(self['outcube'].filename())
             spec     = gammalib.GModelSpectralConst()
             newmodel = gammalib.GModelSky(spat, spec)
             # Give the model a default name
-            newmodel.name( self.cubemodelname() )
+            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() )
+            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 execute(self):
+    def cubemodelname(self):
         """
-        Execute the script
+        Return the name to be given to the cube model in the output XML file
         """
-        # Open logfile
-        self.logFileOpen()
+        # Return
+        return 'FixedSourcesCube'
-        # Read ahead output parameters
-        self._read_ahead(True)
+    def models(self, models):
+        """
+        Set the models
-        # Run the script
-        self.run()
-        self.save()
+        Parameters
+        -------
+        models : `~gammalib.GModels`
+            Model container
+        """
+        # Set models
+        self.m_models = models
         # Return
         return
@@ -28,5 +28,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
@@ -35,7 +35,7 @@ General parameters
     Input model XML file.
 ``outmodel [file]``
-    Output model XML file (not generated if set to 'NONE').
+    Output model XML file (not generated if set to ``NONE``).
 ``outcube [file]``
     Output map cube file.
@@ -48,37 +48,37 @@ 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.
-                 
+
 ``ra [real]``
     Right Ascension of cube centre (J2000, in degrees).
-                 
+
 ``dec [real]``
     Declination of cube centre (J2000, in degrees).
@@ -88,6 +88,7 @@ General parameters
 ``glat [real]``
     Galactic latitude of cube centre (in degrees).
+
 Standard parameters
 -------------------
@@ -97,21 +98,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").
@@ -122,4 +123,4 @@ Standard parameters
 Related tools or scripts
 ------------------------
-ctmapcube
 \ No newline at end of file
+: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
 ------------------------
-csmodelsois
+:ref:`csmodelsois`
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csmodelsois tool.
 #
-# Copyright (C) 2016-2017 Juergen Knoedlseder
+# 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
@@ -71,12 +71,12 @@ class Test(test):
         # 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'
+                          ' 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,
@@ -91,12 +91,12 @@ class Test(test):
         # 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"'
+                          ' 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,
@@ -202,7 +202,8 @@ class Test(test):
         # 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(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')
@@ -147,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.