Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents b8e06a5e 35e553ce
Showing 12 changed files with 238 additions and 133 deletions
cscripts/csiactobs.py
cscripts/cstsmapmerge.py
cscripts/cstsmapsplit.par
doc/source/users/reference_manual/cttsmap.rst
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.hpp
src/cttsmap/cttsmap.par
test/data/tsmap_0.fits
test/data/tsmap_1.fits
test/test_csiactobs.py
test/test_cstsmapmerge.py
test/test_cttsmap.py
@@ -85,6 +85,10 @@ class csiactobs(ctools.cscript):
         # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
+        # Append an observation list to XML instance
+        self._xml.append(gammalib.GXmlElement('observation_list title="observation list"'))
+        self._xml_obslist = self._xml.element('observation_list', 0)
+
         # Return
         return
@@ -559,12 +563,7 @@ class csiactobs(ctools.cscript):
         # Clear models
         self._models.clear()
-
-        # Clear XML file and append an observation list
-        self._xml.clear()
-        self._xml.append(gammalib.GXmlElement('observation_list title="observation list"'))
-        lib = self._xml.element('observation_list', 0)
-
+        
         # Log output
         msg = 'Looping over %d runs' % len(self._runlist)
         self._log_header1(gammalib.TERSE, msg)
@@ -713,7 +712,7 @@ class csiactobs(ctools.cscript):
                                                       run_emax.TeV()))
             # Append observation to XML and set attributes
-            obs = lib.append('observation')
+            obs = self._xml_obslist.append('observation')
             obs.attribute('name', object_name)
             obs.attribute('id', str(obs_id))
             obs.attribute('instrument', telescope)
@@ -768,7 +767,7 @@ class csiactobs(ctools.cscript):
         self._write_summary()
         # Write warning in no observation could be used from runlist
-        if lib.size() == 0:
+        if self._xml_obslist.size() == 0:
             self._log_string(gammalib.NORMAL, 'WARNING: No observation was '
                              'appended from specified runlist')
@@ -30,17 +30,17 @@ import ctools
 # ================== #
 class cstsmapmerge(ctools.cscript):
     """
-    Merge Test Statistic maps
+    Merge Test Statistic maps.
     """
     # Constructor
     def __init__(self, *argv):
         """
-        Constructor
+        Constructor.
         """
         # Set name
         self._name    = 'cstsmapmerge'
-        self._version = '1.2.0'
+        self._version = '1.1.0'
         # Initialise class members
         self._files        = None
@@ -64,7 +64,7 @@ class cstsmapmerge(ctools.cscript):
     # Private methods
     def _get_parameters(self):
         """
-        Get parameters from parfile
+        Get parameters from parfile.
         """
         # Get input maps string
         inmaps = self['inmaps'].string()
@@ -118,7 +118,7 @@ class cstsmapmerge(ctools.cscript):
         if self._read_ahead():
             self['outmap'].filename()
-        #  Write input parameters into logger
+        # Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
         # Return
@@ -126,12 +126,7 @@ class cstsmapmerge(ctools.cscript):
     def _init_ts_map(self, fitsfile):
         """
-        Initialise Test Statistic map
-
-        Parameters
-        ----------
-        fitsfile : str
-            FITS file to be merged
+        Initialise Test Statistic map.
         """
         # Set filename
         self._in_filename = fitsfile
@@ -169,12 +164,10 @@ class cstsmapmerge(ctools.cscript):
     def _merge_ts_map(self, fitsfile):
         """
-        Merge TS map from FITS file into output TS map
+        Merge TS map from FITS file into output TS map.
-        Parameters
-        ----------
-        fitsfile : str
-            FITS file to be merged
+        Args:
+            fitsfile: FITS file to be merged.
         """
         # Open FITS file
         fits = gammalib.GFits(fitsfile)
@@ -204,23 +197,23 @@ class cstsmapmerge(ctools.cscript):
         # Compare size of maps
         if not len(add_maps) == len(self._maps):
-            msg = ('Cannot merge map "%s" into map "%s" since the number of '
-                   'parameters between both maps is different.' %
-                   (fitsfile, self._in_filename))
+            msg = 'Cannot merge map "'+fitsfile+'" into map "'+\
+                  self._in_filename+'" since the number of parameters ' +\
+                  'between both maps is different.'
             raise RuntimeError(msg) 
         # Loop over bins    
         for i in range(self._tsmap.npix()):
             # Consider only bins that have been computed
-            if add_statusmap[i] > 0.5:
+            if add_statusmap[i] > -0.5:
                 # Raise exception if this bin has already been computed
-                if self._statusmap[i] > 0.5 and not self._overwrite:
-                    msg = ('Attempt to merge bin which apparently has already '
-                           'been merged. File "%s" contains already merged '
-                           'bins. Set hidden parameter "overwrite=yes" to '
-                           'avoid this error.' % fitsfile)
+                if self._statusmap[i] > -0.5 and not self._overwrite:
+                    msg = 'Attempt to merge bin which apparently has '+\
+                          'already been merged. File "'+fitsfile+'" '+\
+                          'contains already merged bins. Set hidden '+\
+                          'parameter "overwrite=yes" to avoid this error.'
                     raise RuntimeError(msg)
                 # Copy TS values
@@ -238,17 +231,13 @@ class cstsmapmerge(ctools.cscript):
     def _get_number_of_ts_pixels(self, fitsfile):
         """
-        Return number of pixels with TS values
+        Return number of pixels with TS values.
-        Parameters
-        ----------
-        fitsfile : str
-            FITS file to be merged
+        Args:
+            fitsfile: FITS file to be merged.
-        Returns
-        -------
-        counts : int
-            Number of pixels for which TS has been computed
+        Returns:
+            Number of pixels for which TS has been computed.
         """
         # Get status map for this file
         status = gammalib.GSkyMap(fitsfile+'[STATUS MAP]')
@@ -256,7 +245,7 @@ class cstsmapmerge(ctools.cscript):
         # Count number of pixels with TS status set
         count = 0
         for pix in status:
-            if pix > 0.5:
+            if pix > -0.5:
                 count += 1
         # Return number of pixels
@@ -266,7 +255,7 @@ class cstsmapmerge(ctools.cscript):
     # Public methods
     def run(self):
         """
-        Run the script
+        Run the script.
         """
         # Switch screen logging on in debug mode
         if self._logDebug():
@@ -276,9 +265,7 @@ class cstsmapmerge(ctools.cscript):
         self._get_parameters()
         # Write header into logger
-        if self._logTerse():
-            self._log('\n')
-            self._log.header1('Merge TS maps')
+        self._log_header1(gammalib.TERSE, 'Merge TS maps')
         # Initialise file to start with
         file0              = ''
@@ -290,8 +277,11 @@ class cstsmapmerge(ctools.cscript):
             # Skip file if it's not a FITS file
             if not gammalib.GFilename(fitsfile).is_fits():
-                self._log_value(gammalib.EXPLICIT, 'Skip file',
-                                '%s (not a FITS file)' % fitsfile)
+                
+                # Log message
+                self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)')
+
+                # Continue
                 continue
             # Open FITS file
@@ -299,29 +289,40 @@ class cstsmapmerge(ctools.cscript):
             # If file contains a status map then use it
             if fits.contains('STATUS MAP'):
+                
+                # Close FITS file
                 fits.close()
+                
+                # Set filename to initial filename
                 file0 = fitsfile
-                if self._logTerse():
-                    self._log(gammalib.parformat('Initial TS map file'))
-                    self._log(fitsfile)
-                    count = self._get_number_of_ts_pixels(fitsfile)
-                    self._log(' (%d TS pixels computed)' % count)
-                    self._log('\n')
+                
+                # Log message
+                pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, 'Initial TS map file', fitsfile + pix_info)
+                
+                # Leave loop
                 break
             # ... otherwise signal that file is useless
             else:
-                fits.close()
-                self._log_value(gammalib.EXPLICIT, 'Skip file',
-                                '%s (no "STATUS MAP" extension)' % fitsfile)
+                
+                # Close fits file
+                fits.close()  
+                
+                # Info message             
+                self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)')
+                
+                # Continue
                 continue
         # Signal if no suitable file was found
         if file0 == '':
-            if self._logTerse():
-                self._log('None of the provided files seems to be a sliced '
-                          'TS map file (none has a "STATUS MAP" extension).\n')
-        
+            msg = 'None of the provided files seems to be a sliced ' + \
+                          'TS map file (none has a "STATUS MAP" ' + \
+                          'extension)'
+            raise RuntimeError(msg)
+                          
+
         # ... otherwise merge files
         else:
@@ -342,9 +343,8 @@ class cstsmapmerge(ctools.cscript):
                 # Skip if file is not FITS
                 if not gammalib.GFilename(fitsfile).is_fits():
-                    self._log_value(gammalib.EXPLICIT, 'Skip file',
-                                    '%s (not a FITS file)' % fitsfile)
-                    continue
+                    self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)')
+                    continue   
                 # Open FITS file
                 fits = gammalib.GFits(fitsfile)
@@ -352,20 +352,15 @@ class cstsmapmerge(ctools.cscript):
                 # Skip if file does not contain status map
                 if not fits.contains('STATUS MAP'):
                     fits.close()
-                    self._log_value(gammalib.EXPLICIT, 'Skip file',
-                                    '%s (no "STATUS MAP" extension)' % fitsfile)
+                    self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)')
                     continue
                 # Close FITS file
                 fits.close()
                 # Logging
-                if self._logTerse():
-                    self._log(gammalib.parformat('Merge TS map file'))
-                    self._log(fitsfile)
-                    count = self._get_number_of_ts_pixels(fitsfile)
-                    self._log(' (%d TS pixels computed)' % count)
-                    self._log('\n')
+                pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, 'Merge TS map file', fitsfile + pix_info)
                 # Merge TS map
                 self._merge_ts_map(fitsfile)
@@ -378,16 +373,11 @@ class cstsmapmerge(ctools.cscript):
     def save(self):
         """ 
-        Save TS map and remove slices if requested
+        Save TS map and remove slices if requested.
         """
-        # Write header
-        self._log_header1(gammalib.TERSE, 'Save TS map')
         # Get output filename in case it was not read ahead
         outmap = self['outmap'].filename()
-
-        # Log filename
-        self._log_value(gammalib.TERSE, 'TS map file', outmap.url())
         # Create FITS file
         fits = gammalib.GFits()
@@ -403,33 +393,65 @@ class cstsmapmerge(ctools.cscript):
         for i in range(len(self._mapnames)):   
             fits[i+1].extname(self._mapnames[i])
+        # Treat status map separately 
+        self._statusmap.write(fits)
+        fits[fits.size()-1].extname('STATUS MAP')
+        
         # Check if map is fully done
-        done = True
+        bins_merged = 0
         for pix in self._statusmap:
-            if pix < 0.5:
-                done = False
-                break
+            if pix > -0.5:
+                bins_merged += 1
+
+        # Set boolean if map is fully done
+        done = (bins_merged == self._tsmap.npix())
-        # Write status map if we are not done yet
+        # Log summary header
+        self._log_header2(gammalib.TERSE, 'Merging Summary')
+                
+        # Write computing status if we are not done yet
         if not done:
-            self._statusmap.write(fits)
-            fits[fits.size()-1].extname('STATUS MAP')
+            
+            # Log bins that were computed and those that were missing
+            self._log_value(gammalib.TERSE, 'TS map bins', self._tsmap.npix())
+            self._log_value(gammalib.TERSE, 'Bins computed and merged', bins_merged )
+            self._log_value(gammalib.TERSE, 'Bins missing', self._tsmap.npix() - bins_merged)
+        
+        else:
+            
+            # Write success message
+            self._log_string(gammalib.TERSE, 'TS map was fully computed and successfully merged') 
+        
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save TS map')
+        
+        # Log filename
+        self._log_value(gammalib.TERSE, 'TS map file', outmap.url())
         # Save FITS file
-        fits.saveto(outmap, self._clobber())
+        fits.saveto(outmap, self._clobber()) 
-        # Delete TS input maps if requested
-        if self._delete:
+        # Delete TS input maps if requested and map is fully done
+        if self._delete and done:
+            
+            # Log header
+            self._log_header2(gammalib.TERSE, 'Delete slices TS map files')
+            
+            # Loop over sliced files
             for filename in self._merged_files:
+                
+                # Remove file
                 os.remove(filename)
-                self._log_value(gammalib.TERSE, 'Deleted input file', filename)
+                
+                # Log message about file
+                self._log_value(gammalib.TERSE,'Deleted input file', filename)
         # Return
         return
     def execute(self):
         """
-        Execute the script
+        Execute the script.
         """
         # Open logfile
         self.logFileOpen()
@@ -27,6 +27,7 @@ bkgcube,       f, a, NONE,,, &quot;Input background cube file (only needed for stacke
 caldb,         s, a, prod2,,, "Calibration database"
 irf,           s, a, South_0.5h,,, "Instrument response function"
 edisp,         b, h, no,,, "Apply energy dispersion?"
+errors,        b, h, no,,, "Compute/store parameter errors?"
 outmap,        f, a, tsmap.fits,,, "Output Test Statistic map file"
 bins_per_job,  i, a, 5,,, "Number of TS map bins per task"
 compute_null,  b, a, yes,,,"Compute null hypothesis first?"
@@ -21,7 +21,10 @@ Gaussian sigma).
 cttsmap generates a FITS file comprising a sky map of TS values followed by 
 one extension per free parameter that contains sky maps of the fitted 
-parameter values.
+parameter values. To save computation time, the errors are not computed by
+default. Specifying the hidden parameter ``errors=yes`` will switch on the
+error computation. Accordingly, for each free parameter, the errors are also
+stored in a separate sky map.
 General parameters
@@ -56,6 +59,9 @@ General parameters
 ``outmap [file]``
     Output Test Statistic map file.
+    
+``(errors = no) [boolean]``
+    Compute and store parameter errors?
 ``(usepnt = no) [boolean]``
     Use CTA pointing direction for map centre instead of xref/yref parameters?
@@ -173,9 +173,7 @@ cttsmap&amp; cttsmap::operator=(const cttsmap&amp; app)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear cttsmap tool
- *
- * Clears cttsmap tool.
+ * @brief Clear instance
  ***************************************************************************/
 void cttsmap::clear(void)
 {
@@ -194,9 +192,6 @@ void cttsmap::clear(void)
     this->ctlikelihood::init_members();
     init_members();
-    // Write header into logger
-    log_header();
-
     // Return
     return;
 }
@@ -290,6 +285,14 @@ void cttsmap::run(void)
     	// Optimize observation container
     	m_obs.optimize(m_opt);
+    	// Compute errors if necessary
+    	if (m_errors) {
+    		m_obs.errors(m_opt);
+    	}
+
+    	// Get status of optimization
+    	int status = m_opt.status();
+
     	// Retrieve the Likelihood value
     	double logL1 = -(m_opt.value());
@@ -314,11 +317,27 @@ void cttsmap::run(void)
     	// Extract fitted test source parameters
         for (int j = 0; j < m_mapnames.size(); ++j) {
-            m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
-        }
-    	// Set status of bin to true
-    	m_statusmap(i) = 1.0;
+            // If map name start with "e_" then set the parameter error ...
+            if (m_mapnames[j].substr(0,2) == "e_") {
+
+                // Get parameter name by removing the error prefix
+            	std::string parname = m_mapnames[j].substr(2, m_mapnames[j].size());
+
+                // Get parameter error
+           		m_maps[j](i) = (*testsource)[parname].error();
+
+            }
+
+            // ... otherwise set the parameter value
+            else {
+                m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
+            }
+
+        } // endfor: looped over all maps
+
+        // Set Fit status of the bin
+        m_statusmap(i) = status;
     	// Remove model from container
     	models.remove(m_srcname);
@@ -371,11 +390,9 @@ void cttsmap::save(void)
             fits[i+1]->extname(m_mapnames[i]);
         }
-        // Add computation log if not all bins are computed
-        if (m_binmin != -1 || m_binmax != -1) {
-            m_statusmap.write(fits);
-            fits[m_mapnames.size()+1]->extname("STATUS MAP");
-        }
+        // Write Fit status map
+        m_statusmap.write(fits);
+        fits[fits.size()-1]->extname("STATUS MAP");
         // Save FITS file
         fits.saveto(m_outmap, clobber());
@@ -437,6 +454,7 @@ void cttsmap::init_members(void)
     m_outmap.clear();
     m_apply_edisp = false;
     m_publish     = false;
+    m_errors      = false;
     // Initialise protected members
     m_binmin     = -1;
@@ -465,6 +483,7 @@ void cttsmap::copy_members(const cttsmap&amp; app)
     m_outmap      = app.m_outmap;
     m_apply_edisp = app.m_apply_edisp;
     m_publish     = app.m_publish;
+    m_errors      = app.m_errors;
     // Copy protected members
     m_binmin    = app.m_binmin;
@@ -516,7 +535,7 @@ void cttsmap::get_parameters(void)
 {
     // Setup observations from "inobs" parameter
     setup_observations(m_obs);
-    
+
     // Setup models from "inmodel" parameter
     setup_models(m_obs, (*this)["srcname"].string());
@@ -545,6 +564,9 @@ void cttsmap::get_parameters(void)
     // Create sky map based on task parameters
     GSkyMap map = create_map(m_obs);
+    // Check if errors should be computed
+    m_errors  = (*this)["errors"].boolean();
+
     // Initialise maps from user parameters
     init_maps(map);
@@ -585,14 +607,17 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
     m_tsmap = GSkyMap(map);
     // Initialise map information
-    m_statusmap.clear();
+	m_statusmap.clear();
+
+	// Create status map and progress map
+	m_statusmap = GSkyMap(map);
-    // Create status map
-    m_statusmap = GSkyMap(map);
+	// Set status map to invalid value -1 to signal that bin wasn't computed yet
+	m_statusmap = -1.0;
     // Loop over all model parameters
     for (int i = 0; i < m_testsource->size(); i++) {
-    
+
         // Skip fixed parameters
         if ((*m_testsource)[i].is_fixed()) {
             continue;
@@ -604,7 +629,18 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
         // Store parameter name
         m_mapnames.push_back((*m_testsource)[i].name());
-    
+
+        // Add parameter errors if requested
+        if (m_errors) {
+
+            // Add sky map for fit error of the free parameter
+            m_maps.push_back(map);
+
+            // Store name of parameter error
+            m_mapnames.push_back("e_" + (*m_testsource)[i].name());
+
+        } // endif: error calculation was requested
+        
     } // endfor: looped over all model parameters
     // Return
@@ -82,19 +82,20 @@ protected:
     void init_maps(const GSkyMap& map);
     // User parameters
-    std::string              m_srcname;     //!< Name of source which is moved around
-    GFilename                m_outmap;      //!< Output TS map file name
-    bool                     m_apply_edisp; //!< Apply energy dispersion?
-    bool                     m_publish;     //!< Publish TS map?
+    std::string m_srcname;     //!< Name of source which is moved around
+    GFilename   m_outmap;      //!< Output TS map file name
+    bool        m_apply_edisp; //!< Apply energy dispersion?
+    bool        m_publish;     //!< Publish TS map?
+    bool        m_errors;      //!< Compute and store parameter errors?
     // Parameters to control speed and job splitting
-    int                      m_binmin;      //!< Map bin number from which computation should start
-    int                      m_binmax;      //!< Map bin number where map computation should end
-    double                   m_logL0;       //!< Likelihood value of null hypothesis
+    int         m_binmin;      //!< Map bin number from which computation should start
+    int         m_binmax;      //!< Map bin number where map computation should end
+    double      m_logL0;       //!< Likelihood value of null hypothesis
     // Protected members
     GSkyMap                  m_tsmap;       //!< TS map
-    GSkyMap                  m_statusmap;   //!< Map of computed bins
+    GSkyMap                  m_statusmap;   //!< Map of optimizer fit status
     std::vector<std::string> m_mapnames;    //!< Names of free parameters
     std::vector<GSkyMap>     m_maps;        //!< Sky maps for each free parameter
     GModel*                  m_testsource;  //!< Pointer to test source for TS computation
@@ -5,11 +5,11 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      cttsmap.par                                                   #
-#  Version:   1.1.0                                                         #
+#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Michael Mayer                                                 #
-#             michael.mayer@desy.de                                         #
-#             DESY                                                          #
+#             michael.mayer@physik.hu-berlin.de                             #
+#             HU Berlin                                                     #
 #                                                                           #
 #  Purpose:   Parameter file for the cttsmap application                    #
 #                                                                           #
@@ -30,6 +30,11 @@ edisp,     b, h, no,,, &quot;Apply energy dispersion?&quot;
 outmap,    f, a, "tsmap.fits",,, "Output Test Statistic map file"
 #
+# TS map computation parameters
+#==============================
+errors,    b, h, no,,, "Compute and store parameter errors?"
+
+#
 # Spatial binning parameters
 #===========================
 usepnt,    b, h, no,,, "Use pointing instead of xref/yref parameters?"
@@ -43,7 +48,7 @@ proj,      s, a, &quot;CAR&quot;, AIT|AZP|CAR|MER|MOL|STG|TAN,, &quot;Projection method&quot;
 #
 # Parameters for splitting and speed purpose
-#===========================
+#===========================================
 binmin,    i, h, -1,,, "First bin to compute"
 binmax,    i, h, -1,,, "Last bin to compute"
 logL0,     r, h, 0.0,,, "LogLikelihood value of null hypothesis"
@@ -148,7 +148,7 @@ class Test(test):
         # Check empty observation definition XML file
         self._check_obsdef('csiactobs_obs_py0.xml', 0)
-        
+
         # Check empty model definition XML file
         self._check_moddef('csiactobs_bgd_py0.xml', 0)
@@ -98,7 +98,7 @@ class Test(test):
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
              'Check invalid input file when executed from command line')
-
+        
         # Return
         return
@@ -160,10 +160,10 @@ class Test(test):
         # Check merged TS map
         self._check_result_file('cstsmapmerge_py4.fits')
-
-        # Set-up space-separated cstsmapmerge with invalid files
+        
+        # Set-up cstsmapmerge with incomplete file list
         merge = cscripts.cstsmapmerge()
-        merge['inmaps']  = self._events+' '+self._model
+        merge['inmaps']  = self._map0 + ' ' + self._events
         merge['outmap']  = 'cstsmapmerge_py5.fits'
         merge['logfile'] = 'cstsmapmerge_py5.log'
         merge['chatter'] = 4
@@ -172,13 +172,13 @@ class Test(test):
         merge.execute()
         # Check merged TS map
-        #self._check_result_file('cstsmapmerge_py5.fits')
+        self._check_result_file('cstsmapmerge_py5.fits', False)
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, test_complete = True):
         """
         Check result file
         """
@@ -189,6 +189,7 @@ class Test(test):
         ts        = fits['Primary']
         prefactor = fits['Prefactor']
         index     = fits['Index']
+        status    = fits['STATUS MAP']
         # Check dimensions
         self.test_value(ts.naxis(), 2, 'Check for 2 dimensions')
@@ -200,6 +201,23 @@ class Test(test):
         self.test_value(index.naxis(), 2, 'Check for 2 dimensions')
         self.test_value(index.naxes(0), 2, 'Check for 2 pixels in X')
         self.test_value(index.naxes(1), 1, 'Check for 1 pixel in Y')
+        self.test_value(status.naxis(), 2, 'Check for 2 dimensions')
+        self.test_value(status.naxes(0), 2, 'Check for 2 pixels in X')
+        self.test_value(status.naxes(1), 1, 'Check for 1 pixel in Y')
+        
+        # Initialise flag if map is complete
+        done = True
+        
+        # Loop over status sky map
+        for pix in status:
+            
+            # Check if pix is larger than threshold
+            if pix < -0.5:
+                done = False
+                break
+        
+        # Test for completeness of merged map
+        self.test_assert(done == test_complete, 'Test if map was merged completely')           
         # Return
         return
@@ -98,6 +98,23 @@ class Test(test):
         self.test_assert(self._execute(cmd) != 0,
              'Check invalid input file when executed from command line')
+        # Setup cttsmap command
+        cmd = cttsmap+' inobs="'+self._events+'"'+ \
+                      ' inmodel="'+self._model+'" srcname="Crab"'+ \
+                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                      ' errors=yes'+ \
+                      ' outmap="cttsmap_cmd3.fits"'+ \
+                      ' nxpix=3 nypix=3 binsz=0.02'+ \
+                      ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
+                      ' logfile="cttsmap_cmd3.log" chatter=1'
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check result file
+        self._check_result_file('cttsmap_cmd3.fits')
+
         # Setup cttsmap --help option
         cmd = cttsmap+' --help'
@@ -267,7 +284,7 @@ class Test(test):
             nmaps = 1
         else:
             nmaps = 0
-        
+
         # Check dimensions
         self.test_value(tsmap.nmaps(), nmaps, 'Check number of maps')
         self.test_value(tsmap.nx(), nx, 'Check for number of X pixels')