Reducing sonar issues by adding more comments and cover more code with unit tests

Michael Mayer
1 parent 2b3d392a
Showing 6 changed files with 324 additions and 30 deletions
cscripts/csiactcopy.py
cscripts/csiactdata.py
cscripts/csiactobs.py
test/data/iactdata/master.json
test/test_csiactdata.py
test/test_csiactobs.py
@@ -91,7 +91,18 @@ class csiactcopy(ctools.cscript):
     def _copy(self, source, clobber):     
         """
         Copy file to outpath
+        
+        Parameters
+        ----------
+        source : string
+            Path of file to be copied
+
+        Returns
+        -------
+        filesize : float
+            Size of the file that was copied
         """
+        
         # Get file destination
         destination = os.path.join(self._outpath,
                                    os.path.relpath(source, self._remote_base))
@@ -155,10 +166,21 @@ class csiactcopy(ctools.cscript):
     def _merge(self, localfits, remotefits, hduname, clobber):
         """
         Merge remote and local fits files
-        
         If the local fits file is not present, a new one is created.
+        
+        Parameters
+        ----------
+        localfits : string
+            Path of local index FITS file         
+        remotefits : string
+            Path of remote index FITS file         
+        hduname : string
+            Name of HDU extension to be merged       
+        clobber : bool
+            Flag if remote content should overwrite local content        
         """    
+        # Logging
         if self._logTerse():
             self._log.parformat('Local file')
             self._log(str(localfits))
@@ -229,7 +251,7 @@ class csiactcopy(ctools.cscript):
                 tmp_hdu = remote[hduname]
                 # Get old size of remote HDU
-                size = remote_hdu.size()
+                size = remote_hdu.nrows()
                 # Extend remote HDU
                 remote_hdu.append_rows(tmp_hdu.nrows())
@@ -244,7 +266,8 @@ class csiactcopy(ctools.cscript):
         robs_col   = remote_hdu['OBS_ID']
         for obs_id in robs_col:
             remote_obs.append(obs_id)
-            
+        
+        # Log entries of index files   
         if self._logTerse():
             self._log.parformat('Remote entries')
             self._log(str(len(remote_obs)))
@@ -268,12 +291,18 @@ class csiactcopy(ctools.cscript):
                     # Remove local obs_id entries
                     table_has_obsid = True
+                    
+                    # Loop over obs ids and overwrite with remote content
                     while table_has_obsid:
+                        
+                        # First remove old rows
                         for i in range(local_hdu.nrows()):
                             if remote_obs_id == local_hdu['OBS_ID'][i]:
                                 local_hdu.remove_rows(i,1)
                                 removed_rows += 1
                                 break
+                        
+                        # Replace with remote content
                         table_has_obsid = False
                         for i in range(local_hdu.nrows()):
                             if remote_obs_id == local_hdu['OBS_ID'][i]:
@@ -309,6 +338,7 @@ class csiactcopy(ctools.cscript):
         # Add rmeote HDUs
         local_hdu.insert_rows(old_local_rows, remote_hdu.nrows())
+        # Log actions
         if self._logTerse():
             if clobber:
                 self._log.parformat('Removed local rows')
@@ -362,6 +392,8 @@ class csiactcopy(ctools.cscript):
         # Check if runlist file is available
         if self._runlist.exists():
             runfile = open(self._runlist.url())
+            
+            # Read runlist
             for line in runfile.readlines():
                 if len(line) == 0:
                     continue
@@ -369,6 +401,8 @@ class csiactcopy(ctools.cscript):
                     continue
                 if len(line.split()) > 0:
                     self._runs.append(int(line.split()[0]))
+            
+            # Close runlist file
             runfile.close()   
             # Logging
@@ -452,8 +486,8 @@ class csiactcopy(ctools.cscript):
                 # Initialise flag if SIZE column is present
                 has_size = table.contains('SIZE')
+                # Initialise file size
                 if has_size:
-                    # Initialise file size
                     cp_size = 0
                 # Initialise remote observation IDs
@@ -463,6 +497,8 @@ class csiactcopy(ctools.cscript):
                 # Log runs that are not available remotely
                 for run in self._runs:
+                    
+                    # Check for run not in remote data store
                     if not run in remote_ids:  
                         if self._logNormal():
                             self._log('Skip observation '+str(run)+': ')
@@ -526,6 +562,8 @@ class csiactcopy(ctools.cscript):
                         self._log.parformat('Size')
                         self._log('%.2f'%(float(cp_size)*1e-6)+' MB')
                         self._log('\n')
+                
+                # Log more information if requested
                 if self._logVerbose():
                     self._log('\n')
                     self._log.header3('File names')
@@ -556,20 +594,26 @@ class csiactcopy(ctools.cscript):
             self._log('\n')
             self._log.header1('Copying files')
-        # intialise counter
+        # Intialise counter
         k = 0
-        # initialise values for logging
+        # Initialise values for logging
         last_fraction      =  0.0
         fraction_increment = 20.0
+        
+        # Use 10% step increase
         if self._logNormal():
             fraction_increment = 10.0
+            
+        # Use 5% step increase
         if self._logTerse():
             fraction_increment = 5.0
+        
+        # Use 2% step increase
         if self._logExplicit():
             fraction_increment = 2.0
-        # initialise logging properties
+        # Initialise logging properties
         n_copied   = 0
         total_size = 0.0
@@ -579,20 +623,25 @@ class csiactcopy(ctools.cscript):
             # Log progress
             fraction = float(k) / float(len(files)) * 100.0
             while fraction > last_fraction:
+                
+                # Print status of copying procedure
                 if self._logNormal() and not self._logVerbose():
                     self._log.parformat('Status')
                     self._log(str(int(last_fraction))+'%')
                     self._log('\n')
+                
+                # Increment fraction for next logging step
                 last_fraction += fraction_increment
             # Copy file
             filesize = self._copy(filename, self._clobber())
+            # Increment counter if file size copied
             if filesize > 0.0:
                 total_size += filesize
                 n_copied   += 1
-            # Increment counter
+            # Increment file counter
             k += 1
         # Logging
@@ -600,14 +649,17 @@ class csiactcopy(ctools.cscript):
             self._log.parformat('Status')
             self._log('Finished')
             self._log('\n')
+        
+        # Logging about index files
         if self._logTerse():
             self._log('\n')
             self._log.header1('Updating index files')
-        # Build local index file names
+        # Build local hdu index file name
         local_hdu = os.path.join(self._outpath,
                                  os.path.relpath(remote_hdu,
                                                  self._remote_base))
+        # Build local obs index file name
         local_obs = os.path.join(self._outpath,
                                  os.path.relpath(remote_obs,
                                                  self._remote_base))
@@ -623,6 +675,7 @@ class csiactcopy(ctools.cscript):
             # Merge remote index files with local files
             self._merge(local_hdu, remote_hdu, 'HDU_INDEX', self._clobber())
+            # Logging
             if self._logTerse():
                 self._log('\n')
                 self._log.header3('OBS index')
@@ -634,7 +687,8 @@ class csiactcopy(ctools.cscript):
             # If all files were copied, just copy index files too
             self._copy(remote_hdu, self._clobber())
             self._copy(remote_obs, self._clobber())
-  
+        
+        # Logging
         if self._logTerse():
             self._log('\n')
             self._log.header3('Master index file')
@@ -95,6 +95,12 @@ class csiactdata(ctools.cscript):
     def _log_configuration(self, config):
         """
         Log configuration
+        
+        Parameters
+        ----------
+        config : dict
+            Dictionary containing data configuration
+
         """
         # Log information
         if self._logTerse():
@@ -111,8 +117,12 @@ class csiactdata(ctools.cscript):
             # Print additional information
             for key in config:
+                
+                # Skip keys that were treated above
                 if key in ['name','hduindx','obsindx']:
                     continue
+                
+                # Log residual keys
                 self._log.parformat(str(key))
                 self._log(str(config[key]))
                 self._log('\n')
@@ -174,27 +174,45 @@ class csiactobs(ctools.cscript):
         # Get HDUs
         for config in configs:
+            
+            # Check if prodname is present
             if self._prodname == config['name']:
                 self._hdu_index = str(os.path.join(self._datapath, config['hduindx']))
                 self._obs_index = str(os.path.join(self._datapath, config['obsindx']))
+                
+                # Leave loop if index file names were found
                 break
-        # Check HDUs
+        # Check index files
         if self._hdu_index == '' or self._obs_index == '':
             raise RuntimeError('*** ERROR: FITS data store "'+self._prodname+'" not available. Run csiactdata to get a list of available storage names')
+        
+        # Check HDU names
         filename = gammalib.GFilename(self._hdu_index+'[HDU_INDEX]')
         if not filename.is_fits():
             raise RuntimeError('*** ERROR: HDU index file "'+self._hdu_index+'[HDU_INDEX]" for FITS data store "'+self._prodname+'" not available. Check your master index file or run csiactdata to get a list of available storage names.')
         # Check for existence of 'BKG_SCALE' in the observation index file if required
         if self._use_bkg_scale:
+            
+            # Create filename
             filename = gammalib.GFilename(self._obs_index+'[OBS_INDEX]')
+            
+            # Check if it is a FITS file
             if filename.is_fits():
+                
+                # Open FITS file
                 fits = gammalib.GFits(self._obs_index)
+                
+                # Check if column "BKG_SCALE" is found and signal its possible usage
                 if not fits['OBS_INDEX'].contains('BKG_SCALE'):
                     self._use_bkg_scale = False
+                    
+                # Close FITS file
                 fits.close()
+                
             else:
+                # Signal that there is no background scale
                 self._use_bkg_scale = False
         # Create base data directory from hdu index file location
@@ -212,6 +230,11 @@ class csiactobs(ctools.cscript):
     def runlist(self, runlist):
         """
         Set observation list
+        
+        Parameters
+        ----------
+        runlist : list
+            List of observation IDs
         """
         # Check if convertable to a python list
@@ -228,14 +251,37 @@ class csiactobs(ctools.cscript):
         # Return
         return
-    def _background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0):
+    def _background_spectrum(self, prefactor, index, emin = 0.01, emax = 100.0):
+        """
+        Create a background spectrum model dependent on user parameters
+        
+        Parameters
+        ----------
+        prefactor : float
+            Prefactor parameter of spectral model
+        index : float
+            Index parameter of spectral model
+        emin : float
+            Minimum energy (in case a spectral node function is required)
+        emax : float
+            Maximum energy (in case a spectral node function is required)
+        
+        Returns
+        -------
+        spec : `~gammalib.GModelSpectral()`
+            Spectral model for the background shape     
+        """
         # Handle constant spectral model 
         if index == 0.0 and self._bkgpars <= 1:
             spec = gammalib.GModelSpectralConst()
+            
+            # Set parameter range
             spec['Normalization'].min(prefactor / self._bkg_range_factor)
             spec['Normalization'].max(prefactor * self._bkg_range_factor)
             spec['Normalization'].value(prefactor)
+            
+            # Freeze or release normalisation parameter
             if self._bkgpars == 0:
                 spec['Normalization'].fix()
             else:
@@ -245,8 +291,10 @@ class csiactobs(ctools.cscript):
             # Create power law model 
             if self._bkgpars <= 2:
-                e    = gammalib.GEnergy(1.0,'TeV')
-                spec = gammalib.GModelSpectralPlaw(prefactor, index, e)  
+                
+                # Set Power Law model with arbitrary pivot energy
+                pivot = gammalib.GEnergy(1.0,'TeV')
+                spec  = gammalib.GModelSpectralPlaw(prefactor, index, pivot)  
                 # Set parameter ranges
                 spec[0].min(prefactor / self._bkg_range_factor)
@@ -269,19 +317,31 @@ class csiactobs(ctools.cscript):
             else:
                 # Create reference powerlaw
-                plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,'TeV')) 
+                pivot = gammalib.GEnergy(1.0,'TeV')
+                plaw = gammalib.GModelSpectralPlaw(prefactor, index, pivot) 
                 # Create spectral model and energy values
                 spec = gammalib.GModelSpectralNodes()
+                
+                # Create logarithmic energy nodes
                 bounds = gammalib.GEbounds(self._bkgpars,gammalib.GEnergy(emin,'TeV'),gammalib.GEnergy(emax,'TeV'), True)
+                
+                # Loop over bounds and set intensity value
                 for i in range(bounds.size()):     
                     energy = bounds.elogmean(i)
-                    value = plaw.eval(energy, gammalib.GTime())
+                    value = plaw.eval(energy)
+                    
+                    # Append energy, value - tuple to node function
                     spec.append(energy, value)
+                
+                # Loop over parameters
                 for par in spec:
+                    # Fix energy nodes
                     if 'Energy' in par.name():
                         par.fix()
+                    
+                    # Release intensity nodes
                     elif 'Intensity' in par.name():   
                         value = par.value() 
                         par.scale(value)    
@@ -292,13 +352,36 @@ class csiactobs(ctools.cscript):
         return spec
     def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100):
+        """
+        Create an IACT background model
+        
+        Parameters
+        ----------
+        telescope : string
+            Name of telescope
+        obs_id : string
+            Observation ID
+        bkg_scale : float
+            Background scaling factor
+        bkgtype : string
+            Type of background (irf,aeff, or gauss)
+        Returns
+        -------
+        model : `~gammalib.GModelData()`
+            Background model for IACT observation
+        """
+              
         # handle IrfBackground
         if bkgtype == 'irf':
+            
+            # Set parameters to have a constant spectral model
             prefactor  = 1.0
             index      = 0.0
             prefactor *= bkg_scale
-            spec       = self._background_spectrum(obs_id, prefactor, index, emin, emax)
+            
+            # Create background spectrum
+            spec       = self._background_spectrum(prefactor, index, emin, emax)
             # Create background model instance
             bck = gammalib.GCTAModelIrfBackground(spec)
@@ -306,7 +389,7 @@ class csiactobs(ctools.cscript):
         # Set AeffBackground   
         elif bkgtype == 'aeff':
             prefactor = bkg_scale * self._bkg_aeff_norm
-            spec      = self._background_spectrum(obs_id, prefactor, self._bkg_aeff_index, emin, emax)
+            spec      = self._background_spectrum(prefactor, self._bkg_aeff_index, emin, emax)
             # Create background model instance
             bck = gammalib.GCTAModelAeffBackground(spec)
@@ -314,12 +397,14 @@ class csiactobs(ctools.cscript):
         # Set Gaussian Background
         elif bkgtype == 'gauss':
             prefactor = bkg_scale * self._bkg_gauss_norm
-            spec      = self._background_spectrum(obs_id, prefactor, self._bkg_gauss_index, emin, emax)
+            
+            # Create background spectrum
+            spec      = self._background_spectrum(prefactor, self._bkg_gauss_index, emin, emax)
             radial    = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma)
             bck       = gammalib.GCTAModelRadialAcceptance(radial, spec)
         else:
-            msg = 'Background type "'+self._bkgtype+'" unsupported'
+            msg = 'Background type "'+bkgtype+'" unsupported'
             raise RuntimeError(msg)
         # Copy model
@@ -383,6 +468,7 @@ class csiactobs(ctools.cscript):
             eventfile = aefffile = psffile = edispfile = bkgfile = ''
             eventhdu  = aeffhdu  = ''
+            # Initialise types and formats
             types = []
             formats = []
             for i in range(hduindx_hdu.nrows()):  
@@ -402,7 +488,11 @@ class csiactobs(ctools.cscript):
                                           hduindx_hdu['FILE_NAME'][index])
                 eventhdu   = hduindx_hdu['HDU_NAME'][index]
                 eventfile += '['+eventhdu+']'
+            
+            # Check if psf file is available
             if not gammalib.GFilename(eventfile).is_fits():
+                
+                # Print warning that observation will be skipped since no events are present
                 if self._logTerse():
                     self._log('Skipping observation '+str(obs_id)+
                               ': eventfile "'+eventfile+'" not found\n')
@@ -421,7 +511,11 @@ class csiactobs(ctools.cscript):
                                          hduindx_hdu['FILE_NAME'][index])
                 aeffhdu   = hduindx_hdu['HDU_NAME'][index]
                 aefffile += '['+aeffhdu+']'
+            
+            # Check if psf file is available
             if not gammalib.GFilename(aefffile).is_fits():
+                
+                # Print warning that observation will be skipped since no effective area is present
                 if self._logTerse():
                     self._log('Skipping observation '+str(obs_id)+
                               ': effective area "'+aefffile+'" not found\n')
@@ -439,7 +533,11 @@ class csiactobs(ctools.cscript):
                                         hduindx_hdu['FILE_DIR'][index],
                                         hduindx_hdu['FILE_NAME'][index])
                 psffile += '['+hduindx_hdu['HDU_NAME'][index]+']'
+            
+            # Check if psf file is available
             if not gammalib.GFilename(psffile).is_fits():
+                
+                # Print warning that observation will be skipped since no PSF is present
                 if self._logTerse():
                     self._log('Skipping observation '+str(obs_id)+
                               ': point spread function "'+psffile+
@@ -458,7 +556,11 @@ class csiactobs(ctools.cscript):
                                           hduindx_hdu['FILE_DIR'][index],
                                           hduindx_hdu['FILE_NAME'][index])
                 edispfile += '['+hduindx_hdu['HDU_NAME'][index]+']'
+            
+            # Check if edisp file is available
             if not gammalib.GFilename(edispfile).is_fits():
+                
+                # Print warning that edisp cannot be used for this observation
                 if self._logTerse():
                     self._log('Warning: observation '+str(obs_id)+
                               ' has no energy dispersion "'+edispfile+
@@ -480,17 +582,22 @@ class csiactobs(ctools.cscript):
                 bkgfile += '['+hduindx_hdu['HDU_NAME'][index]+']'
             if 'irf' in bkg_mod_hierarchy and not gammalib.GFilename(bkgfile).is_fits():
                 bkg_mod_hierarchy.remove('irf')
+                
+                # Throw warning if IRF background cannot be used
                 if self._logTerse():
                     self._log('Warning: observation '+str(obs_id)+
                               ' has no background information (file="'+
                               bkgfile+'"). IRF background cannot be used\n')
                     bkgfile = ''       
                     if len(bkg_mod_hierarchy) == 0:
+                        
+                        # Skip observation if no background can be used
                         if self._logTerse():
                             self._log('Skipping observation '+str(obs_id)+
                                       ': No background can be used\n')
                         continue
                     else:
+                        # Log if we fall back to next background approach
                         if self._logTerse():
                             self._log('Observation '+str(obs_id)+
                                       ': Falling back to background "'+
@@ -499,15 +606,21 @@ class csiactobs(ctools.cscript):
             # Close hdu index file
             hduindx.close()
-            # Handle background scale information if available
+            # Initialise background scale
             bkg_scale = 1.0
+            
+            # Check if background scale should be used
             if self._use_bkg_scale:
+                
+                # Read background scale from fits file
                 obsindx   = gammalib.GFits(self._obs_index+'[OBS_INDEX]'+obs_selection)
                 bkg_scale = obsindx['OBS_INDEX']['BKG_SCALE'][0]
                 obsindx.close()
-            # Open fits file to determine the observation name
+            # Open event FITS file
             fits        = gammalib.GFits(eventfile)
+            
+            # Get object and telescope strings from event hdu
             events      = fits[eventhdu]
             object_name = events.string('OBJECT')
             telescope   = events.string('TELESCOP')
@@ -550,6 +663,7 @@ class csiactobs(ctools.cscript):
                 self._log('Adding observation '+str(obs_id)+
                           ' ("'+object_name+'")\n')
+            # log more details
             if self._logExplicit():
                 self._log(' Event file: '+eventfile+'\n')
                 self._log(' Effective area: '+aefffile+'\n')
@@ -615,10 +729,13 @@ class csiactobs(ctools.cscript):
                 self._log('\n')
         else:
+            # Print warning if no observation could be be used from runlist
             self._log.header2('WARNING: No observation from given runlist available')
         # Append models provided by 'inmodels' if necessary
-        if not self._inmodels == None:
+        if self._inmodels.size() > 0:
+            
+            # Logging
             if self._logTerse():
                 self._log('\n')
                 self._log.header1('Appending models')
@@ -3,7 +3,14 @@
         {
             "name": "unit-test",
 	    "obsindx": "obs-index.fits.gz",
-            "hduindx": "hdu-index.fits.gz"
+            "hduindx": "hdu-index.fits.gz",
+	    "comment": "this is a comment"
+        },
+	{
+            "name": "unit-test-that-doesnt-exist",
+	    "obsindx": "badfile.fits.gz",
+            "hduindx": "unavailable_file.fits.gz",
+	    "comment": "this is a comment"
         }
     ]
 }
@@ -108,6 +108,13 @@ class Test(test):
         # Run csiactdata script and save run list
         iactdata.logFileOpen()   # Make sure we get a log file
         iactdata.run()
+        
+        # Get available configs
+        names = iactdata.names()
+        
+        # Check if there is at least one config available
+        self.test_assert(len(names) > 0,
+             'Check successful execution from Python')
         # Return
         return
@@ -77,7 +77,7 @@ class Test(test):
                         ' infile="runlist_cmd1.dat"'+ \
                         ' bkgpars=1'+\
                         ' outobs="obs_iactobs_cmd1.xml"'+ \
-                        ' outmodel="bgd_iactobs_cmd1.xml"'+ \
+                        ' outmodel="bkg_iactobs_cmd1.xml"'+ \
                         ' logfile="csiactobs_cmd1.log" chatter=1'
         # Check if execution of wrong command fails
@@ -92,7 +92,7 @@ class Test(test):
         self._check_obsdef('obs_iactobs_cmd1.xml')
         # Check model definition XML file
-        self._check_moddef('bgd_iactobs_cmd1.xml')
+        self._check_moddef('bkg_iactobs_cmd1.xml')
         # Setup csiactobs command
         cmd = csiactobs+' datapath="data_path_that_does_not_exist"'+ \
@@ -100,12 +100,25 @@ class Test(test):
                         ' infile="runlist_cmd1.dat"'+ \
                         ' bkgpars=1'+\
                         ' outobs="obs_iactobs_cmd2.xml"'+ \
-                        ' outmodel="bgd_iactobs_cmd2.xml"'+ \
+                        ' outmodel="bkg_iactobs_cmd2.xml"'+ \
                         ' logfile="csiactobs_cmd2.log" chatter=1'
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
-             'Check invalid input file when executed from command line')
+             'Check invalid input datapath when executed from command line')
+        
+        # Setup csiactobs command
+        cmd = csiactobs+' datapath="'+self._datapath+'"'+ \
+                        ' prodname="unit-test-doesnt-exist"'+ \
+                        ' infile="runlist_cmd1.dat"'+ \
+                        ' bkgpars=1'+\
+                        ' outobs="obs_iactobs_cmd3.xml"'+ \
+                        ' outmodel="bkg_iactobs_cmd3.xml"'+ \
+                        ' logfile="csiactobs_cmd3.log" chatter=1'       
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check invalid input prodname when executed from command line')
         # Return
         return
@@ -122,7 +135,7 @@ class Test(test):
         iactobs['infile']   = 'runlist_py1.dat'
         iactobs['bkgpars']  = 1
         iactobs['outobs']   = 'obs_iactobs_py1.xml'
-        iactobs['outmodel'] = 'bgd_iactobs_py1.xml'
+        iactobs['outmodel'] = 'bkg_iactobs_py1.xml'
         iactobs['logfile']  = 'csiactobs_py1.log'
         iactobs['chatter']  = 2
@@ -135,7 +148,93 @@ class Test(test):
         self._check_obsdef('obs_iactobs_py1.xml')
         # Check model definition XML file
-        self._check_moddef('bgd_iactobs_py1.xml')
+        self._check_moddef('bkg_iactobs_py1.xml')
+        
+        # Create test runlist
+        runlist = ['0']
+           
+        # Set-up csiactobs - part 2 (set runlist via function, 2 background pars)
+        iactobs = cscripts.csiactobs()
+        iactobs['datapath'] = self._datapath
+        iactobs['prodname'] = 'unit-test'
+        iactobs['bkgpars']  = 2
+        iactobs['outobs']   = 'obs_iactobs_py2.xml'
+        iactobs['outmodel'] = 'bkg_iactobs_py2.xml'
+        iactobs['logfile']  = 'csiactobs_py2.log'
+        iactobs['chatter']  = 3
+        iactobs.runlist(runlist)
+   
+        # Run csiactobs script and save run list
+        iactobs.logFileOpen()   # Make sure we get a log file
+        iactobs.run()
+        iactobs.save()
+           
+        # Test return functions
+        obs     = iactobs.obs()
+        ebounds = iactobs.ebounds() 
+           
+        self.test_assert(obs.size() == 1,
+                         'Check output observation container')
+        self.test_assert(ebounds.size() == 0,
+                         'Check output energy boundaries container')
+   
+        # Check observation definition XML file
+        self._check_obsdef('obs_iactobs_py2.xml')
+           
+        # Check model definition XML file
+        self._check_moddef('bkg_iactobs_py2.xml')
+           
+        # Set-up csiactobs - part 3 (use a large number of free parameters, aeff background)
+        iactobs = cscripts.csiactobs()
+        iactobs['datapath']      = self._datapath
+        iactobs['prodname']      = 'unit-test'
+        iactobs['infile']        = 'runlist_py1.dat'
+        iactobs['bkgpars']       = 8
+        iactobs['bkg_mod_hiera'] = 'aeff'
+        iactobs['outobs']        = 'obs_iactobs_py3.xml'
+        iactobs['outmodel']      = 'bkg_iactobs_py3.xml'
+        iactobs['logfile']       = 'csiactobs_py3.log'
+        iactobs['chatter']       = 4
+    
+        # Run csiactobs script and save run list
+        iactobs.logFileOpen()   # Make sure we get a log file
+        iactobs.run()   
+        iactobs.save()
+        
+        # Check observation definition XML file
+        self._check_obsdef('obs_iactobs_py3.xml')
+            
+        # Check model definition XML file
+        self._check_moddef('bkg_iactobs_py3.xml')
+            
+        # Set-up csiactobs - part 4 (use a Gauss background and inmodel parameter)
+        iactobs = cscripts.csiactobs()
+        iactobs['datapath']      = self._datapath
+        iactobs['inmodel']       = 'results.xml'
+        iactobs['prodname']      = 'unit-test'
+        iactobs['infile']        = 'runlist_py1.dat'
+        iactobs['bkgpars']       = 1
+        iactobs['bkg_mod_hiera'] = 'gauss'
+        iactobs['outobs']        = 'NONE'
+        iactobs['outmodel']      = 'NONE'
+        iactobs['logfile']       = 'csiactobs_py4.log'
+        iactobs['chatter']       = 2
+   
+        # Run csiactobs script
+        iactobs.logFileOpen()   # Make sure we get a log file
+        iactobs.run()
+        
+        # Get observations and models
+        observations      = iactobs.obs()
+        models            = observations.models()
+         
+        # Check number of models
+        self.test_assert(models.size() == 3,
+                       'Check for three models in container')
+
+        # Check number of observations
+        self.test_assert(observations.size() == 1,
+                       'Check for one observation in container')
         # Return
         return