Merge branch '1828-enhance-csiactdata-tests' into integration

Jürgen Knödlseder
2 parents 73847eac 4622f176
Showing 24 changed files with 496 additions and 274 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,30 +155,44 @@ 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_2.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/obs-index.fits.gz
+
+# Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata/events
+testiactevents = $(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
 # it is)
-binscriptsdir      = $(bindir)
-binscripts_DATA    = $(binscripts)
-pythonscriptsdir   = $(pyexecdir)/cscripts
-pythonscripts_DATA = $(pythonscripts)
-pfilesdir          = $(prefix)/syspfiles
-pfiles_DATA        = $(parfiles)
-testscriptsdir     = $(pyexecdir)/cscripts/tests
-testscripts_DATA   = $(testscripts)
-testdatadir        = $(pyexecdir)/cscripts/tests/data
-testdata_DATA      = $(testdata)
-testiactdatadir    = $(pyexecdir)/cscripts/tests/data/iactdata
-testiactdata_DATA  = $(testiactdata)
+binscriptsdir       = $(bindir)
+binscripts_DATA     = $(binscripts)
+pythonscriptsdir    = $(pyexecdir)/cscripts
+pythonscripts_DATA  = $(pythonscripts)
+pfilesdir           = $(prefix)/syspfiles
+pfiles_DATA         = $(parfiles)
+testscriptsdir      = $(pyexecdir)/cscripts/tests
+testscripts_DATA    = $(testscripts)
+testdatadir         = $(pyexecdir)/cscripts/tests/data
+testdata_DATA       = $(testdata)
+testiactdatadir     = $(pyexecdir)/cscripts/tests/data/iactdata
+testiactdata_DATA   = $(testiactdata)
+testiacteventsdir   = $(pyexecdir)/cscripts/tests/data/iactdata/events
+testiactevents_DATA = $(testiactdata)
 # Files to be distributed in the tar file
 EXTRA_DIST = $(binscripts) $(pythonscripts) $(parfiles) \
@@ -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:
@@ -172,8 +172,10 @@ class csiactobs(ctools.cscript):
             # 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']))
+                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
@@ -292,7 +294,9 @@ class csiactobs(ctools.cscript):
                 spec = gammalib.GModelSpectralNodes()
                 # Create logarithmic energy nodes
-                bounds = gammalib.GEbounds(self._bkgpars,gammalib.GEnergy(emin,'TeV'),gammalib.GEnergy(emax,'TeV'), True)
+                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()):     
@@ -319,7 +323,8 @@ class csiactobs(ctools.cscript):
         # Return spectrum
         return spec
-    def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100):
+    def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype,
+                         emin=0.01, emax=100):
         """
         Create an IACT background model
@@ -348,27 +353,35 @@ class csiactobs(ctools.cscript):
             prefactor *= bkg_scale
             # Create background spectrum
-            spec       = self._background_spectrum(prefactor, index, emin, emax)
+            spec = self._background_spectrum(prefactor, index, emin, emax)
-            # Create background model instance
+            # Create background model
             bck = gammalib.GCTAModelIrfBackground(spec)
         # Set AeffBackground   
         elif bkgtype == 'aeff':
+
+            # Set background components
             prefactor = bkg_scale * self._bkg_aeff_norm
-            spec      = self._background_spectrum(prefactor, self._bkg_aeff_index, emin, emax)
+            spec      = self._background_spectrum(prefactor,
+                                                  self._bkg_aeff_index,
+                                                  emin, emax)
-            # Create background model instance
+            # Create background model
             bck = gammalib.GCTAModelAeffBackground(spec)
         # Set Gaussian Background
         elif bkgtype == 'gauss':
+
+            # Set background components
             prefactor = bkg_scale * self._bkg_gauss_norm
-            
-            # Create background spectrum
-            spec      = self._background_spectrum(prefactor, self._bkg_gauss_index, emin, emax)
+            spec      = self._background_spectrum(prefactor,
+                                                  self._bkg_gauss_index,
+                                                  emin, emax)
             radial    = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma)
-            bck       = gammalib.GCTAModelRadialAcceptance(radial, spec)
+
+            # Create background model
+            bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
         else:
             msg = 'Background type "'+bkgtype+'" unsupported'
@@ -422,7 +435,7 @@ class csiactobs(ctools.cscript):
         Write observation summary
         """
         # Log header
-        self._log_header1(gammalib.TERSE, 'Observation summary')
+        self._log_header1(gammalib.NORMAL, 'Observation summary')
         # Set energy range dependent on whether boundaries exist or
         # not
@@ -438,6 +451,95 @@ 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 str
+            List of strings containing the hierarchy how to look for files
+        formats : list of str
+            File formats available for this observation
+        
+        Returns
+        -------
+        filename : str
+            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 : str
+            Filename
+        obs_id : int
+            Observation ID
+        filetype : str
+            Type of file
+        
+        Returns
+        -------
+        present : bool
+            True if filename is present, False otherwise
+        """
+                
+        # Initialise return value
+        present = True
+        
+        # Create filename instance
+        fname = gammalib.GFilename(filename)
+        
+        # Check if file is empty
+        if fname.is_empty():
+            msg = 'Skipping observation "%s": No %s found' % (obs_id, filetype)
+            self._log_string(gammalib.NORMAL, msg)
+            present = False
+        
+        # Check if file exists
+        elif not fname.exists():
+            msg = ('Skipping observation "%s": %s "%s" does not exist' %
+                   (obs_id, filetype, filename))
+            self._log_string(gammalib.NORMAL, msg)
+            present = False
+        
+        # Return present flag
+        return present    
+
     # Public methods
     def run(self):
@@ -460,9 +562,8 @@ class csiactobs(ctools.cscript):
         lib = self._xml.element('observation_list', 0)
         # Log output
-        if self._logTerse():
-            self._log('\n')
-            self._log.header1('Looping over '+str(len(self._runlist))+' runs')
+        msg = 'Looping over %d runs' % len(self._runlist)
+        self._log_header1(gammalib.TERSE, msg)
         # Initialise energy range values for logging
         self._ebounds.clear()
@@ -474,152 +575,86 @@ class csiactobs(ctools.cscript):
             obs_selection = '[OBS_ID=='+str(obs_id)+']'
             # Open HDU index file
-            hduindx     = gammalib.GFits(self._hdu_index+'[HDU_INDEX]'+obs_selection)
+            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 = ''
+                msg = ('Warning: observation "%s" has no energy dispersion '
+                       'information' % obs_id)
+                self._log_string(gammalib.NORMAL, msg)
-            # 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
+                    msg = ('Warning: observation "%s" has no background '
+                           'information (IRF background cannot be used)' %
+                           obs_id)
+                    self._log_string(gammalib.NORMAL, msg)
+                
+                # Check if background model hierarchy is empty
+                if len(run_bkg_mod_hierarchy) == 0:
+                    
+                    # Skip observation if no background can be used
+                    msg = ('Skipping observation "%s": No background can be '
+                           'used' % obs_id)
+                    self._log_string(gammalib.NORMAL, msg)
+                    continue
+                
+                else:
+                    # Log if we fall back to next background approach
+                    msg = ('Observation "%s": Falling back to background "%s"' %
+                           (obs_id, run_bkg_mod_hierarchy[0]))
+                    self._log_string(gammalib.NORMAL, msg)
             # Close hdu index file
             hduindx.close()
-
+            
             # Initialise background scale
             bkg_scale = 1.0
@@ -627,7 +662,8 @@ class csiactobs(ctools.cscript):
             if self._use_bkg_scale:
                 # Read background scale from fits file
-                obsindx   = gammalib.GFits(self._obs_index+'[OBS_INDEX]'+obs_selection)
+                obsindx   = gammalib.GFits(self._obs_index+'[OBS_INDEX]'+
+                                           obs_selection)
                 bkg_scale = obsindx['OBS_INDEX']['BKG_SCALE'][0]
                 obsindx.close()
@@ -635,7 +671,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 +680,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 +704,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 +720,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)
@@ -706,17 +742,20 @@ class csiactobs(ctools.cscript):
             obs.append(bck)
             # Log the observation ID and object name that has been appended
-            self._log_value(gammalib.TERSE, 'Append observation', '%s ("%s")' %
-                                            (obs_id, object_name))
+            self._log_value(gammalib.NORMAL, 'Append observation', '%s ("%s")' %
+                            (obs_id, object_name))
             # Log the file names of the event and response files
             self._log_value(gammalib.EXPLICIT, ' Event file', eventfile)
             self._log_value(gammalib.EXPLICIT, ' Effective area file', aefffile)
-            self._log_value(gammalib.EXPLICIT, ' Point spread function file', psffile)
-            self._log_value(gammalib.EXPLICIT, ' Energy dispersion file', edispfile)
+            self._log_value(gammalib.EXPLICIT, ' Point spread function file',
+                            psffile)
+            self._log_value(gammalib.EXPLICIT, ' Energy dispersion file',
+                            edispfile)
             self._log_value(gammalib.EXPLICIT, ' Background file', bkgfile)
             if self._use_bkg_scale:
-                self._log_value(gammalib.EXPLICIT, ' Background scale', bkg_scale)
+                self._log_value(gammalib.EXPLICIT, ' Background scale',
+                                bkg_scale)
         # Append models provided by "inmodel" to the model container
         self._append_inmodels()
@@ -726,8 +765,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.NORMAL, 'WARNING: No observation was '
+                             'appended from specified runlist')
         # Return
         return       
@@ -756,28 +795,22 @@ class csiactobs(ctools.cscript):
         Save observation definition and model definition XML file
         """
         # Write header
-        if self._logTerse():
-            self._log('\n')
-            self._log.header1('Save XML files')
+        self._log_header1(gammalib.TERSE, 'Save XML files')
         # Get filenames
         outobs   = self['outobs'].filename()
         outmodel = self['outmodel'].filename()
-        # Log file name
-        if self._logTerse():
-            self._log(gammalib.parformat('Obs. definition XML file'))
-            self._log(outobs.url())
-            self._log('\n')
+        # Log observation definition file name
+        self._log_value(gammalib.NORMAL, 'Obs. definition XML file',
+                        outobs.url())
         # Save observation definition XML file
         self._xml.save(outobs.url())
-        # Log file name
-        if self._logTerse():
-            self._log(gammalib.parformat('Model definiton XML file'))
-            self._log(outmodel.url())
-            self._log('\n')
+        # Log model definition file name
+        self._log_value(gammalib.NORMAL, 'Model definiton XML file',
+                        outmodel.url())
         # Save model XML file
         self._models.save(outmodel.url())
@@ -336,7 +336,7 @@ class csobsdef(ctools.cscript):
             # Write observation into logger
             name  = obs.instrument()+' observation'
             value = 'Name="%s" ID="%s"' % (obs.name(), obs.id())
-            self._log_value(gammalib.TERSE, name, value)
+            self._log_value(gammalib.NORMAL, name, value)
             self._log_string(gammalib.EXPLICIT, str(obs)+'\n')
             # Append observation
@@ -386,6 +386,17 @@ 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