Enhance ctools unit tests

Jürgen Knödlseder
1 parent ce48d505
Showing 21 changed files with 335 additions and 99 deletions
src/ctbkgcube/ctbkgcube.cpp
src/ctbutterfly/ctbutterfly.cpp
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
src/ctedispcube/ctedispcube.cpp
src/cterror/cterror.cpp
src/ctexpcube/ctexpcube.cpp
src/ctlike/ctlike.cpp
src/ctmodel/ctmodel.cpp
src/ctobssim/ctobssim.cpp
src/ctpsfcube/ctpsfcube.cpp
src/cttsmap/cttsmap.cpp
src/ctulimit/ctulimit.cpp
src/support/ctool.cpp
test/test_ctbin.py
test/test_ctbkgcube.py
test/test_ctcubemask.py
test/test_ctedispcube.py
test/test_ctexpcube.py
test/test_ctpsfcube.py
@@ -488,6 +488,12 @@ void ctbkgcube::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -573,9 +573,8 @@ void ctbutterfly::free_members(void)
 /***********************************************************************//**
  * @brief Get application parameters
  *
- * @exception GException::invalid_value
- *            Test source not found or no RA/DEC parameters found for test
- *            source.
+ * @exception GException::feature_not_implemented
+ *            Loading of covariance matrix is not yet implemented.
  *
  * Get all task parameters from parameter file or (if required) by querying
  * the user. Most parameters are only required if no observation exists so
@@ -597,6 +596,11 @@ void ctbutterfly::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
     // If there is are no models associated with the observations then
     // load now the model definition
@@ -301,7 +301,7 @@ void ctcubemask::save(void)
     // Get counts cube filename
     m_outcube = (*this)["outcube"].filename();
-    // Case A: Save event file(s) and XML metadata information
+    // Case A: Save counts cube(s) and XML metadata information
     if (m_use_xml) {
         // Get prefix
@@ -311,7 +311,7 @@ void ctcubemask::save(void)
         save_xml();
     }
-    // Case B: Save event file as FITS file
+    // Case B: Save counts cube as FITS file
     else {
         save_fits();
     }
@@ -404,9 +404,9 @@ void ctcubemask::init_members(void)
 void ctcubemask::copy_members(const ctcubemask& app)
 {
     // Copy parameters
-	m_regfile = app.m_regfile;
+    m_regfile = app.m_regfile;
     m_outcube = app.m_outcube;
-	m_prefix  = app.m_prefix;
+    m_prefix  = app.m_prefix;
     m_usepnt  = app.m_usepnt;
     m_ra      = app.m_ra;
     m_dec     = app.m_dec;
@@ -471,7 +471,7 @@ void ctcubemask::get_parameters(void)
     } // endif: there was no observation in the container
     // Get parameters
-	m_regfile = (*this)["regfile"].filename();
+    m_regfile = (*this)["regfile"].filename();
     m_usepnt  = (*this)["usepnt"].boolean();
     if (!m_usepnt) {
@@ -725,18 +725,23 @@ std::string ctcubemask::set_outfile_name(const std::string&amp; filename) const
  ***************************************************************************/
 void ctcubemask::save_fits(void)
 {
-    // Save only if filename is non-empty
-    if (!m_outcube.is_empty()) {
+    // Save only if filename is non-empty and if there are observations
+    if (!m_outcube.is_empty() && m_obs.size() > 0) {
         // Get CTA observation from observation container
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
-        // Log counts cube file name
-        log_value(NORMAL, "Counts cube file", m_outcube.url());
+        // Handle only CTA observations
+        if (obs != NULL) {
+
+            // Log counts cube file name
+            log_value(NORMAL, "Counts cube file", m_outcube.url());
-        // Save event list
-        save_counts_map(obs, m_outcube);
+            // Save counts cube
+            obs->save(m_outcube, clobber());
+        } // endif: observation was a CTA observation
+
     }
     // Return
@@ -745,30 +750,20 @@ void ctcubemask::save_fits(void)
 /***********************************************************************//**
- * @brief Save counts map(s) in XML format.
+ * @brief Save counts cube(s) in XML format.
  *
- * Save the counts map(s) into FITS files and write the file path information
- * into a XML file. The filename of the XML file is specified by the
- * m_outfile member, the filename(s) of the counts map(s) are built by
+ * Save the counts cube(s) into FITS files and write the file path
+ * information into an XML file. The filename of the XML file is specified by
+ * the m_outfile member, the filename(s) of the counts cube(s) are built by
  * prepending the prefix given by the m_prefix member to the input counts
- * map(s) filenames. Any path present in the input filename will be stripped,
- * i.e. the counts map(s) will be written in the local working directory
- * (unless a path is specified in the m_prefix member).
+ * cube(s) filenames. Any path present in the input filename will be
+ * stripped, i.e. the counts cube(s) will be written in the local working
+ * directory (unless a path is specified in the m_prefix member).
  ***************************************************************************/
 void ctcubemask::save_xml(void)
 {
     // Issue warning if output filename has no .xml suffix
-    std::string suffix = gammalib::tolower(m_outcube.url().substr(
-                                                    m_outcube.length()-4,4));
-    if (suffix != ".xml") {
-        std::string msg = "\nWARNING: Name of observation definition output "
-                          "file \""+m_outcube+"\" does not terminate with\n"
-                          "         \".xml\". This is not an error, but may "
-                          "be misleading. It is recommended to use the\n"
-                          "         suffix \".xml\" for observation definition "
-                          "files.";
-        log_string(TERSE, msg);
-    }
+    log_string(TERSE, warn_xml_suffix(m_outcube.url()));
     // Loop over all observation in the container
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -788,8 +783,8 @@ void ctcubemask::save_xml(void)
             // Log counts cube file name
             log_value(NORMAL, "Counts cube file", outfile);
-            // Save event list
-            save_counts_map(obs, outfile);
+            // Save counts cube
+            obs->save(outfile, clobber());
         } // endif: observation was a CTA observations
@@ -801,30 +796,3 @@ void ctcubemask::save_xml(void)
     // Return
     return;
 }
-
-
-
-/***********************************************************************//**
- * @brief Save a single counts map into a FITS file
- *
- * @param[in] obs Pointer to CTA observation.
- * @param[in] outfile Output file name.
- *
- * This method saves a single counts map into a FITS file. The method does
- * nothing if the observation pointer is not valid.
- ***************************************************************************/
-void ctcubemask::save_counts_map(const GCTAObservation* obs,
-                                 const std::string&     outfile) const
-{
-    // Save only if observation is valid
-    if (obs != NULL) {
-
-        // Save observation into FITS file
-        obs->save(outfile, clobber());
-
-    } // endif: observation was valid
-
-    // Return
-    return;
-}
-
@@ -36,7 +36,7 @@
 /* __Definitions _________________________________________________________ */
 #define CTCUBEMASK_NAME    "ctcubemask"
-#define CTCUBEMASK_VERSION "1.1.0"
+#define CTCUBEMASK_VERSION "1.2.0"
 /***********************************************************************//**
@@ -75,8 +75,6 @@ protected:
     std::string set_outfile_name(const std::string& filename) const;
     void        save_fits(void);
     void        save_xml(void);
-    void        save_counts_map(const GCTAObservation* obs,
-                                const std::string&     outfile) const;
     // User parameters
 	GFilename   m_regfile;    //!< ds9 region file
@@ -361,6 +361,12 @@ void ctedispcube::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -475,6 +475,11 @@ void cterror::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
     // If there are no models associated with the observations then
     // load now the model definition
@@ -386,6 +386,12 @@ void ctexpcube::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -373,6 +373,12 @@ void ctlike::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // If only single observation is used, read statistics parameter
     if (!m_use_xml) {
@@ -456,6 +456,12 @@ void ctmodel::get_parameters(void)
         get_obs();
     }
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // If we have now excactly one CTA observation (but no cube has yet been
     // appended to the observation) then check whether this observation
     // is a binned observation, and if yes, extract the counts cube for
@@ -561,9 +561,10 @@ void ctobssim::get_parameters(void)
     }
-    // ... otherwise make sure that observation boundaries are set
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
     else {
-        set_obs_bounds(m_obs);
+        setup_observations(m_obs);
     }
     // Read model definition file if required
@@ -346,6 +346,12 @@ void ctpsfcube::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // Get the incube filename
     std::string incube = (*this)["incube"].filename();
@@ -523,6 +523,12 @@ void cttsmap::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
     // If there is are no models associated with the observations then
     // load now the model definition
     if (m_obs.models().size() == 0) {
@@ -481,6 +481,11 @@ void ctulimit::get_parameters(void)
     } // endif: there was no observation in the container
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
     // If there is are no models associated with the observations then
     // load now the model definition
@@ -1646,7 +1646,7 @@ std::string ctool::warn_xml_suffix(const GFilename&amp; filename) const
     // Warn if there are not enough energy bins
     if (suffix != ".xml") {
-        warning.append("\nWARNING: Name of observation definition output file "
+        warning.append("WARNING: Name of observation definition output file "
                        "\""+filename+"\" does not terminate with \".xml\"."
                        "\n         This is not an error, but may be misleading. "
                        "It is recommended to use the suffix"
@@ -194,15 +194,9 @@ class Test(test):
         self.test_value(bin.cube().size(), 0,
              'Check that empty ctbin has an empty counts cube')
-        # Prepare observation container for stacked analysis
-        cta = gammalib.GCTAObservation(self._events)
-        obs = gammalib.GObservations()
-        cta.id('0001')
-        obs.append(cta)
-        cta.id('0002')
-        obs.append(cta)
-        cta.id('0003')
-        obs.append(cta)
+        # Prepare observation container for stacking of events into a
+        # single counts cube
+        obs = self._obs_events()
         # Set-up ctbin using an observation container
         bin = ctools.ctbin(obs)
@@ -213,18 +213,8 @@ class Test(test):
         self._check_cube(cpy_bkgcube.bkgcube(), nenergies=0)
         self._check_models(cpy_bkgcube.models(), nmodels=0)
-        # Prepare observation container
-        obs1 = gammalib.GCTAObservation(self._events)
-        obs2 = gammalib.GCTAObservation(self._cntcube)
-        obs3 = gammalib.GLATObservation()
-        obs1.id('0001')
-        obs2.id('0002')
-        obs3.id('0001')
-        obs3.events(gammalib.GLATEventList())
-        obs = gammalib.GObservations()
-        obs.append(obs1)
-        obs.append(obs2)
-        obs.append(obs3)
+        # Get mixel observation container
+        obs = self._obs_mixed()
         obs.models(gammalib.GModels(self._model))
         # Set-up ctbkgcube from observation container
@@ -18,6 +18,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
 import gammalib
 import ctools
 from testing import test
@@ -106,8 +107,24 @@ class Test(test):
         """
         Test ctcubemask from Python
         """
-        # Set-up ctcubemask
+        # Allocate ctcubemask
         mask = ctools.ctcubemask()
+
+        # Check that empty ctcubemask tool holds an empty observation
+        self._check_obs(mask.obs(), nobs=0)
+
+        # Check that saving saves an empty counts cube
+        mask['outcube'] = 'ctcubemask_py0.fits'
+        mask['logfile'] = 'ctcubemask_py0.log'
+        mask.logFileOpen()
+        mask.save()
+        self.test_assert(not os.path.isfile('ctcubemask_py0.fits'),
+             'Check that no counts cube has been created')
+
+        # Check that clearing does not lead to an exception or segfault
+        mask.clear()
+
+        # Now set ctcubemask parameters
         mask['inobs']   = self._cntcube
         mask['regfile'] = self._exclusion
         mask['ra']      = 83.63
@@ -127,14 +144,15 @@ class Test(test):
         # Check result file
         self._check_result_file('ctcubemask_py1.fits')
-        # Set-up ctcubemask without exclusion regions
+        # Set-up ctcubemask without exclusion regions but tighter energy
+        # selection
         mask = ctools.ctcubemask()
         mask['inobs']   = self._cntcube
         mask['regfile'] = 'NONE'
         mask['ra']      = 83.63
         mask['dec']     = 22.01
         mask['rad']     = 3.0
-        mask['emin']    = 0.1
+        mask['emin']    = 1.0
         mask['emax']    = 100.0
         mask['outcube'] = 'ctcubemask_py2.fits'
         mask['logfile'] = 'ctcubemask_py2.log'
@@ -146,22 +164,54 @@ class Test(test):
         mask.execute()
         # Check result file
-        self._check_result_file('ctcubemask_py2.fits', events=5542)
+        self._check_result_file('ctcubemask_py2.fits', events=573)
+
+        # Copy ctcubemask tool
+        cpy_mask = mask.copy()
-        # Copy ctcubemask tool and execute copy
-        cpy_mask = mask
+        # Check that ctcubemask tool holds one observation
+        self._check_obs(cpy_mask.obs())
+
+        # Execute copy of ctcubemask tool again, now with a higher chatter
+        # level than before
         cpy_mask['outcube'] = 'ctcubemask_py3.fits'
         cpy_mask['logfile'] = 'ctcubemask_py3.log'
-        cpy_mask['chatter']  = 4
+        cpy_mask['chatter'] = 4
         cpy_mask.execute()
         # Check result file
-        self._check_result_file('ctcubemask_py3.fits', events=5542)
+        self._check_result_file('ctcubemask_py3.fits', events=573)
         # Clear ctcubemask tool
-        mask.clear()
+        cpy_mask.clear()
-        # TODO: Do some test after clearing
+        # Check that empty ctcubemask tool holds an empty observation
+        self._check_obs(cpy_mask.obs(), nobs=0)
+
+        # Prepare observation container
+        obs = self._obs_mixed()
+        obs.models(gammalib.GModels(self._model))
+
+        # Set-up ctcubemask from observation container, don't perform any
+        # energy selection
+        mask = ctools.ctcubemask(obs)
+        mask['regfile'] = 'NONE'
+        mask['ra']      = 83.63
+        mask['dec']     = 22.01
+        mask['rad']     = 3.0
+        mask['emin']    = 'NONE'
+        mask['emax']    = 'NONE'
+        mask['outcube'] = 'ctcubemask_py4.xml'
+        mask['logfile'] = 'ctcubemask_py4.log'
+        mask['chatter'] = 3
+        mask['debug']   = True
+
+        # Execute ctcubemask tool
+        mask.logFileOpen()   # Make sure we get a log file
+        mask.execute()
+
+        # Check result file
+        self._check_result_file('filtered_crab_cntmap.fits', events=5542)
         # Return
         return
@@ -180,3 +230,21 @@ class Test(test):
         # Return
         return
+
+    # Check observation container
+    def _check_obs(self, obs, nobs=1):
+        """
+        Check observation container
+
+        Parameters
+        ----------
+        obs : `~gammalib.GObservations`
+            Observation container
+        nobs : int, optional
+            Expected number of observations
+        """
+        # Check size of observation container
+        self.test_value(obs.size(), nobs, 'Check size of observation container')
+
+        # Return
+        return
@@ -185,6 +185,38 @@ class Test(test):
         # cube
         self._check_cube(cpy_edispcube.edispcube(), nenergies=0, nmigras=0)
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Set-up ctedispcube from observation container
+        edispcube = ctools.ctedispcube(obs)
+        edispcube['incube']    = 'NONE'
+        edispcube['caldb']     = self._caldb
+        edispcube['irf']       = self._irf
+        edispcube['ebinalg']   = 'LOG'
+        edispcube['emin']      = 0.1
+        edispcube['emax']      = 100.0
+        edispcube['enumbins']  = 20
+        edispcube['nxpix']     = 10
+        edispcube['nypix']     = 10
+        edispcube['binsz']     = 0.4
+        edispcube['coordsys']  = 'CEL'
+        edispcube['proj']      = 'CAR'
+        edispcube['xref']      = 83.63
+        edispcube['yref']      = 22.01
+        edispcube['migramax']  = 2.0
+        edispcube['migrabins'] = 10
+        edispcube['outcube']   = 'ctedispcube_py3.fits'
+        edispcube['logfile']   = 'ctedispcube_py3.log'
+        edispcube['chatter']   = 4
+
+        # Execute ctedispcube tool
+        edispcube.logFileOpen()  # Make sure we get a log file
+        edispcube.execute()
+
+        # Check result file
+        self._check_result_file('ctedispcube_py3.fits')
+
         # Return
         return
@@ -181,6 +181,36 @@ class Test(test):
         # Check that the cleared copy has also cleared the exposure cube
         self._check_cube(cpy_expcube.expcube(), nenergies=0)
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Set-up ctexpcube from observation container
+        expcube = ctools.ctexpcube(obs)
+        expcube['incube']   = 'NONE'
+        expcube['caldb']    = self._caldb
+        expcube['irf']      = self._irf
+        expcube['ebinalg']  = 'LOG'
+        expcube['emin']     = 0.1
+        expcube['emax']     = 100.0
+        expcube['enumbins'] = 20
+        expcube['nxpix']    = 200
+        expcube['nypix']    = 200
+        expcube['binsz']    = 0.02
+        expcube['coordsys'] = 'CEL'
+        expcube['proj']     = 'CAR'
+        expcube['xref']     = 83.63
+        expcube['yref']     = 22.01
+        expcube['logfile']  = 'ctexpcube_py3.log'
+        expcube['outcube']  = 'ctexpcube_py3.fits'
+        expcube['chatter']  = 4
+
+        # Execute ctexpcube tool
+        expcube.logFileOpen()  # Needed to get a new log file
+        expcube.execute()
+
+        # Check result file
+        self._check_result_file('ctexpcube_py3.fits')
+
         # Return
         return
@@ -131,7 +131,6 @@ class Test(test):
         # Now set ctpsfcube parameters
         psfcube['inobs']    = self._events
         psfcube['incube']   = 'NONE'
-        psfcube['outcube']  = 'ctpsfcube_py1.fits'
         psfcube['caldb']    = self._caldb
         psfcube['irf']      = self._irf
         psfcube['ebinalg']  = 'LOG'
@@ -147,6 +146,7 @@ class Test(test):
         psfcube['yref']     = 22.01
         psfcube['amax']     = 0.3
         psfcube['anumbins'] = 10
+        psfcube['outcube']  = 'ctpsfcube_py1.fits'
         psfcube['logfile']  = 'ctpsfcube_py1.log'
         psfcube['chatter']  = 2
@@ -184,6 +184,38 @@ class Test(test):
         # Check that the cleared copy has also cleared the PSF cube
         self._check_cube(cpy_psfcube.psfcube(), nenergies=0)
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Set-up ctpsfcube from observation container
+        psfcube = ctools.ctpsfcube(obs)
+        psfcube['incube']   = 'NONE'
+        psfcube['caldb']    = self._caldb
+        psfcube['irf']      = self._irf
+        psfcube['ebinalg']  = 'LOG'
+        psfcube['emin']     = 0.1
+        psfcube['emax']     = 100
+        psfcube['enumbins'] = 20
+        psfcube['nxpix']    = 10
+        psfcube['nypix']    = 10
+        psfcube['binsz']    = 0.4
+        psfcube['coordsys'] = 'CEL'
+        psfcube['proj']     = 'CAR'
+        psfcube['xref']     = 83.63
+        psfcube['yref']     = 22.01
+        psfcube['amax']     = 0.3
+        psfcube['anumbins'] = 10
+        psfcube['outcube']  = 'ctpsfcube_py3.fits'
+        psfcube['logfile']  = 'ctpsfcube_py3.log'
+        psfcube['chatter']  = 4
+
+        # Execute ctpsfcube tool
+        psfcube.logFileOpen()   # Make sure we get a log file
+        psfcube.execute()
+
+        # Check result file
+        self._check_result_file('ctpsfcube_py3.fits')
+
         # Return
         return
@@ -97,3 +97,70 @@ class test(gammalib.GPythonTestSuite):
         # Return return code
         return rc
+
+    # Return observation container with 3 identical event lists
+    def _obs_events(self):
+        """
+        Return observation container with 3 identical event lists
+
+        Returns
+        -------
+        obs : `~gammalib.gammalib.GObservations`
+            Observations container
+        """
+        # Load CTA event list
+        cta = gammalib.GCTAObservation(self._events)
+
+        # Allocate observation container
+        obs = gammalib.GObservations()
+
+        # Append CTA observation with ID "0001"
+        cta.id('0001')
+        obs.append(cta)
+
+        # Append CTA observation with ID "0002"
+        cta.id('0002')
+        obs.append(cta)
+
+        # Append CTA observation with ID "0003"
+        cta.id('0003')
+        obs.append(cta)
+
+        # Return observation container
+        return obs
+
+    # Return observation container with event list, counts cube and LAT
+    # observation
+    def _obs_mixed(self):
+        """
+        Return observation container with event list, counts cube and LAT
+        observation
+
+        Returns
+        -------
+        obs : `~gammalib.gammalib.GObservations`
+            Observations container
+        """
+        # Load CTA event list
+        obs1 = gammalib.GCTAObservation(self._events)
+        obs1.id('0001')
+
+        # Load CTA counts cube
+        obs2 = gammalib.GCTAObservation(self._cntcube)
+        obs2.id('0002')
+
+        # Create empty LAT observation
+        obs3 = gammalib.GLATObservation()
+        obs3.id('0001')
+        obs3.events(gammalib.GLATEventList())
+
+        # Allocate observation container
+        obs = gammalib.GObservations()
+
+        # Append observations to container
+        obs.append(obs1)
+        obs.append(obs2)
+        obs.append(obs3)
+
+        # Return observation container
+        return obs