Added a more realistic and larger IACT test database and added various IACT test…

… cases, restructured csiactobs to avoid code duplication and added csobsdef::obs() function to retrieve observation container

Added a more realistic and larger IACT test database and added various IACT test…
… cases, restructured csiactobs to avoid code duplication and added csobsdef::obs() function to retrieve observation container
Michael Mayer
1 parent ff0d52ee
Showing 24 changed files with 427 additions and 224 deletions
cscripts/Makefile.am
cscripts/csiactcopy.py
cscripts/csiactobs.py
cscripts/csobsdef.py
test/data/iact_runlist.dat
test/data/iact_runlist_1.dat
test/data/iact_runlist_2.dat
test/data/iactdata/events.fits.gz
test/data/iactdata/events/events_0.fits.gz
test/data/iactdata/events/events_1.fits.gz
test/data/iactdata/events/events_2.fits.gz
test/data/iactdata/events/events_3.fits.gz
test/data/iactdata/events/events_4.fits.gz
test/data/iactdata/events/events_5.fits.gz
test/data/iactdata/events/events_6.fits.gz
test/data/iactdata/events/events_7.fits.gz
test/data/iactdata/events/events_8.fits.gz
test/data/iactdata/irf_file.fits.gz → test/data/iactdata/events/irf_file.fits.gz
test/data/iactdata/hdu-index.fits.gz
test/data/iactdata/master.json
@@ -155,14 +155,24 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/model_filefunction.txt \
            $(top_srcdir)/test/data/model_map.fits.gz \
            $(top_srcdir)/test/data/model_mapcube.fits.gz \
-           $(top_srcdir)/test/data/iact_runlist.dat
+           $(top_srcdir)/test/data/iact_runlist.dat \
+           $(top_srcdir)/test/data/iact_runlist_1.dat \
+           $(top_srcdir)/test/data/iact_runlist_1.dat
  
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
                $(top_srcdir)/test/data/iactdata/master.json \
                $(top_srcdir)/test/data/iactdata/obs-index.fits.gz \
-               $(top_srcdir)/test/data/iactdata/irf_file.fits.gz \
-               $(top_srcdir)/test/data/iactdata/events.fits.gz
+               $(top_srcdir)/test/data/iactdata/events/irf_file.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_0.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_1.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_2.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_3.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_4.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_5.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_6.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_7.fits.gz \
+               $(top_srcdir)/test/data/iactdata/events/events_8.fits.gz              
  
 # Install scripts, parfiles, test scripts and test data. Note that we don't
 # need to prepend $(DESTDIR) here (not fully clear why, but that's the way
@@ -377,9 +377,7 @@ class csiactcopy(ctools.cscript):
  
         # Copy all data if runlist file is 'NONE'
         elif self._runlist == 'NONE':
-            if self._logTerse():
-                self._log('Copy all available data')
-                self._log('\n')
+            self._log_string(gammalib.TERSE, 'Copy all available data')
  
         # Raise exception if file not valid
         else:
@@ -438,6 +438,94 @@ class csiactobs(ctools.cscript):
         # Return
         return
  
+    def _get_filename(self, hdu, hierarchy, formats): 
+        """
+        Retrieves a filename from a hdu index
+        
+        Parameters
+        ----------
+        hdu : `gammalib.GFitsTable`
+            HDU FITS table
+        hierarchy : list of strings
+            List of strings containing the hierarchy how to look for files
+        formats : list of strings
+            File formats available for this observation
+        
+        Returns
+        -------
+        filename : string
+            Filename of requested file
+        """
+        
+        # Initialise filename
+        filename = ''
+        
+        # Handle hierarchies
+        index = -1
+        
+        # Loop over possible formats
+        for version in hierarchy:
+            
+            # Count number of same available  versions
+            n = formats.count(version) 
+            
+            # Leave loop and store index of first occurence
+            if n > 0:
+                index = formats.index(version)
+                break
+        
+        # Build filename from information stored in hduindx_hdu
+        if index >= 0:
+            filename  = os.path.join(self._subdir,
+                                      hdu['FILE_DIR'][index],
+                                      hdu['FILE_NAME'][index])
+            filehdu   = hdu['HDU_NAME'][index]
+            filename += '['+filehdu+']'
+
+        # Return filename
+        return filename
+    
+    
+    def _is_present(self, filename, obs_id, filetype):
+        """
+        Checks if a filename is present
+        
+        Parameters
+        ----------
+        filename : string
+            Filename
+        obs_id : int
+            Observation ID
+        filetype : string
+            Type of file
+        
+        Returns
+        -------
+        present : bool
+            Flag signaling if file is present
+        """
+                
+        # Initialise return value
+        present = True
+        
+        # Create filename instance
+        fname = gammalib.GFilename(filename)
+        
+        # Check if file is empty
+        if fname.is_empty():
+            self._log_string(gammalib.TERSE, 'Skipping observation '+str(obs_id)+
+                          ': No '+filetype+' found\n')
+            present = False
+        
+        # Check if file exists
+        elif not fname.exists():
+            self._log_string(gammalib.TERSE, 'Skipping observation '+str(obs_id)+
+                          ': '+filetype+' "'+filename+'" does not exist\n')
+            present = False
+        
+        # Return boolean
+        return present    
+
  
     # Public methods
     def run(self):
@@ -477,149 +565,77 @@ class csiactobs(ctools.cscript):
             hduindx     = gammalib.GFits(self._hdu_index+'[HDU_INDEX]'+obs_selection)
             hduindx_hdu = hduindx['HDU_INDEX']
  
-            # Initialise files and hdu names
-            eventfile = aefffile = psffile = edispfile = bkgfile = ''
-            eventhdu  = aeffhdu  = ''
-            
             # Initialise types and formats
-            types = []
             formats = []
             for i in range(hduindx_hdu.nrows()):  
-                types.append(hduindx_hdu['HDU_TYPE'][i])
                 formats.append(hduindx_hdu['HDU_CLASS'][i])
  
-            # Handle events
-            index = -1
-            for version in self._ev_hiera:
-                n = formats.count(version) 
-                if n > 0:
-                    index = formats.index(version)
-                    break
-            if index >= 0:
-                eventfile  = os.path.join(self._subdir,
-                                          hduindx_hdu['FILE_DIR'][index],
-                                          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')
+            # Retrieve filenames            
+            eventfile = self._get_filename(hduindx_hdu, self._ev_hiera, formats)
+            aefffile  = self._get_filename(hduindx_hdu, self._aeff_hiera, formats)
+            psffile   = self._get_filename(hduindx_hdu, self._psf_hiera, formats)
+            edispfile = self._get_filename(hduindx_hdu, self._edisp_hiera, formats)
+            bkgfile   = self._get_filename(hduindx_hdu, self._bkg_hiera, formats)
+            
+            # Check for presence of required event file
+            if not self._is_present(eventfile, obs_id, "event file"):
                 continue
  
-            # Handle Aeff
-            index = -1
-            for version in self._aeff_hiera:
-                n = formats.count(version) 
-                if n > 0:
-                    index = formats.index(version)
-                    break
-            if index >= 0:
-                aefffile  = os.path.join(self._subdir,
-                                         hduindx_hdu['FILE_DIR'][index],
-                                         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')
+            # Check for presence of required effective area file
+            if not self._is_present(aefffile, obs_id, "effective area file"):
                 continue
-                                
-            # Handle psf
-            index = -1
-            for version in self._psf_hiera:
-                n = formats.count(version) 
-                if n > 0:
-                    index = formats.index(version)
-                    break
-            if index >= 0:
-                psffile  = os.path.join(self._subdir,
-                                        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+
-                              '" not found\n')
+            
+            # Check for presence of required psf file
+            if not self._is_present(psffile, obs_id, "PSF file"):
                 continue
-                
-            # Handle edisp
-            index = -1
-            for version in self._edisp_hiera:
-                n = formats.count(version) 
-                if n > 0:
-                    index = formats.index(version)
-                    break
-            if index >= 0:
-                edispfile  = os.path.join(self._subdir,
-                                          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():
+            
+            # Check for optional energy dispersion file
+            if not gammalib.GFilename(edispfile).exists():
  
                 # 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+
-                              '" information\n')
-                    edispfile = ''
+                self._log_string(gammalib.TERSE, 'Warning: observation '+str(obs_id)+
+                              ' has no energy dispersion information\n')
  
-            # Handle background
-            index = -1
-            bkg_mod_hierarchy = list(self._bkg_mod_hiera)
-            for version in self._bkg_hiera:
-                n = formats.count(version) 
-                if n > 0:
-                    index = formats.index(version)
-                    break
-            if index >= 0:
-                bkgfile  = os.path.join(self._subdir,
-                                        hduindx_hdu['FILE_DIR'][index],
-                                        hduindx_hdu['FILE_NAME'][index])
-                bkgfile += '['+hduindx_hdu['HDU_NAME'][index]+']'
-            if 'irf' in bkg_mod_hierarchy and not gammalib.GFilename(bkgfile).is_fits():
-                bkg_mod_hierarchy.remove('irf')
+                # Set energy dispersion file as empty
+                edispfile = ''
+            
+            # Create a copy of background model hierarchy for this run
+            run_bkg_mod_hierarchy = list(self._bkg_mod_hiera)
+            
+            # Check if background file is available
+            # Remove IRF background from hierarchy if file doesnt exist
+            if not gammalib.GFilename(bkgfile).exists():
  
-                # 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 "'+
-                                      bkg_mod_hierarchy[0]+'"\n')
+                # Set background file as empty
+                bkgfile = ''
+                
+                # Check for IRF background in background model hierarchy
+                if 'irf' in run_bkg_mod_hierarchy:
+                    
+                    # Remove IRF background if file doesnt exist
+                    run_bkg_mod_hierarchy.remove('irf')
+                
+                    # Print warning that edisp cannot be used for this observation
+                    self._log_string(gammalib.TERSE, 'Warning: observation '+str(obs_id)+
+                                  ' has no background information (IRF background cannot be used)\n')
+                
+                # Check if background model hierarchy is empty
+                if len(run_bkg_mod_hierarchy) == 0:
+                    
+                    # Skip observation if no background can be used
+                    self._log_string(gammalib.TERSE, 'Skipping observation '+str(obs_id)+
+                                  ': No background can be used\n')
+                    continue
+                
+                else:
+                    # Log if we fall back to next background approach
+                    self._log_string(gammalib.TERSE, 'Observation '+str(obs_id)+
+                                  ': Falling back to background "'+
+                                  run_bkg_mod_hierarchy[0]+'"\n')
  
             # Close hdu index file
             hduindx.close()
-
+            
             # Initialise background scale
             bkg_scale = 1.0
  
@@ -635,7 +651,7 @@ class csiactobs(ctools.cscript):
             fits        = gammalib.GFits(eventfile)
  
             # Get object and telescope strings from event hdu
-            events      = fits[eventhdu]
+            events      = fits[gammalib.GFilename(eventfile).extname()]
             object_name = events.string('OBJECT')
             telescope   = events.string('TELESCOP')
  
@@ -644,7 +660,7 @@ class csiactobs(ctools.cscript):
  
             # Open effective area file to look for threshold
             aeff_fits  = gammalib.GFits(aefffile)
-            aeff_table = aeff_fits[aeffhdu]
+            aeff_table = aeff_fits[gammalib.GFilename(aefffile).extname()]
  
             # Set energy range from header keyword if present
             if aeff_table.has_card('LO_THRES') and aeff_table.has_card('HI_THRES'):
@@ -668,7 +684,7 @@ class csiactobs(ctools.cscript):
             self._models.append(self._iact_background(telescope,
                                                       obs_id,
                                                       bkg_scale,
-                                                      bkg_mod_hierarchy[0],
+                                                      run_bkg_mod_hierarchy[0],
                                                       run_emin.TeV(),
                                                       run_emax.TeV()))
  
@@ -684,19 +700,19 @@ class csiactobs(ctools.cscript):
  
             # Append effective area
             aeff = gammalib.GXmlElement('parameter name="EffectiveArea"')
-            aeff.attribute('file',aefffile)
+            aeff.attribute('file', aefffile)
  
             # Append PSF
             psf = gammalib.GXmlElement('parameter name="PointSpreadFunction"')
-            psf.attribute('file',psffile)
+            psf.attribute('file', psffile)
  
             # Append energy dispersion
             edisp = gammalib.GXmlElement('parameter name="EnergyDispersion"')
-            edisp.attribute('file',edispfile)
+            edisp.attribute('file', edispfile)
  
             # Append background
             bck = gammalib.GXmlElement('parameter name="Background"')
-            bck.attribute('file',bkgfile)
+            bck.attribute('file', bkgfile)
  
             # Assign events and IRFs to observation
             obs.append(ev)
@@ -726,8 +742,8 @@ class csiactobs(ctools.cscript):
  
         # Write warning in no observation could be used from runlist
         if lib.size() == 0:
-            self._log('WARNING: No observation was appended from specified '
-                      'runlist')
+            self._log_string(gammalib.TERSE, 'WARNING: No observation was appended from specified '
+                      'runlist\n')
  
         # Return
         return       
@@ -386,6 +386,19 @@ class csobsdef(ctools.cscript):
         # Return
         return
  
+    def obs(self):
+        """
+        Returns observation container
+        -----
+        Returns
+        obs : `~gammalib.GObservations`
+            Observation container
+        """
+        
+        # Return container
+        return self._obs
+
+    
     def pntdef(self, csv):
         """
         Set pointing definition from a CSV table
-0 
+15000
+15001
+15002
+15003
+15004
+28000
+28001
+28002
+28003
+15000
+15001
+15002
+15003
+15004
+28000
+28001
+28002
+28003
 {
-    "datasets": [
-        {
-            "name": "unit-test",
-	    "obsindx": "obs-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"
-        }
-    ]
-}
+  "datasets": [
+    {
+      "name": "unit-test", 
+      "hduindx": "hdu-index.fits.gz", 
+      "obsindx": "obs-index.fits.gz", 
+      "Total observations": 9, 
+      "Total time": "0.07 hours", 
+      "Total events": 533, 
+      "Disk usage": "0.43 MB", 
+      "Date created": "2016-08-03"
+    }
+  ]
+}
 \ No newline at end of file
@@ -74,7 +74,7 @@ class Test(test):
         # Setup csfindobs command
         cmd = csfindobs+' datapath="'+self._datapath+'"'+ \
                         ' prodname="unit-test"'+ \
-                        ' ra=83.63 dec=22.01 rad=1.0'+ \
+                        ' ra=83.63 dec=22.01 rad=1.5'+ \
                         ' outfile="runlist_cmd1.dat"'+ \
                         ' logfile="csfindobs_cmd1.log" chatter=1'
  
@@ -87,12 +87,12 @@ class Test(test):
              'Check successful execution from command line')
  
         # Check run list
-        self._check_runlist('runlist_cmd1.dat')
+        self._check_runlist('runlist_cmd1.dat',5)
  
         # Setup csfindobs command
         cmd = csfindobs+' datapath="path_that_does_not_exist"'+ \
                         ' prodname="unit-test"'+ \
-                        ' ra=83.63 dec=22.01 rad=1.0'+ \
+                        ' ra=83.63 dec=22.01 rad=1.5'+ \
                         ' outfile="runlist_cmd2.dat"'+ \
                         ' logfile="csfindobs_cmd2.log"'
  
@@ -114,10 +114,10 @@ class Test(test):
         findobs['prodname'] = 'unit-test'
         findobs['ra']       = 83.63
         findobs['dec']      = 22.01
-        findobs['rad']      = 1.0
+        findobs['rad']      = 1.5
         findobs['outfile']  = 'runlist_py1.dat'
         findobs['logfile']  = 'csfindobs_py1.log'
-        findobs['chatter']  = 2
+        findobs['chatter']  = 4
  
         # Run csfindobs script
         findobs.logFileOpen()   # Make sure we get a log file
@@ -125,16 +125,16 @@ class Test(test):
  
         # Check list of runs
         runs = findobs.runs()
-        self.test_value(len(runs), 1,
+        self.test_value(len(runs), 5,
                         'Check for number of runs in runlist')
-        self.test_value(str(runs[0]), '0',
+        self.test_value(str(runs[0]), '15000',
                         'Check run ID in runlist')
  
         # Save run list
         findobs.save()
  
         # Check run list
-        self._check_runlist('runlist_py1.dat')
+        self._check_runlist('runlist_py1.dat', runs=5)
  
         # Set-up csfindobs for undefined spatial parameters
         findobs = cscripts.csfindobs()
@@ -151,16 +151,16 @@ class Test(test):
         findobs.execute()
  
         # Check run list
-        self._check_runlist('runlist_py2.dat')
+        self._check_runlist('runlist_py2.dat', runs=9)
  
         # Set-up csfindobs for user expression
         findobs = cscripts.csfindobs()
         findobs['datapath']   = self._datapath
         findobs['prodname']   = 'unit-test'
-        findobs['ra']         = 83.63
-        findobs['dec']        = 22.01
-        findobs['rad']        = 0.0
-        findobs['expression'] = 'ONTIME<10.0'
+        findobs['ra']         = 329.71
+        findobs['dec']        = -30.2339
+        findobs['rad']        = 2.0
+        findobs['expression'] = 'ONTIME<5.0'
         findobs['outfile']    = 'runlist_py3.dat'
         findobs['logfile']    = 'csfindobs_py3.log'
         findobs['chatter']    = 4
@@ -195,7 +195,7 @@ class Test(test):
         if runs > 0:
             self.test_value(runlist.ncols(), 1,
                             'Check for single column in runlist file')
-            self.test_value(runlist[0,0], '0',
+            self.test_value(runlist[0,0], '15000',
                             'Check for run ID in runlist file')
  
         # Return
@@ -18,6 +18,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
 import gammalib
 import cscripts
 from testing import test
@@ -43,7 +44,9 @@ class Test(test):
         test.__init__(self)
  
         # Set data members
-        self._remote = self._datadir + '/iactdata/master.json'
+        self._remote    = self._datadir + '/iactdata/master.json'
+        self._runlist1  = self._datadir + '/iact_runlist_1.dat'
+        self._runlist2  = self._datadir + '/iact_runlist_2.dat'
  
         # Return
         return
@@ -107,20 +110,115 @@ class Test(test):
         Test csiactcopy from Python
         """
         # Set-up csiactcopy
-        findobs = cscripts.csiactcopy()
-        findobs['remote_master'] = self._remote
-        findobs['prodname']      = 'unit-test'
-        findobs['outpath']       = 'iactdata_py1'
-        findobs['logfile']       = 'csiactcopy_py1.log'
-        findobs['chatter']       = 2
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py1'
+        iactcopy['logfile']       = 'csiactcopy_py1.log'
+        iactcopy['chatter']       = 4
  
         # Run csiactcopy script and save run list
-        findobs.logFileOpen()   # Make sure we get a log file
-        findobs.run()
-        findobs.save()
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+        
+        # Check IACT data store
+        self._check_copy('iactdata_py1')
+        self._check_n_obs('iactdata_py1', 9)
+        
+        # Set-up csiactcopy
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py1'
+        iactcopy['logfile']       = 'csiactcopy_py2.log'
+        iactcopy['chatter']       = 4
+        iactcopy['clobber']       = True
  
+        # Run csiactcopy script and save run list
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+        
         # Check copy
         self._check_copy('iactdata_py1')
+        self._check_n_obs('iactdata_py1', 9)
+        
+        # Set-up csiactcopy
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py1'
+        iactcopy['logfile']       = 'csiactcopy_py3.log'
+        iactcopy['chatter']       = 4
+        iactcopy['clobber']       = False
+
+        # Run csiactcopy script and save run list
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+        
+        # Check copy
+        self._check_copy('iactdata_py1')
+        self._check_n_obs('iactdata_py1', 9)
+        
+        if os.path.isdir('iactdata_py2'):
+            os.system("rm -r iactdata_py2")
+        
+        # Set-up csiactcopy
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['runlist']       = self._runlist1
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py2'
+        iactcopy['logfile']       = 'csiactcopy_py4.log'
+        iactcopy['chatter']       = 3       
+
+        # Run csiactcopy script and save run list
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+
+        # Check copy
+        self._check_copy('iactdata_py2')
+        self._check_n_obs('iactdata_py2', 5)
+
+        # Set-up csiactcopy
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['runlist']       = self._runlist2
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py2'
+        iactcopy['logfile']       = 'csiactcopy_py5.log'
+        iactcopy['chatter']       = 3
+        
+        # Run csiactcopy script and save run list
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+
+        # Check copy
+        self._check_copy('iactdata_py2')
+        self._check_n_obs('iactdata_py2', 9)
+        
+        # Set-up csiactcopy
+        iactcopy = cscripts.csiactcopy()
+        iactcopy['runlist']       = self._runlist2
+        iactcopy['remote_master'] = self._remote
+        iactcopy['prodname']      = 'unit-test'
+        iactcopy['outpath']       = 'iactdata_py2'
+        iactcopy['logfile']       = 'csiactcopy_py6.log'
+        iactcopy['chatter']       = 3
+        iactcopy['clobber']       = True
+        
+        # Run csiactcopy script and save run list
+        iactcopy.logFileOpen()   # Make sure we get a log file
+        iactcopy.run()
+        iactcopy.save()
+
+        # Check copy
+        self._check_copy('iactdata_py2')
+        self._check_n_obs('iactdata_py2',9)
  
         # Return
         return
@@ -151,3 +249,27 @@ class Test(test):
  
         # Return
         return
+
+    # Check copy
+    def _check_n_obs(self, pathname, n_expected):
+        """
+        Check number of available observations.
+        """
+        # Set file name
+        obs_index_name = gammalib.GFilename(pathname+'/obs-index.fits[OBS_INDEX]')
+        
+        # Open index file
+        fits = gammalib.GFits(obs_index_name)
+        
+        # Get number of observations
+        n_obs = fits[obs_index_name.extname()].nrows()
+        
+        # Close FITS file
+        fits.close()
+
+        # Check for existence of observations
+        self.test_assert(n_obs == n_expected,
+                         'Check for '+str(n_expected)+' observations.')
+                
+        # Return
+        return
 \ No newline at end of file
@@ -90,10 +90,10 @@ class Test(test):
              'Check successful execution from command line')
  
         # Check observation definition XML file
-        self._check_obsdef('csiactobs_obs_cmd1.xml')
+        self._check_obsdef('csiactobs_obs_cmd1.xml', 6)
  
         # Check model definition XML file
-        self._check_moddef('csiactobs_bgd_cmd1.xml')
+        self._check_moddef('csiactobs_bgd_cmd1.xml', 6)
  
         # Setup csiactobs command
         cmd = csiactobs+' datapath="data_path_that_does_not_exist"'+ \
@@ -146,13 +146,13 @@ class Test(test):
         iactobs.save()
  
         # Check observation definition XML file
-        self._check_obsdef('csiactobs_obs_py1.xml')
+        self._check_obsdef('csiactobs_obs_py1.xml', 6)
  
         # Check model definition XML file
-        self._check_moddef('csiactobs_bgd_py1.xml')
+        self._check_moddef('csiactobs_bgd_py1.xml', 6)
  
         # Create test runlist
-        runlist = ['0']
+        runlist = ['15000','15001']
  
         # Set-up csiactobs using a runlist with 2 background parameters
         iactobs = cscripts.csiactobs()
@@ -171,16 +171,16 @@ class Test(test):
         iactobs.save()
  
         # Test return functions
-        self.test_value(iactobs.obs().size(), 1,
+        self.test_value(iactobs.obs().size(), 2,
                         'Check number of observations in container')
         self.test_value(iactobs.ebounds().size(), 0,
                         'Check number of energy boundaries')
  
         # Check observation definition XML file
-        self._check_obsdef('csiactobs_obs_py2.xml')
+        self._check_obsdef('csiactobs_obs_py2.xml',2)
  
         # Check model definition XML file
-        self._check_moddef('csiactobs_bgd_py2.xml')
+        self._check_moddef('csiactobs_bgd_py2.xml',2)
  
         # Set-up csiactobs with a large number of free parameters and "aeff"
         # background
@@ -199,10 +199,10 @@ class Test(test):
         iactobs.execute()
  
         # Check observation definition XML file
-        self._check_obsdef('csiactobs_obs_py3.xml')
+        self._check_obsdef('csiactobs_obs_py3.xml',6)
  
         # Check model definition XML file
-        self._check_moddef('csiactobs_bgd_py3.xml')
+        self._check_moddef('csiactobs_bgd_py3.xml',6)
  
         # Set-up csiactobs with a "gauss" background and "inmodel" parameter
         iactobs = cscripts.csiactobs()
@@ -222,56 +222,85 @@ class Test(test):
         iactobs.run()
  
         # Check number of observations
-        self.test_value(iactobs.obs().size(), 1,
+        self.test_value(iactobs.obs().size(), 6,
                         'Check number of observations in container')
  
         # Check number of models
-        self.test_value(iactobs.obs().models().size(), 3,
-                        'Check number of models in container')
+        self.test_value(iactobs.obs().models().size(), 8,
+                        'Check number of models in container')        
+        
+        # Set-up csiactobs with a "gauss" background and "inmodel" parameter
+        iactobs = cscripts.csiactobs()
+        iactobs['datapath']      = self._datapath
+        iactobs['inmodel']       = self._model
+        iactobs['prodname']      = 'unit-test'
+        iactobs['infile']        = self._runlist
+        iactobs['bkgpars']       = 1
+        iactobs['bkg_mod_hiera'] = 'irf'
+        iactobs['outobs']        = 'NONE'
+        iactobs['outmodel']      = 'NONE'
+        iactobs['logfile']       = 'csiactobs_py4.log'
+        iactobs['chatter']       = 4
+   
+        # Run csiactobs script
+        iactobs.logFileOpen()   # Make sure we get a log file
+        iactobs.run()
+        
+        # Check number of observations
+        self.test_value(iactobs.obs().size(), 5,
+                        'Check number of observations in container')
+
+        # Check number of models
+        self.test_value(iactobs.obs().models().size(), 7,
+                        'Check number of models in container')     
  
         # Return
         return
  
     # Check observation definition XML file
-    def _check_obsdef(self, filename):
+    def _check_obsdef(self, filename, obs_expected):
         """
         Check observation definition XML file
         """
         # Load observation definition XML file
         obs = gammalib.GObservations(filename)
  
-        # Get response
-        rsp = obs[0].response()
-
         # Check number of observations
-        self.test_value(obs.size(), 1,
-                        'Check for single observation in XML file')
-        self.test_value(obs[0].eventfile().file(), 'events.fits.gz',
-                        'Check event file name')
-        self.test_value(obs[0].eventfile().extname(), 'EVENTS',
-                        'Check event extension name')
-        self.test_value(rsp.aeff().filename().file(), 'irf_file.fits.gz',
-                        'Check effective area file name')
-        self.test_value(rsp.aeff().filename().extname(), 'EFFECTIVE AREA',
-                        'Check effective area extension name')
-        self.test_value(rsp.psf().filename().file(), 'irf_file.fits.gz',
-                        'Check point spread function file name')
-        self.test_value(rsp.psf().filename().extname(), 'POINT SPREAD FUNCTION',
-                        'Check point spread function extension name')
-        self.test_value(rsp.edisp().filename().file(), 'irf_file.fits.gz',
-                        'Check energy dispersion file name')
-        self.test_value(rsp.edisp().filename().extname(), 'ENERGY DISPERSION',
-                        'Check energy dispersion extension name')
-        self.test_value(rsp.background().filename().file(), 'irf_file.fits.gz',
-                        'Check background file name')
-        self.test_value(rsp.background().filename().extname(), 'BACKGROUND',
-                        'Check background extension name')
+        self.test_value(obs.size(), obs_expected,
+                        'Check for '+str(obs_expected)+' observations in XML file')
+        
+        if obs_expected > 0:
+            
+            # Get response
+            rsp = obs[0].response()
+            
+            # Test response
+            self.test_value(obs[0].eventfile().file(), 'events_0.fits.gz',
+                            'Check event file name')
+            self.test_value(obs[0].eventfile().extname(), 'EVENTS',
+                            'Check event extension name')
+            self.test_value(rsp.aeff().filename().file(), 'irf_file.fits.gz',
+                            'Check effective area file name')
+            self.test_value(rsp.aeff().filename().extname(), 'EFFECTIVE AREA',
+                            'Check effective area extension name')
+            self.test_value(rsp.psf().filename().file(), 'irf_file.fits.gz',
+                            'Check point spread function file name')
+            self.test_value(rsp.psf().filename().extname(), 'POINT SPREAD FUNCTION',
+                            'Check point spread function extension name')
+            self.test_value(rsp.edisp().filename().file(), 'irf_file.fits.gz',
+                            'Check energy dispersion file name')
+            self.test_value(rsp.edisp().filename().extname(), 'ENERGY DISPERSION',
+                            'Check energy dispersion extension name')
+            self.test_value(rsp.background().filename().file(), 'irf_file.fits.gz',
+                            'Check background file name')
+            self.test_value(rsp.background().filename().extname(), 'BACKGROUND',
+                            'Check background extension name')
  
         # Return
         return
  
     # Check model XML file
-    def _check_moddef(self, filename):
+    def _check_moddef(self, filename, models_expected):
         """
         Check model definition XML file
         """
@@ -279,8 +308,8 @@ class Test(test):
         models = gammalib.GModels(filename)
  
         # Check number of models
-        self.test_value(models.size(), 1,
-                        'Check for single model in XML file')
+        self.test_value(models.size(), models_expected,
+                        'Check for '+str(models_expected)+' models in XML file')
  
         # Return
         return