Increase ctool branch coverage (#1828)

Jürgen Knödlseder
1 parent 1e00ea3d
Showing 12 changed files with 404 additions and 122 deletions
src/cterror/cterror.cpp
src/cterror/cterror.par
src/ctulimit/ctulimit.cpp
test/test_ctedispcube.py
test/test_cterror.py
test/test_ctlike.py
test/test_ctmapcube.py
test/test_ctmodel.py
test/test_ctobssim.py
test/test_ctpsfcube.py
test/test_ctskymap.py
test/test_ctulimit.py
@@ -250,85 +250,80 @@ void cterror::run(void)
     log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_best_logL,3));
     log_string(NORMAL, m_obs.models().print(m_chatter));
-    // Continue only if source model exists
-    if (m_obs.models().contains(m_srcname)) {
+    // Save best fitting models
+    GModels models_best = m_obs.models();
-        // Save best fitting models
-        GModels models_best = m_obs.models();
+    // Get pointer on model
+    GModel* model = models_best[m_srcname];
-        // Get pointer on model
-        GModel* model = models_best[m_srcname];
+    // Get number of parameters
+    int npars = model->size();
-        // Get number of parameters
-        int npars = model->size();
+    // Loop over parameters of sky model
+    for (int i = 0; i < npars; ++i) {
-        // Loop over parameters of sky model
-        for (int i = 0; i < npars; ++i) {
-
-            // Skip parameter if it is fixed
-            if (model->at(i).is_fixed()) {
-                continue;
-            }
+        // Skip parameter if it is fixed
+        if (model->at(i).is_fixed()) {
+            continue;
+        }
-            // Initialise with best fitting models
-            m_obs.models(models_best);
-
-            // Get pointer on model parameter
-            GModels& current_models = const_cast<GModels&>(m_obs.models());
-            m_model_par             = &(current_models[m_srcname]->at(i));
-
-            // Extract current value
-            m_value = m_model_par->factor_value();
-
-            // Compute parameter bracketing
-            double parmin = std::max(m_model_par->factor_min(),
-                                     m_value - 10.0*m_model_par->factor_error());
-            double parmax = std::min(m_model_par->factor_max(),
-                                     m_value + 10.0*m_model_par->factor_error());
-
-            // Write header and initial parameters into logger
-            log_header1(TERSE, "Compute error for source \""+m_srcname+"\""
-                               " parameter \""+m_model_par->name()+"\"");
-            log_value(NORMAL, "Confidence level",
-                      gammalib::str(m_confidence*100.0)+" %");
-            log_value(NORMAL, "Log-likelihood difference", m_dlogL);
-            log_value(NORMAL, "Initial factor range",
-                      "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
-
-            // Compute lower and upper boundaries
-            double value_lo = error_bisection(parmin, m_value);
-            double value_hi = error_bisection(m_value, parmax);
-
-            // Compute errors
-            double error           = 0.5 * (value_hi - value_lo);
-            double error_neg       = m_value  - value_lo;
-            double error_pos       = value_hi - m_value;
-            double error_value     = std::abs(error*m_model_par->scale());
-            double error_value_neg = std::abs(error_neg*m_model_par->scale());
-            double error_value_pos = std::abs(error_pos*m_model_par->scale());
-
-            // Write results into logger
-            std::string unit = " " + m_model_par->unit();
-            log_value(NORMAL, "Lower parameter factor", value_lo);
-            log_value(NORMAL, "Upper parameter factor", value_hi);
-            log_value(NORMAL, "Error from curvature",
-                      gammalib::str(m_model_par->error()) + unit);
-            log_value(NORMAL, "Error from profile",
-                      gammalib::str(error_value) + unit);
-            log_value(NORMAL, "Negative profile error",
-                      gammalib::str(error_value_neg) + unit);
-            log_value(NORMAL, "Positive profile error",
-                      gammalib::str(error_value_pos) + unit);
-
-            // Save error result
-            model->at(i).factor_error(error);
-
-        } // endfor: looped over spectral parameters
-
-        // Restore best fitting models (now with new errors computed)
+        // Initialise with best fitting models
         m_obs.models(models_best);
-    } // endif: source model exists
+        // Get pointer on model parameter
+        GModels& current_models = const_cast<GModels&>(m_obs.models());
+        m_model_par             = &(current_models[m_srcname]->at(i));
+
+        // Extract current value
+        m_value = m_model_par->factor_value();
+
+        // Compute parameter bracketing
+        double parmin = std::max(m_model_par->factor_min(),
+                                 m_value - 10.0*m_model_par->factor_error());
+        double parmax = std::min(m_model_par->factor_max(),
+                                 m_value + 10.0*m_model_par->factor_error());
+
+        // Write header and initial parameters into logger
+        log_header1(TERSE, "Compute error for source \""+m_srcname+"\""
+                           " parameter \""+m_model_par->name()+"\"");
+        log_value(NORMAL, "Confidence level",
+                  gammalib::str(m_confidence*100.0)+" %");
+        log_value(NORMAL, "Log-likelihood difference", m_dlogL);
+        log_value(NORMAL, "Initial factor range",
+                  "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
+
+        // Compute lower and upper boundaries
+        double value_lo = error_bisection(parmin, m_value);
+        double value_hi = error_bisection(m_value, parmax);
+
+        // Compute errors
+        double error           = 0.5 * (value_hi - value_lo);
+        double error_neg       = m_value  - value_lo;
+        double error_pos       = value_hi - m_value;
+        double error_value     = std::abs(error*m_model_par->scale());
+        double error_value_neg = std::abs(error_neg*m_model_par->scale());
+        double error_value_pos = std::abs(error_pos*m_model_par->scale());
+
+        // Write results into logger
+        std::string unit = " " + m_model_par->unit();
+        log_value(NORMAL, "Lower parameter factor", value_lo);
+        log_value(NORMAL, "Upper parameter factor", value_hi);
+        log_value(NORMAL, "Error from curvature",
+                  gammalib::str(m_model_par->error()) + unit);
+        log_value(NORMAL, "Error from profile",
+                  gammalib::str(error_value) + unit);
+        log_value(NORMAL, "Negative profile error",
+                  gammalib::str(error_value_neg) + unit);
+        log_value(NORMAL, "Positive profile error",
+                  gammalib::str(error_value_pos) + unit);
+
+        // Save error result
+        model->at(i).factor_error(error);
+
+    } // endfor: looped over spectral parameters
+
+    // Restore best fitting models (now with new errors computed)
+    m_obs.models(models_best);
     // Recover optimizer
     m_opt = best_opt;
@@ -35,7 +35,7 @@ outmodel, f, a, &quot;results.xml&quot;,,, &quot;Output model XML file&quot;
 #=============================
 confidence, r, h, 0.68,0.0,1.0, "Confidence level"
 tol,        r, h, 1e-3,,, "Computation tolerance"
-max_iter,   i, h, 50,,, "Maximum number of iterations"
+max_iter,   i, h, 50,1,1000, "Maximum number of iterations"
 #
 # Standard parameters
@@ -228,17 +228,9 @@ void ctulimit::run(void)
     // Get task parameters
     get_parameters();
-    // Set energy dispersion flag for all CTA observations and save old
+    // Set energy dispersion flags of all CTA observations and save old
     // values in save_edisp vector
-    std::vector<bool> save_edisp;
-    save_edisp.assign(m_obs.size(), false);
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            save_edisp[i] = obs->response()->apply_edisp();
-            obs->response()->apply_edisp(m_apply_edisp);
-        }
-    }
+    std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
     // Write input observation container into logger
     log_observations(NORMAL, m_obs, "Input observation");
@@ -338,13 +330,8 @@ void ctulimit::run(void)
     // Recover optimizer
     m_opt = best_opt;
-    // Restore energy dispersion flag for all CTA observations
-    for (int i = 0; i < m_obs.size(); ++i) {
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-        if (obs != NULL) {
-            obs->response()->apply_edisp(save_edisp[i]);
-        }
-    }
+    // Restore energy dispersion flags of all CTA observations
+    restore_edisp(m_obs, save_edisp);
     // Return
     return;
@@ -206,6 +206,7 @@ class Test(test):
         edispcube['yref']      = 22.01
         edispcube['migramax']  = 2.0
         edispcube['migrabins'] = 10
+        edispcube['addbounds'] = True
         edispcube['outcube']   = 'ctedispcube_py3.fits'
         edispcube['logfile']   = 'ctedispcube_py3.log'
         edispcube['chatter']   = 4
@@ -161,6 +161,41 @@ class Test(test):
         # Check that the cleared copy has also cleared the observations
         self._check_obs(cpy_error.obs(), nobs=0, nmodels=0)
+        # Prepare observation container with a single event list
+        obs = gammalib.GObservations()
+        obs.append(gammalib.GCTAObservation(self._events))
+        obs.models(gammalib.GModels(self._model))
+
+        # Setup cterror tool from observation container
+        error = ctools.cterror(obs)
+        error['srcname']  = 'Crab'
+        error['caldb']    = self._caldb
+        error['irf']      = self._irf
+        error['outmodel'] = 'cterror_py3.xml'
+        error['logfile']  = 'cterror_py3.log'
+        error['chatter']  = 4
+
+        # Execute cterror tool
+        error.logFileOpen()   # Make sure we get a log file
+        error.execute()
+
+        # Check result file
+        self._check_result_file('cterror_py3.xml')
+
+        # And now a run with not enough iterations
+        error['max_iter'] = 1
+        error['outmodel'] = 'cterror_py4.xml'
+        error['logfile']  = 'cterror_py4.log'
+        
+        # Execute cterror tool and catch the exception
+        error.logFileOpen()   # Make sure we get a log file
+        self.test_try('Test cterror with not enough iterations')
+        try:
+            error.execute()
+            self.test_try_failure('Exception not thrown')
+        except ValueError:
+            self.test_try_success()
+
         # Return
         return
@@ -160,6 +160,32 @@ class Test(test):
         # Check that the cleared copy has also cleared the observations
         self._check_obs(cpy_like.obs(), nobs=0, nmodels=0)
+        # Prepare observation container with a single event list
+        obs = gammalib.GObservations()
+        obs.append(gammalib.GCTAObservation(self._events))
+
+        # Prepare a special model to enforce TS computation
+        models = gammalib.GModels(self._model)
+        models['Crab'].tscalc(True)
+        obs.models(models)
+
+        # Allocate ctlike tool from observation container
+        like = ctools.ctlike(obs)
+        like['caldb']           = self._caldb
+        like['irf']             = self._irf
+        like['refit']           = True
+        like['fix_spat_for_ts'] = True
+        like['outmodel']        = 'ctlike_py3.xml'
+        like['logfile']         = 'ctlike_py3.log'
+        like['chatter']         = 4
+
+        # Execute ctlike
+        like.logFileOpen()  # Needed to get a new log file
+        like.execute()
+
+        # Check result file
+        self._check_result_file('ctlike_py3.xml')
+
         # Return
         return
@@ -177,6 +177,36 @@ class Test(test):
         # Check that the cleared copy has also cleared the map cube
         self._check_cube(cpy_mapcube.mapcube(), nmaps=0, npixels=0)
+        # Allocate ctmapcube tools, set models and do now a linear
+        # binning; also do not use Gaussian point sources
+        mapcube = ctools.ctmapcube()
+        mapcube.models(gammalib.GModels(self._model))
+        mapcube['ebinalg']  = 'LIN'
+        mapcube['emin']     = 0.1
+        mapcube['emax']     = 100.0
+        mapcube['enumbins'] = 10
+        mapcube['nxpix']    = 100
+        mapcube['nypix']    = 100
+        mapcube['binsz']    = 0.04
+        mapcube['coordsys'] = 'CEL'
+        mapcube['proj']     = 'CAR'
+        mapcube['xref']     = 83.63
+        mapcube['yref']     = 22.01
+        mapcube['ptsrcsig'] = 0.0
+        mapcube['outcube']  = 'ctmapcube_py3.fits'
+        mapcube['logfile']  = 'ctmapcube_py3.log'
+        mapcube['chatter']  = 4
+
+        # Execute mapcube
+        mapcube.logFileOpen()  # Needed to get a new log file
+        mapcube.execute()
+
+        # Check result file
+        self._check_result_file('ctmapcube_py3.fits')
+
+        # Publish with a name
+        mapcube.publish('My map cube')
+
         # Return
         return
@@ -189,11 +189,38 @@ class Test(test):
         # Check that the cleared copy has also cleared the model cube
         self._check_cube(cpy_model.cube(), nx=0, ny=0, nebins=0)
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Setup ctmodel tool from observation container and get model cube
+        # definition from counts cube
+        model = ctools.ctmodel(obs)
+        model.models(gammalib.GModels(self._model))
+        model['incube']   = self._cntcube
+        model['expcube']  = 'NONE'
+        model['psfcube']  = 'NONE'
+        model['bkgcube']  = 'NONE'
+        model['caldb']    = self._caldb
+        model['irf']      = self._irf
+        model['outcube']  = 'ctmodel_py3.fits'
+        model['logfile']  = 'ctmodel_py3.log'
+        model['chatter']  = 4
+
+        # Execute ctmodel tool
+        model.logFileOpen()  # Needed to get a new log file
+        model.execute()
+
+        # Check result file
+        self._check_result_file('ctmodel_py3.fits', nx=200, ny=200, nebins=20)
+
+        # Publish with name
+        model.publish('My model')
+
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, nx=40, ny=40, nebins=10):
         """
         Check result file
@@ -201,12 +228,18 @@ class Test(test):
         ----------
         filename : str
             Model cube file name
+        nx : int, optional
+            Number of pixels in X direction
+        ny : int, optional
+            Number of pixels in Y direction
+        nebins : int, optional
+            Number of energy bins
         """
         # Open result file
         cube = gammalib.GCTAEventCube(filename)
         # Check cube
-        self._check_cube(cube)
+        self._check_cube(cube, nx=nx, ny=ny, nebins=nebins)
         # Return
         return
@@ -154,6 +154,12 @@ class Test(test):
         sim['logfile']   = 'ctobssim_py2.log'
         sim['chatter']   = 3
+        # Double maximum event range
+        max_rate = sim.max_rate()
+        sim.max_rate(2*max_rate)
+        self.test_value(sim.max_rate(), 2*max_rate, 1.0e-3,
+                        'Check setting of maximum event rate')
+
         # Run tool
         sim.logFileOpen()
         sim.run()
@@ -187,27 +187,29 @@ class Test(test):
         # Get mixed observation container
         obs = self._obs_mixed()
-        # Set-up ctpsfcube from observation container
+        # Set-up ctpsfcube from observation container and input counts
+        # cube
         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
+        psfcube['incube']    = self._cntcube
+        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['addbounds'] = True
+        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
@@ -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
@@ -103,10 +104,25 @@ class Test(test):
         """
         Test ctskymap from Python
         """
-        # Set-up ctskymap
+        # Allocate ctskymap
         skymap = ctools.ctskymap()
+
+        # Check that empty ctskymap tool holds a map that has no pixels
+        self._check_map(skymap.map(), nx=0, ny=0)
+
+        # Check that saving does not nothing
+        skymap['outmap']  = 'ctskymap_py0.fits'
+        skymap['logfile'] = 'ctskymap_py0.log'
+        skymap.logFileOpen()
+        skymap.save()
+        self.test_assert(not os.path.isfile('ctskymap_py0.fits'),
+             'Check that no sky map has been created')
+
+        # Check that clearing does not lead to an exception or segfault
+        skymap.clear()
+
+        # Now set ctskymap parameters
         skymap['inobs']    = self._events
-        skymap['outmap']   = 'ctskymap_py1.fits'
         skymap['emin']     = 0.1
         skymap['emax']     = 100
         skymap['nxpix']    = 200
@@ -116,6 +132,7 @@ class Test(test):
         skymap['proj']     = 'CAR'
         skymap['xref']     = 83.63
         skymap['yref']     = 22.01
+        skymap['outmap']   = 'ctskymap_py1.fits'
         skymap['logfile']  = 'ctskymap_py1.log'
         skymap['chatter']  = 2
@@ -127,21 +144,111 @@ class Test(test):
         # Check result file
         self._check_result_file('ctskymap_py1.fits')
+        # Copy ctskymap tool
+        cpy_skymap = skymap.copy()
+
+        # Check sky map of ctskymap copy
+        self._check_map(cpy_skymap.map())
+
+        # Execute copy of ctskymap tool again, now with a higher chatter
+        # level than before
+        cpy_skymap['usepnt']  = True
+        cpy_skymap['emin']    = 0.2
+        cpy_skymap['emax']    = 150.0
+        cpy_skymap['outmap']  = 'ctskymap_py2.fits'
+        cpy_skymap['logfile'] = 'ctskymap_py2.log'
+        cpy_skymap['publish'] = True
+        cpy_skymap['chatter'] = 3
+        cpy_skymap.logFileOpen()  # Needed to get a new log file
+        cpy_skymap.execute()
+
+        # Check result file
+        self._check_result_file('ctskymap_py2.fits')
+
+        # Now clear copy of ctskymap tool
+        cpy_skymap.clear()
+
+        # Check that cleared ctskymap tool holds a map that has no pixels
+        self._check_map(cpy_skymap.map(), nx=0, ny=0)
+
+        # Get mixed observation container
+        obs = self._obs_mixed()
+
+        # Allocate ctskymap tool from observation container
+        skymap = ctools.ctskymap(obs)
+        skymap['emin']     = 0.1
+        skymap['emax']     = 100
+        skymap['nxpix']    = 200
+        skymap['nypix']    = 200
+        skymap['binsz']    = 0.02
+        skymap['coordsys'] = 'CEL'
+        skymap['proj']     = 'CAR'
+        skymap['xref']     = 83.63
+        skymap['yref']     = 22.01
+        skymap['outmap']   = 'ctskymap_py3.fits'
+        skymap['logfile']  = 'ctskymap_py3.log'
+        skymap['chatter']  = 4
+
+        # Execute tool
+        skymap.logFileOpen()  # Needed to get a new log file
+        skymap.execute()
+
+        # Check result file
+        self._check_result_file('ctskymap_py3.fits')
+
+        # Publish with name
+        skymap.publish('My sky map')
+
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, nx=200, ny=200):
         """
         Check result file
+
+        Parameters
+        ----------
+        filename : str
+            Sky map file name
+        nx : int, optional
+            Number of X pixels
+        ny : int, optional
+            Number of Y pixels
         """
         # Open result file
-        result = gammalib.GSkyMap(filename)
+        skymap = gammalib.GSkyMap(filename)
+        # Check sky map
+        self._check_map(skymap, nx=nx, ny=ny)
+
+        # Return
+        return
+
+    # Check sky map
+    def _check_map(self, skymap, nx=200, ny=200):
+        """
+        Check sky map
+
+        Parameters
+        ----------
+        skymap : `~gammalib.GSkyMap`
+            Sky map
+        nx : int, optional
+            Number of X pixels
+        ny : int, optional
+            Number of Y pixels
+        """
+        # Determine number of maps
+        if nx > 0 and ny > 0:
+            nmaps = 1
+        else:
+            nmaps = 0
+        
         # Check dimensions
-        self.test_value(result.nmaps(), 1, 'Check for one map')
-        self.test_value(result.nx(), 200, 'Check for 200 pixels in X')
-        self.test_value(result.ny(), 200, 'Check for 200 pixels in Y')
+        self.test_value(skymap.nmaps(), nmaps, 'Check number of maps')
+        self.test_value(skymap.nx(), nx, 'Check for number of X pixels')
+        self.test_value(skymap.ny(), ny, 'Check for number of Y pixels')
         # Return
         return
@@ -99,7 +99,25 @@ class Test(test):
         """
         Test ctulimit from Python
         """
-        # Set-up ctulimit
+        # Allocate ctulimit
+        ulimit = ctools.ctulimit()
+
+        # Check that empty ctulimit tool holds zero upper limits
+        self.test_value(ulimit.diff_ulimit(), 0.0, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(ulimit.flux_ulimit(), 0.0, 1.0e-17,
+                        'Check upper limit on photon flux')
+        self.test_value(ulimit.eflux_ulimit(), 0.0, 1.0e-16,
+                        'Check upper limit on energy flux')
+
+        # Check that saving does not nothing
+        ulimit.logFileOpen()
+        ulimit.save()
+
+        # Check that clearing does not lead to an exception or segfault
+        ulimit.clear()
+
+        # Now set ctulimit parameters
         ulimit = ctools.ctulimit()
         ulimit['inobs']   = self._events
         ulimit['inmodel'] = self._model
@@ -114,8 +132,50 @@ class Test(test):
         ulimit.run()
         ulimit.save()
-        # Check result file
-        self._check_result_file('ctulimit_py1.log')
+        # Check results
+        self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
+                        'Check upper limit on photon flux')
+        self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
+                        'Check upper limit on energy flux')
+
+        # Copy ctulimit tool
+        cpy_ulimit = ulimit.copy()
+
+        # Check results of copy
+        self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
+                        'ulimit upper limit on photon flux')
+        self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
+                        'Check upper limit on energy flux')
+
+        # Execute copy of ctulimit tool again, now with a higher chatter
+        # level than before
+        cpy_ulimit['logfile'] = 'ctulimit_py2.log'
+        cpy_ulimit['chatter'] = 3
+        cpy_ulimit.logFileOpen()  # Needed to get a new log file
+        cpy_ulimit.execute()
+
+        # Check results
+        self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
+                        'ulimit upper limit on photon flux')
+        self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
+                        'Check upper limit on energy flux')
+
+        # Now clear copy of ctulimit tool
+        cpy_ulimit.clear()
+
+        # Check that empty ctulimit tool holds zero upper limits
+        self.test_value(cpy_ulimit.diff_ulimit(), 0.0, 1.0e-23,
+                        'Check differential upper limit')
+        self.test_value(cpy_ulimit.flux_ulimit(), 0.0, 1.0e-17,
+                        'Check upper limit on photon flux')
+        self.test_value(cpy_ulimit.eflux_ulimit(), 0.0, 1.0e-16,
+                        'Check upper limit on energy flux')
         # Return
         return