diff --git a/ChangeLog b/ChangeLog
index 65bf1cd..2960fe1 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,11 @@
-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)
diff --git a/NEWS b/NEWS
index 2940dcf..49b1fc0 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
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).
diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index 2adc17d..d0bd9d1 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -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 \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index e079b96..09f8982 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -1,7 +1,7 @@
# ==========================================================================
# 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
diff --git a/cscripts/csmodelsois.par b/cscripts/csmodelsois.par
new file mode 100644
index 0000000..62a14bc
--- /dev/null
+++ b/cscripts/csmodelsois.par
@@ -0,0 +1,55 @@
+#############################################################################
+# #
+# 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"
diff --git a/cscripts/csmodelsois.py b/cscripts/csmodelsois.py
new file mode 100755
index 0000000..37eacbb
--- /dev/null
+++ b/cscripts/csmodelsois.py
@@ -0,0 +1,292 @@
+#! /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 .
+#
+# ==========================================================================
+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()
diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst
index 9e4ab6b..1debf07 100644
--- a/doc/source/admin/release_history/1.5.rst
+++ b/doc/source/admin/release_history/1.5.rst
@@ -22,11 +22,16 @@ Bug fixes
Improvements
------------
-* TBW
+* [`2194 `_] -
+ Add keywords to sky maps produced by :ref:`ctskymap`
+* [`2189 `_] -
+ Enhance efficiency of computations in :ref:`ctmodel`
New features
------------
+* [`2180 `_] -
+ Add :ref:`csmodelsois` script
* [`1864 `_] -
Use t-type parameters for all ``tmin`` and ``tmax`` parameters
diff --git a/doc/source/users/reference_manual/csmodelsois.rst b/doc/source/users/reference_manual/csmodelsois.rst
new file mode 100644
index 0000000..3dfff45
--- /dev/null
+++ b/doc/source/users/reference_manual/csmodelsois.rst
@@ -0,0 +1,126 @@
+.. _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 [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 [string]``
+ Coordinate system (CEL - celestial, GAL - galactic).
+
+``proj [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`
diff --git a/doc/source/users/reference_manual/ctmapcube.rst b/doc/source/users/reference_manual/ctmapcube.rst
index e94f81e..741b17b 100644
--- a/doc/source/users/reference_manual/ctmapcube.rst
+++ b/doc/source/users/reference_manual/ctmapcube.rst
@@ -36,40 +36,40 @@ General parameters
``ebinalg [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 [string]``
Coordinate system (CEL - celestial, GAL - galactic).
-
+
``proj [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`
diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py
index 18f18be..5a5a6f9 100755
--- a/examples/show_spectrum.py
+++ b/examples/show_spectrum.py
@@ -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}$)')
diff --git a/pyext/ctobservation.i b/pyext/ctobservation.i
index 5791958..67e9966 100644
--- a/pyext/ctobservation.i
+++ b/pyext/ctobservation.i
@@ -1,7 +1,7 @@
/***************************************************************************
* 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;
diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp
index 332e0a0..5be5ff0 100644
--- a/src/ctmodel/ctmodel.cpp
+++ b/src/ctmodel/ctmodel.cpp
@@ -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(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(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(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(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;
+}
diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp
index 40b925b..2ce9564 100644
--- a/src/ctmodel/ctmodel.hpp
+++ b/src/ctmodel/ctmodel.hpp
@@ -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
diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp
index fd0a134..bc4fdce 100644
--- a/src/ctskymap/ctskymap.cpp
+++ b/src/ctskymap/ctskymap.cpp
@@ -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
diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp
index 65ac034..72facae 100644
--- a/src/ctskymap/ctskymap.hpp
+++ b/src/ctskymap/ctskymap.hpp
@@ -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 {
diff --git a/src/support/ctobservation.cpp b/src/support/ctobservation.cpp
index 4006a9b..2e03c6d 100644
--- a/src/support/ctobservation.cpp
+++ b/src/support/ctobservation.cpp
@@ -192,6 +192,53 @@ ctobservation& ctobservation::operator=(const ctobservation& 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& 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;
}
diff --git a/src/support/ctobservation.hpp b/src/support/ctobservation.hpp
index 3e6dc7a..7b0bee0 100644
--- a/src/support/ctobservation.hpp
+++ b/src/support/ctobservation.hpp
@@ -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
};
diff --git a/test/test_csmodelsois.py b/test/test_csmodelsois.py
new file mode 100755
index 0000000..79ed9dd
--- /dev/null
+++ b/test/test_csmodelsois.py
@@ -0,0 +1,257 @@
+#! /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 .
+#
+# ==========================================================================
+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
diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py
index 80e5ca9..7721701 100755
--- a/test/test_ctskymap.py
+++ b/test/test_ctskymap.py
@@ -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
diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py
index b674df8..0ffbfe0 100755
--- a/test/test_python_cscripts.py
+++ b/test/test_python_cscripts.py
@@ -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.