Add Sphinx tutorials test script (#3510)

The current script covers all 1DC tutorials.

Add Sphinx tutorials test script (#3510)
The current script covers all 1DC tutorials.
Jürgen Knödlseder
1 parent e1129812
Showing 12 changed files with 398 additions and 149 deletions
.gitignore
ChangeLog
Makefile.am
NEWS
README.md
cscripts/csspec.py
doc/source/users/tutorials/1dc/CasA_model.xml
doc/source/users/tutorials/1dc/first_fitting.rst
doc/source/users/tutorials/1dc/first_onoff.rst
doc/source/users/tutorials/1dc/results_stacked_eplaw.xml
doc/source/users/tutorials/1dc/show_onoff_spectrum.py
test/tutorials.py
@@ -99,6 +99,7 @@ test/butterfly.txt
 test/.coverage
 test/test_python_ctools.sh
 test/test_python_cscripts.sh
+test/tutorials
 examples/*.fits
 examples/*.xml
 examples/*.dat
-2020-12-03
+2020-02-03
         * Version 2.0.0 released
         ========================
+        Add Sphinx tutorials test script (#3510)
         Add "refit_if_failed" parameter to ctlike
         Add models() method to ctobssim (#3441)
         Add ctlike results to XML file (#3431)
 #############################################################################
 # Makefile for ctools package                                               #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2016-2017 Juergen Knoedlseder                               #
+# Copyright (C) 2016-2021 Juergen Knoedlseder                               #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
@@ -211,13 +211,23 @@ valgrind:
 # Add science verification rule
 science-verification:
 	@if type test/science_verification.py >/dev/null 2>&1; then \
-	    cd test; \
-	    ./science_verification.py; \
-	    cd ..; \
+		cd test; \
+		./science_verification.py; \
+		cd ..; \
 	else \
 		echo "science_verification.py not found in test."; \
 	fi
+# Add tutorials verification rule
+check-tutorials:
+	@if type test/tutorials.py >/dev/null 2>&1; then \
+		cd test; \
+		./tutorials.py; \
+		cd ..; \
+	else \
+		echo "tutorials.py not found in test."; \
+	fi
+
 # Add Sphinx rule
 sphinx:
 	@if type sphinx-build >/dev/null 2>&1; then \
 New Features and Important Changes in ctools 2.0.0
-3 December 2020
+3 February 2021
 Introduction
@@ -16,7 +16,8 @@ None
 Documentation
 -------------
-None
+Add test/tutorials.py script and check-tutorial target to test the code in the
+Sphinx tutorials (#3510).
 Python interface
 ctools information
 ==================
-* Version:             2.0.0.dev (3 December 2020)
+* Version:             2.0.0.dev (3 February 2021)
 * GammaLib dependency: 2.0.0.dev
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a spectrum.
 #
-# Copyright (C) 2014-2020 Michael Mayer
+# Copyright (C) 2014-2021 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -1071,7 +1071,7 @@ class csspec(ctools.csobservation):
         energy_high  = gammalib.GFitsTableDoubleCol('e_max', nrows)
         norm         = gammalib.GFitsTableDoubleCol('norm', nrows)
         norm_err     = gammalib.GFitsTableDoubleCol('norm_err', nrows)
-        ulim_values  = gammalib.GFitsTableDoubleCol('norm_ul', nrows)
+        norm_ul      = gammalib.GFitsTableDoubleCol('norm_ul', nrows)
         e2dnde       = gammalib.GFitsTableDoubleCol('ref_e2dnde', nrows)
         dnde         = gammalib.GFitsTableDoubleCol('ref_dnde', nrows)
         Npred_values = gammalib.GFitsTableDoubleCol('ref_npred', nrows)
@@ -1084,7 +1084,7 @@ class csspec(ctools.csobservation):
         energy_high.unit('TeV')
         norm.unit('')
         norm_err.unit('')
-        ulim_values.unit('')
+        norm_ul.unit('')
         e2dnde.unit('erg/cm2/s')
         dnde.unit('counts/MeV/cm2/s')
         Npred_values.unit('counts')
@@ -1098,8 +1098,12 @@ class csspec(ctools.csobservation):
             energy_low[i]   = result['energy_low']
             energy_high[i]  = result['energy_high']
             norm[i]         = 1.0
-            norm_err[i]     = result['flux_err'] / result['flux']
-            ulim_values[i]  = result['ulimit'] / result['flux']
+            if result['flux'] != 0.0:
+                norm_err[i] = result['flux_err'] / result['flux']
+                norm_ul[i]  = result['ulimit']   / result['flux']
+            else:
+                norm_err[i] = 0.0
+                norm_ul[i]  = 0.0
             e2dnde[i]       = result['e2dnde']
             dnde[i]         = result['flux']
             Npred_values[i] = result['Npred']
@@ -1109,7 +1113,10 @@ class csspec(ctools.csobservation):
             # Add likelihood scan values
             for fbin in range(ncols):
                 dloglike_scan[i,fbin] = result['dloglike'][fbin]
-                norm_scan[i,fbin]     = result['norm_scan'][fbin] / result['flux']
+                if result['flux'] != 0.0:
+                    norm_scan[i,fbin] = result['norm_scan'][fbin] / result['flux']
+                else:
+                    norm_scan[i,fbin] = 0.0
         # Initialise FITS Table with extension "SPECTRUM"
         table = gammalib.GFitsBinTable(nrows)
@@ -1125,7 +1132,7 @@ class csspec(ctools.csobservation):
         table.append(energy_high)
         table.append(norm)
         table.append(norm_err)
-        table.append(ulim_values)
+        table.append(norm_ul)
         table.append(e2dnde)
         table.append(dnde)
         table.append(TSvalues)
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<source_library title="source library">
+  <source name="Cassiopeia A" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   value="1.45" scale="1e-18" min="0"            free="1"/>
+      <parameter name="Index"       value="2.75" scale="-1"    min="-10" max="10" free="1"/>
+      <parameter name="PivotEnergy" value="1"    scale="1e6"                      free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA"  value="350.8500" scale="1" free="0"/>
+      <parameter name="DEC" value="58.8150"  scale="1" free="0"/>
+    </spatialModel>
+  </source>
+  <source name="Background model" type="CTAIrfBackground" instrument="CTA">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" scale="1"   min="0.001" max="1000" free="1"/>
+      <parameter name="Index"     value="0" scale="1"   min="-5"    max="5"    free="1"/>
+      <parameter name="Scale"     value="1" scale="1e6" min="0.01"  max="1000" free="0"/>
+    </spectrum>
+  </source>
+</source_library>
@@ -39,7 +39,7 @@ You run the :ref:`ctlike` tool as follows:
    Input model definition XML file [$CTOOLS/share/models/crab.xml] bkgcube.xml
    Output model definition XML file [crab_results.xml] results_stacked.xml
-The tool will take a few minutes (on Mac OS X) to perform the model fitting,
+The tool will take a few minutes (on Mac OS) to perform the model fitting,
 and will write the results into an updated
 :ref:`model definition file <glossary_moddef>`
 containing the fitted model parameters and their statistical uncertainties.
@@ -124,7 +124,7 @@ You create the butterfly diagram for ``Src001`` using
    Stop value for last energy bin in TeV [100.0]
    Output ASCII file [butterfly.txt] butterfly_src001.txt
-and for `Src003`` using
+and for ``Src003`` using
 .. code-block:: bash
@@ -192,28 +192,14 @@ a maximum likelihood model fit.
 Obviously, ``Src001`` has a spectral cut-off (red flux points) and hence is not
 adequately described by a power law model. You should therefore replace the
-power law in the
-:ref:`model definition file <glossary_moddef>`
-by an exponentially cutoff power law, as shown below:
-
-.. code-block:: xml
-
-   <?xml version="1.0" encoding="UTF-8" standalone="no"?>
-   <source_library title="source library">
-     <source name="Src001" type="PointSource">
-       <spectrum type="ExponentialCutoffPowerLaw">
-         <parameter name="Prefactor"    scale="1e-18" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-         <parameter name="Index"        scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-         <parameter name="CutoffEnergy" scale="1e7"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
-         <parameter name="PivotEnergy"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
-       </spectrum>
-       <spatialModel type="PointSource">
-         <parameter name="RA"  value="266.424004498437"  scale="1" free="1" />
-         <parameter name="DEC" value="-29.0049010253548" scale="1" free="1" />
-       </spatialModel>
-     </source>
-     ...
-   </source_library>
+power law spectrum in the :ref:`model definition file <glossary_moddef>`
+by an exponentially cutoff power law spectrum, as shown below (only the first
+lines of the :ref:`model definition file <glossary_moddef>` are shown for
+illustration):
+
+.. literalinclude:: results_stacked_eplaw.xml
+   :language: xml
+   :lines: 1-16
 Fitting this model to the data improves the fit and the resulting butterfly
 diagram follows now reasonably well the spectral points:
@@ -208,6 +208,7 @@ follows:
    Output observation definition XML file [onoff_obs.xml]
    Output model definition XML file [onoff_model.xml]
    Method for background estimation (REFLECTED|CUSTOM) [REFLECTED]
+   Shape of the source region (CIRCLE|RECT) [CIRCLE]
    Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
    Right Ascension of source region centre (deg) (0-360) [83.63] 350.85
    Declination of source region centre (deg) (-90-90) [22.01] 58.815
@@ -355,29 +356,8 @@ The fit result can be inspected by peeking the log file:
    provided as input model to :ref:`csphagen`. Here an example for an
    input model:
-   .. code-block:: xml
-
-      <?xml version="1.0" encoding="UTF-8" standalone="no"?>
-      <source_library title="source library">
-        <source name="Cassiopeia A" type="PointSource">
-          <spectrum type="PowerLaw">
-            <parameter name="Prefactor"   value="1.45" scale="1e-18" min="0"            free="1"/>
-            <parameter name="Index"       value="2.75" scale="-1"    min="-10" max="10" free="1"/>
-            <parameter name="PivotEnergy" value="1"    scale="1e6"                      free="0"/>
-          </spectrum>
-          <spatialModel type="PointSource">
-            <parameter name="RA"  value="350.8500" scale="1" free="0"/>
-            <parameter name="DEC" value="58.8150"  scale="1" free="0"/>
-          </spatialModel>
-        </source>
-        <source name="Background model" type="CTAIrfBackground" instrument="CTA">
-          <spectrum type="PowerLaw">
-            <parameter name="Prefactor" value="1" scale="1"   min="0.001" max="1000" free="1"/>
-            <parameter name="Index"     value="0" scale="1"   min="-5"    max="5"    free="1"/>
-            <parameter name="Scale"     value="1" scale="1e6" min="0.01"  max="1000" free="0"/>
-          </spectrum>
-        </source>
-      </source_library>
+   .. literalinclude:: CasA_model.xml
+      :language: xml
    Now rerun :ref:`csphagen` as follows:
@@ -395,6 +375,7 @@ The fit result can be inspected by peeking the log file:
       Output observation definition XML file [onoff_obs.xml] onoff_obs_cstat.xml
       Output model definition XML file [onoff_model.xml] onoff_model_cstat.xml
       Method for background estimation (REFLECTED|CUSTOM) [REFLECTED]
+      Shape of the source region (CIRCLE|RECT) [CIRCLE]
       Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
       Right Ascension of source region centre (deg) (0-360) [350.85]
       Declination of source region centre (deg) (-90-90) [58.815]
@@ -452,7 +433,7 @@ Like for a binned or unbinned analysis, you can create a butterfly diagram by ty
    $ ctbutterfly
    Input event list, counts cube or observation definition XML file [events.fits] onoff_obs.xml
-   Source of interest [Crab] Cassiopeia A
+   Source of interest [Crab] Dummy
    Input model definition XML file [$CTOOLS/share/models/crab.xml] CasA_results.xml
    Lower energy limit (TeV) [0.1]
    Upper energy limit (TeV) [100.0] 50.0
@@ -465,7 +446,7 @@ and the SED by typing
    $ csspec
    Input event list, counts cube, or observation definition XML file [events.fits] onoff_obs.xml
    Input model definition XML file [$CTOOLS/share/models/crab.xml] CasA_results.xml
-   Source name [Crab] Cassiopeia A
+   Source name [Crab] Dummy
    Spectrum generation method (SLICE|NODES|AUTO) [AUTO]
    Binning algorithm (FILE|LIN|LOG|POW) [LOG]
    Lower energy limit (TeV) [0.1]
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<source_library title="source library">
+  <source name="Src001" type="PointSource">
+    <spectrum type="ExponentialCutoffPowerLaw">
+      <parameter name="Prefactor"    scale="1e-18" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"        scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="CutoffEnergy" scale="1e7"   value="1.0"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="PivotEnergy"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA"  value="266.424004498437"  scale="1" free="1"/>
+      <parameter name="DEC" value="-29.0049010253548" scale="1" free="1"/>
+    </spatialModel>
+  </source>
+  <source name="Src002" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   value="1" scale="5.7e-18" min="0" free="1" />
+      <parameter name="Index"       value="1" scale="-2.48"   min="-4.03225806451613" max="4.03225806451613" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA" value="264.802012181975" error="0" scale="1" free="1" />
+      <parameter name="DEC" value="-31.0050720432975" error="0" scale="1" free="1" />
+    </spatialModel>
+  </source>
+  <source name="Src003" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" error="0" scale="5.7e-18" min="0" free="1" />
+      <parameter name="Index" value="1" error="-0" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA" value="266.831945177213" error="0" scale="1" free="1" />
+      <parameter name="DEC" value="-28.1460284439951" error="0" scale="1" free="1" />
+    </spatialModel>
+  </source>
+  <source name="Src004" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" error="0" scale="5.7e-18" min="0" free="1" />
+      <parameter name="Index" value="1" error="-0" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA" value="268.244867737546" error="0" scale="1" free="1" />
+      <parameter name="DEC" value="-30.3728286753247" error="0" scale="1" free="1" />
+    </spatialModel>
+  </source>
+  <source name="Src005" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" error="0" scale="5.7e-18" min="0" free="1" />
+      <parameter name="Index" value="1" error="-0" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA" value="270.424492131939" error="0" scale="1" free="1" />
+      <parameter name="DEC" value="-29.6624142779832" error="0" scale="1" free="1" />
+    </spatialModel>
+  </source>
+  <source name="Src006" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" error="0" scale="5.7e-18" min="0" free="1" />
+      <parameter name="Index" value="1" error="-0" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA" value="263.809469183313" error="0" scale="1" free="1" />
+      <parameter name="DEC" value="-31.2230130593475" error="0" scale="1" free="1" />
+    </spatialModel>
+  </source>
+  <source name="BackgroundModel" type="CTACubeBackground" instrument="CTA,HESS,MAGIC,VERITAS">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor" value="1" error="0" scale="1" min="0.01" max="100" free="1" />
+      <parameter name="Index" value="0" error="0" scale="1" min="-5" max="5" free="1" />
+      <parameter name="PivotEnergy" value="1" scale="1000000" free="0" />
+    </spectrum>
+  </source>
+</source_library>
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Display butterfly generated by ctbutterfly and spectrum generated by csspec
 #
-# Copyright (C) 2018 Juergen Knoedlseder
+# Copyright (C) 2018-2021 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -51,13 +51,14 @@ def plot_spectrum_butterfly(specfile, butterfile, fmt=&#39;o&#39;, color=&#39;red&#39;):
     # Read spectrum file    
     fits     = gammalib.GFits(specfile)
     table    = fits.table(1)
-    c_energy = table['Energy']
-    c_ed     = table['ed_Energy']
-    c_eu     = table['eu_Energy']
-    c_flux   = table['Flux']
-    c_eflux  = table['e_Flux']
-    c_ts     = table['TS']
-    c_upper  = table['UpperLimit']
+    c_energy = table['e_ref']
+    c_ed     = table['e_min']
+    c_eu     = table['e_max']
+    c_flux   = table['ref_e2dnde']
+    c_norm   = table['norm']
+    c_eflux  = table['norm_err']
+    c_upper  = table['norm_ul']
+    c_ts     = table['ts']
     # Read butterfly file
     csv = gammalib.GCsv(butterfile)
@@ -83,22 +84,23 @@ def plot_spectrum_butterfly(specfile, butterfile, fmt=&#39;o&#39;, color=&#39;red&#39;):
         # Get Test Statistic, flux and flux error
         ts    = c_ts.real(row)
-        flx   = c_flux.real(row)
-        e_flx = c_eflux.real(row)
+        norm  = c_norm.real(row)
+        flx   = norm * c_flux.real(row)
+        e_flx = flx * c_eflux.real(row)
         # If Test Statistic is larger than 9 and flux error is smaller than
         # flux then append flux plots ...
         if ts > 9.0 and e_flx < flx:
             energies.append(c_energy.real(row))
-            flux.append(c_flux.real(row))
+            flux.append(flx)
             ed_engs.append(c_ed.real(row))
             eu_engs.append(c_eu.real(row))
-            e_flux.append(c_eflux.real(row))
+            e_flux.append(e_flx)
         # ... otherwise append upper limit
         else:
             ul_energies.append(c_energy.real(row))
-            ul_flux.append(c_upper.real(row))
+            ul_flux.append(flx*c_upper.real(row))
             ul_ed_engs.append(c_ed.real(row))
             ul_eu_engs.append(c_eu.real(row))
@@ -46,47 +46,128 @@ class tutorials(gammalib.GPythonTestSuite):
         # Call base class constructor
         gammalib.GPythonTestSuite.__init__(self)
-        # Initialise results
-        self.results = None
+        # Initialise private members
+        self._rst_step = 0     # Step counter
+        self._rst_path = ''    # Path where the tested RST file resides
         # Return
         return
-    # Set test functions
-    def set(self):
+    # Clean pfiles directory
+    def _clean_pfiles(self):
         """
-        Set all test functions
+        Clean pfiles directory
         """
-        # Set test name
-        self.name('Tutorials Verification')
-
-        # Append tutorials
-        self.append(self.tutorials_1dc, 'Test 1DC tutorial')
+        # Remove all files from pfiles directory
+        os.system('rm -rf %s/*.par' % os.environ['PFILES'])
         # Return
         return
-    # Test line execution
-    def _execute_cmd(self, cmd, args):
+    # Set result directory
+    def _set_result_dir(self, path):
         """
-        Execute command
+        Set result directory
+
+        Parameters
+        ----------
+        path : str
+            Result directory path
         """
-        # Create subprocess for tool
-        p = subprocess.Popen(cmd.split(' '), stdout=subprocess.PIPE,
-                                             stdin=subprocess.PIPE,
-                                             stderr=subprocess.PIPE)
+        # Create result directory
+        try:
+            os.makedirs(path)
+        except:
+            pass
-        # Write command arguments parameters
-        if len(args) > 0:
-            p.stdin.write(b'%s' % args)
-            p.stdin.close()
+        # Clear result directory
+        os.system('rm -rf %s/*' % path)
-        # Get results
-        res = p.communicate()
+        # Step in result directory
+        os.chdir(path)
-        # Set error flag and text
-        name = 'Test execution of command %s.' % cmd
-        self.test_assert((len(res[1]) == 0), name, res[1])
+        # Return
+        return
+
+    # Test command execution
+    def _test_execute_cmd(self, cmd, args):
+        """
+        Test command execution
+
+        Tests the execution of a command. Currently supported commands are:
+        - cd
+        - mkdir
+        - nano
+        - ./xxxxxx (any command xxxxxx that resides in the directory of the RST file)
+
+        Parameters
+        ----------
+        cmd : str
+            Command line
+        args : list of str
+            Additional lines following the command line
+        """
+        # Split command line
+        cmdline = cmd.split(' ')
+
+        # Check command line
+        self.test_assert((len(cmdline) > 0), 'Test command.', cmd)
+
+        # Continue only if valid
+        if len(cmdline) > 0:
+
+            # Extract command
+            command = gammalib.strip_whitespace(cmdline[0])
+
+            # Handle commands
+            if command == 'mkdir':
+                self.test_value(len(cmdline), 2, 'Test command "mkdir".')
+                self.test_value(len(args), 0, 'Test that no lines follow "mkdir".')
+                if len(cmdline) > 1:
+                    dirname = gammalib.strip_whitespace(cmdline[1])
+                    os.system('mkdir -p %s' % dirname)
+
+            elif command == 'cd':
+                self.test_value(len(cmdline), 2, 'Test command "cd".')
+                self.test_value(len(args), 0, 'Test that no lines follow "cd".')
+                if len(cmdline) > 1:
+                    dirname = gammalib.strip_whitespace(cmdline[1])
+                    os.chdir(dirname)
+
+            elif command == 'nano':
+                self.test_value(len(cmdline), 2, 'Test command "nano".')
+                self.test_assert((len(args) > 0), 'Test that lines follow "nano".')
+                if len(cmdline) > 1:
+                    fname = gammalib.strip_whitespace(cmdline[1])
+                    f = open(fname, 'wb')
+                    for arg in args:
+                        f.write(arg)
+                    f.close()
+
+            elif command[0:2] == './':
+                scriptname = command[2:]
+                self.test_assert((len(args) == 0), 'Test script "%s".' % scriptname)
+                tool = '%s/%s' % (self._rst_path, scriptname)
+
+                # Debug: show tool
+                #print tool.split(' '),
+
+                # Create subprocess for tool
+                p = subprocess.Popen(tool.split(' '), stdout=subprocess.PIPE,
+                                                      stdin=subprocess.PIPE,
+                                                      stderr=subprocess.PIPE)
+
+                # Get results
+                res = p.communicate()
+
+                # Set error flag and text
+                name = 'Test execution of %s.' % tool
+                self.test_assert((len(res[1]) == 0), name, res[1])
+
+            else:
+                name = 'Test command "%s".'   % command
+                msg  = 'Unknown command "%s"' % command
+                self.test_assert(False, name, msg)
         # Return
         return
@@ -113,8 +194,7 @@ class tutorials(gammalib.GPythonTestSuite):
                 # If we have a command then execute it
                 if cmd != '':
-                    print(cmd, args)
-                    self._execute_cmd(cmd, args)
+                    self._test_execute_cmd(cmd, args)
                 # Set new command
                 cmd  = line[2:]
@@ -126,8 +206,7 @@ class tutorials(gammalib.GPythonTestSuite):
         # If there is a pending command then execute it
         if cmd != '':
-            print(cmd, args)
-            self._execute_cmd(cmd, args)
+            self._test_execute_cmd(cmd, args)
         # Return
         return
@@ -149,23 +228,32 @@ class tutorials(gammalib.GPythonTestSuite):
         lines : list of str
             Command lines
         """
-        # Initialise command
-        cmd = ''
+        # Increment step counter
+        self._rst_step += 1
+        
+        # Initialise arguments
+        args = ''
+
+        # Set command
+        cmd = tool.split(' ')
         # Append lines
         for line in lines:
             pos = line.find(']')
             if pos != -1:
                 arg  = gammalib.strip_whitespace(line[pos+1:])
-                cmd += '%s\n' % arg
+                args += '%s\n' % arg
+
+        # Debug: show tool
+        #print cmd,
         # Create subprocess for tool
-        p = subprocess.Popen([tool], stdout=subprocess.PIPE,
-                                     stdin=subprocess.PIPE,
-                                     stderr=subprocess.PIPE)
+        p = subprocess.Popen(cmd, stdout=subprocess.PIPE,
+                                   stdin=subprocess.PIPE,
+                                   stderr=subprocess.PIPE)
         # Write command line parameters
-        p.stdin.write(b'%s' % cmd)
+        p.stdin.write(b'%s' % args)
         # Get results
         res = p.communicate()
@@ -177,13 +265,20 @@ class tutorials(gammalib.GPythonTestSuite):
         name = 'Test execution of %s.' % tool
         self.test_assert((len(res[1]) == 0), name, res[1])
+        # Preprend step counter to logfile
+        logfile = gammalib.GFilename('%s.log' % (cmd[0]))
+        if logfile.exists():
+            src = logfile.url()
+            dst = '%2.2d_%s' % (self._rst_step, src)
+            os.system('mv %s %s' % (src, dst))
+
         # Return
         return
-    # Add workflow step
-    def _add_workflow_step(self, tool, lines):
+    # Test workflow step
+    def _test_step(self, tool, lines):
         """
-        Add workflow step
+        Test workflow step
         Parameters
         ----------
@@ -191,18 +286,12 @@ class tutorials(gammalib.GPythonTestSuite):
             Tool name
         lines : list of str
             Command lines
-
-        Returns
-        -------
-        cmd : str
-            Workflow command
         """
-        # If tool is a ctool or cscript then add a ctool or cscript workflow
-        # step
+        # If tool is a ctool or cscript then test ctool or cscript execution
         if tool[0:2] == 'ct' or tool[0:2] == 'cs':
             self._test_ctool_execution(tool, lines)
-        # ... otherwise execute lines
+        # ... otherwise test line execution
         else:
             all_lines = ['$ %s' % tool]
             all_lines.extend(lines)
@@ -211,21 +300,33 @@ class tutorials(gammalib.GPythonTestSuite):
         # Return
         return
-    # Extract workflow
-    def _get_workflow(self, filename):
+    # Copy file
+    def _copy_file(self, filename):
         """
-        Extract workflow from Sphinx rst file
+        Copy file from rst directory into working directory
+        """
+        # Set source filename
+        src = self._rst_path + '/' + filename
+
+        # Copy file
+        os.system('cp %s %s' % (src, filename))
+
+        # Return
+        return
+
+    # Test Sphinx rst file
+    def _test_rst_file(self, filename):
+        """
+        Test Sphinx rst file
         Parameters
         ----------
         filename : str
-            RST filename
-
-        Returns
-        -------
-        workflow : `~gammalib.GXml`
-            Workflow XML file
+            Sphinx rst filename
         """
+        # Set RST path
+        self._rst_path = os.path.abspath(os.path.dirname(filename))
+
         # Open RST file
         f = open(filename, 'r')
@@ -238,6 +339,8 @@ class tutorials(gammalib.GPythonTestSuite):
             # If we are not in a code block then check whether a code block
             # starts
             if not code:
+
+                # Check for start of a code block
                 pos = line.find('code-block:: bash')
                 if pos != -1:
                     code  = True
@@ -246,6 +349,14 @@ class tutorials(gammalib.GPythonTestSuite):
                     lines = []
                     continue
+                # ... otherwise check for literalinclude, and if found, copy
+                # relevant file
+                pos = line.find('literalinclude::')
+                if pos != -1:
+                    filename = gammalib.strip_whitespace(line[pos+16:-1])
+                    self._copy_file(filename)
+                    continue
+
             # ... otherwise
             else:
@@ -257,7 +368,7 @@ class tutorials(gammalib.GPythonTestSuite):
                         tool  = line[pos+3:-1]
                         continue
-                # ... otherwise append lines or add workflow step
+                # ... otherwise append lines or test step
                 else:
                     blank = gammalib.strip_whitespace(line[0:start])
                     cmd   = line[start:-1]
@@ -265,13 +376,27 @@ class tutorials(gammalib.GPythonTestSuite):
                         lines.append(cmd)
                         continue
                     else:
-                        self._add_workflow_step(tool, lines)
+                        self._test_step(tool, lines)
                         code = False
                         continue
-        # Add step if some code is pending
+        # Test step if some code is pending
         if code:
-            self._add_workflow_step(tool, lines)
+            self._test_step(tool, lines)
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions
+        """
+        # Set test name
+        self.name('Tutorials Verification')
+
+        # Append tutorials
+        self.append(self.tutorials_1dc, 'Test 1DC tutorial')
         # Return
         return
@@ -281,15 +406,49 @@ class tutorials(gammalib.GPythonTestSuite):
         """
         Test 1DC tutorials
         """
-        # Set environment variables
-        os.environ['CTADATA'] = '/project-data/cta/data/1dc'
-        os.environ['CALDB']   = '/project-data/cta/data/1dc/caldb'
+        # Continue only if CTADATA1DC environment variable is set
+        if 'CTADATA1DC' in os.environ:
+        
+            # Set environment variables
+            os.environ['CTADATA'] = os.environ['CTADATA1DC']
+            os.environ['CALDB']   = os.environ['CTADATA1DC']+'/caldb'
+
+            # Set Sphinx rst file path
+            path = os.path.abspath('../doc/source/users/tutorials/1dc')
+
+            # Set result directory
+            self._set_result_dir('tutorials/1dc')
-        # Extract workflow
-        #workflow = self._get_workflow('../doc/source/users/tutorials/1dc/first_onoff.rst')
-        workflow = self._get_workflow('first_onoff.rst')
+            # Get current working directory
+            cwd = os.getcwd()
-        # Run workflow
+            # Clean pfiles
+            self._clean_pfiles()
+
+            # Initialise step counter
+            self._rst_step = 0
+
+            # Test Sphinx rst files
+            self._test_rst_file('%s/first_select_obs.rst' % path)
+            self._test_rst_file('%s/first_select_events.rst' % path)
+            self._test_rst_file('%s/first_skymap.rst' % path)
+            self._test_rst_file('%s/first_detect.rst' % path)
+            self._test_rst_file('%s/first_stacked.rst' % path)
+            self._test_rst_file('%s/first_fitting.rst' % path)
+            self._test_rst_file('%s/first_improving.rst' % path)
+            self._test_rst_file('%s/first_unbinned.rst' % path)
+
+            # Reset working directory
+            os.chdir(cwd)
+
+            # Clean pfiles
+            self._clean_pfiles()
+
+            # Initialise step counter
+            self._rst_step = 0
+
+            # Test Sphinx rst file for On/Off analysis
+            self._test_rst_file('%s/first_onoff.rst' % path)
         # Return
         return
@@ -318,15 +477,18 @@ if __name__ == &#39;__main__&#39;:
     except:
         pass
-    # Copy ctools parameter files into pfiles directory
-    os.system('cp -r ../src/*/*.par pfiles/')
-
     # Set PFILES environment variable
-    os.environ['PFILES'] = 'pfiles'
+    os.environ['PFILES'] = os.path.abspath('pfiles')
+
+    # Get current working directory
+    cwd = os.getcwd()
     # Run test suite
     success = suites.run()
+    # Reset working directory
+    os.chdir(cwd)
+
     # Save test results
     suites.save('reports/tutorials.xml')