Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents e86c3204 cba4f988
Showing 12 changed files with 328 additions and 86 deletions
AUTHORS
ChangeLog
NEWS
cscripts/Makefile.am
cscripts/csphagen.par
cscripts/csphagen.py
cscripts/obsutils.py
doc/source/admin/release_history/1.5.rst
doc/source/users/reference_manual/csphagen.rst
test/data/crab_bkg_reg.reg
test/data/crab_src_reg.reg
test/test_csphagen.py
@@ -15,5 +15,6 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
 Stefan Ohm (stefan.ohm81@gmail.com)
 David Sanchez (david.sanchez@lapp.in2p3.fr)
 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
+Andreas Specovius (andreas.specovius@fau.de)
 Luigi Tibaldo (luigi.tibaldo@irap.omp.eu)
 Alexander Ziegler (alexander.ziegler@fau.de)
-2018-01-09
+2018-01-12
         * Version 1.5.0 released
         ========================
+        Add CUSTOM background method to csphagen script (#2284)
         Add csresspec script (#2289)
         Add ctools::time_reference constant and remove ctools::time_reference() method (#2295)
         Correctly strip .gz suffix from output filenames (#2287)
-New Features and Important Changes in ctools 1.5.0
-
-9 January 2018
+12 January 2018
 Introduction
@@ -228,6 +226,11 @@ None
 csphagen - Generate On/Off observations
 ---------------------------------------
+The CUSTOM background method was added to allow specifying On and Off regions
+from a region file or FITS WCS map. For multiple observations, a custom Off
+region per observation can be specified in the observation definition XML file
+using the OffRegions parameter (#2284).
+
 The script was added to allow the generation of PHA, ARF and RMF datasets for classical On/Off IACT analysis (#2208).
 #############################################################################
 # Makefile for cscripts scripts                                             #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2016-2017 Juergen Knoedlseder                               #
+# Copyright (C) 2016-2018 Juergen Knoedlseder                               #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
@@ -172,6 +172,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/crab_offaxis1.fits.gz \
            $(top_srcdir)/test/data/crab_offaxis2.fits.gz \
            $(top_srcdir)/test/data/crab_exclusion.fits.gz \
+           $(top_srcdir)/test/data/crab_src_reg.reg \
+           $(top_srcdir)/test/data/crab_bkg_reg.reg \
            $(top_srcdir)/test/data/tsmap_0.fits \
            $(top_srcdir)/test/data/tsmap_1.fits \
            $(top_srcdir)/test/data/exclusion.reg \
@@ -26,25 +26,27 @@ prefix,      s, h, onoff,,, &quot;Prefix for output files&quot;
 #
 # Script parameters
 #==================
-emin,      r, a, 0.1,,, "Lower energy limit (TeV)"
-emax,      r, a, 100.0,,, "Upper energy limit (TeV)"
-enumbins,  i, a, 120,,, "Number of energy bins"
-ebinalg,   s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
-ebinfile,  f, a, NONE,,, "Name of the file containing the energy bin definition"
-srcshape,  s, h, CIRCLE,CIRCLE,, "Shape of the source region"
-coordsys,  s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
-ra,        r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
-dec,       r, a, 22.01,-90,90, "Declination of source region centre (deg)"
-glon,      r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
-glat,      r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
-rad,       r, a, 0.2,0,180, "Radius of source region circle (deg)"
-bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
-bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
-maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
-stack,     b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?"
-etruemin,  r, h, 0.01,,, "Minimum true energy (TeV)"
-etruemax,  r, h, 300.0,,, "Maximum true energy (TeV)"
-etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
+emin,       r, a, 0.1,,, "Lower energy limit (TeV)"
+emax,       r, a, 100.0,,, "Upper energy limit (TeV)"
+enumbins,   i, a, 120,,, "Number of energy bins"
+ebinalg,    s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
+ebinfile,   f, a, NONE,,, "Name of the file containing the energy bin definition"
+srcshape,   s, h, CIRCLE,CIRCLE,, "Shape of the source region"
+coordsys,   s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
+ra,         r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
+dec,        r, a, 22.01,-90,90, "Declination of source region centre (deg)"
+glon,       r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
+glat,       r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
+rad,        r, a, 0.2,0,180, "Radius of source region circle (deg)"
+srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)"
+bkgmethod,  s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation"
+bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)"
+bkgregmin,  i, h, 2,1,1000, "Minimum number of background regions"
+maxoffset,  r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
+stack,      b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?"
+etruemin,   r, h, 0.01,,, "Minimum true energy (TeV)"
+etruemax,   r, h, 300.0,,, "Maximum true energy (TeV)"
+etruebins,  i, h, 30,,, "Number of bins per decade for true energy bins"
 #
 # Standard parameters
@@ -3,7 +3,7 @@
 # Computes the PHA spectra for source/background and ARF/RMF files using the
 # reflected region method
 #
-# Copyright (C) 2017 Luigi Tibaldo
+# Copyright (C) 2017-2018 Luigi Tibaldo
 #
 # 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
@@ -54,21 +54,14 @@ class csphagen(ctools.csobservation):
         # Return
         return
+
     # Private methods
-    def _get_parameters(self):
+    def _query_src_direction(self):
         """
-        Get parameters from parfile and setup observations
+        Set up the source direction parameter.
+        Query relevant parameters.
         """
-        # Setup observations (require response and allow event list, don't
-        # allow counts cube)
-        self._setup_observations(self.obs(), True, True, False)
-
-        # Set energy bounds
-        self._ebounds = self._create_ebounds()
-
-        # Initialise source position/region querying relevant parameters
         self._src_dir = gammalib.GSkyDir()
-        self._src_reg = gammalib.GSkyRegions()
         coordsys = self['coordsys'].string()
         if coordsys == 'CEL':
             ra  = self['ra'].real()
@@ -78,17 +71,100 @@ class csphagen(ctools.csobservation):
             glon = self['glon'].real()
             glat = self['glat'].real()
             self._src_dir.lb_deg(glon, glat)
+
+    def _get_parameters_bkgmethod_reflected(self):
+        """
+        Get parameters for REFLECTED background method
+        """
+        # Get source shape
         self._srcshape = self['srcshape'].string()
+
+        # Set source region (so far only CIRCLE is supported)
         if self._srcshape == 'CIRCLE':
+
+            # Query source direction
+            self._query_src_direction()
+
+            # Query minimum number of background regions
+            self['bkgregmin'].integer()
+
+            # Set circular source region
             self._rad = self['rad'].real()
             self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
-        # Query background estimation method and parameters
+        # Return
+        return
+
+    def _get_parameters_bkgmethod_custom(self):
+        """
+        Get parameters for CUSTOM background method
+
+        Raises
+        ------
+        RuntimeError
+            Only one On region is allowed
+        """
+        # Set up source region
+        filename      = self['srcregfile'].filename()
+        self._src_reg = gammalib.GSkyRegions(filename)
+
+        # Raise an exception if there is more than one source region
+        if len(self._src_reg) != 1:
+            raise RuntimeError('Only one On region is allowed')
+
+        # Set up source direction. Query parameters if neccessary.
+        if isinstance(self._src_reg[0], gammalib.GSkyRegionCircle):
+            self._src_dir = self._src_reg[0].centre()
+            self._rad     = self._src_reg[0].radius()
+        else:
+            self._query_src_direction()
+
+        # Make sure that all CTA observations have an Off region by loading the
+        # Off region region the parameter 'bkgregfile' for all CTA observations
+        # without Off region
+        for obs in self.obs():
+            if obs.instrument() == 'CTA':
+                if obs.off_regions().is_empty():
+                    filename = self['bkgregfile'].filename()
+                    regions  = gammalib.GSkyRegions(filename)
+                    obs.off_regions(regions)
+
+        # Return
+        return
+
+    def _get_parameters_bkgmethod(self):
+        """
+        Get background method parameters
+        """
+        # Get background method
         bkgmethod = self['bkgmethod'].string()
+
+        # Get background method dependent parameters
         if bkgmethod == 'REFLECTED':
-            self['bkgregmin'].integer()
+            self._get_parameters_bkgmethod_reflected()
+        elif bkgmethod == 'CUSTOM':
+            self._get_parameters_bkgmethod_custom()
+
+        # Query parameters that are needed for all background methods
         self['maxoffset'].real()
+        # Return
+        return
+
+    def _get_parameters(self):
+        """
+        Get parameters from parfile and setup observations
+        """
+        # Setup observations (require response and allow event list, don't
+        # allow counts cube)
+        self._setup_observations(self.obs(), True, True, False)
+
+        # Set energy bounds
+        self._ebounds = self._create_ebounds()
+
+        # Initialize empty src regions container
+        self._src_reg = gammalib.GSkyRegions()
+
         # Exclusion map
         if self['inexclusion'].filename().is_fits():
             self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename())
@@ -109,6 +185,10 @@ class csphagen(ctools.csobservation):
         if self.obs().is_empty():
             self.obs(self._get_observations(False))
+        # Get background method parameters (have to come after setting up of
+        # observations)
+        self._get_parameters_bkgmethod()
+
         # Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -247,6 +327,38 @@ class csphagen(ctools.csobservation):
         # Return energy boundaries
         return ebounds
+    def _set_background_regions(self, obs):
+        """
+        Set background regions for an observation
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCTAObservation()`
+            CTA observation
+
+        Returns
+        -------
+        regions : `~gammalib.GSkyRegions()`
+            Background regions
+        """
+        # Initialise empty background regions for this observation
+        bkg_reg = gammalib.GSkyRegions()
+
+        # If reflected background is requested then created reflected
+        # background regions
+        if self['bkgmethod'].string() == 'REFLECTED':
+            bkg_reg = self._reflected_regions(obs)
+
+        # ... otherwise if custom background is requested then get the
+        # background regions from the observation. We use a copy here since
+        # otherwise the background regions go out of scope once the observations
+        # are replaced by the On/Off observations.
+        elif self['bkgmethod'].string() == 'CUSTOM':
+            bkg_reg = obs.off_regions().copy()
+
+        # Return background regions
+        return bkg_reg
+
     # Public methods
     def run(self):
@@ -286,15 +398,20 @@ class csphagen(ctools.csobservation):
         # Loop through observations and generate pha, arf, rmf files
         for obs in self.obs():
-            # Initialise background regions for this observation
-            bkg_reg = gammalib.GSkyRegions()
+            # Skip non CTA observations
+            if obs.instrument() != 'CTA':
+                self._log_string(gammalib.NORMAL, 'Skip %s observation "%s"' % \
+                                 (obs.instrument(), obs.id()))
+                continue
+
+            # Log current observation
+            self._log_string(gammalib.NORMAL, 'Process observation "%s"' % \
+                             (obs.id()))
-            # If reflected background is requested then created reflected
-            # background regions
-            if self['bkgmethod'].string() == 'REFLECTED':
-                bkg_reg = self._reflected_regions(obs)
+            # Set background regions for this observation
+            bkg_reg = self._set_background_regions(obs)
-            # If there are reflected regions then create On/Off observation
+            # If there are background regions then create On/Off observation
             # and append it to the output container
             if bkg_reg.size() >= self['bkgregmin'].integer():
                 onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir,
@@ -130,6 +130,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
             phagen['rad']         = onrad
             phagen['stack']       = False
             phagen['inexclusion'] = 'NONE'
+            phagen['bkgmethod']   = 'REFLECTED'
             # Optionally open the log file
             if log:
@@ -29,6 +29,8 @@ Improvements
 * [`2295 <https://cta-redmine.irap.omp.eu/issues/2295>`_] -
   Add ``ctools::time_reference`` constant and remove ``ctools::time_reference()``
   method
+* [`2284 <https://cta-redmine.irap.omp.eu/issues/2284>`_] -
+  Add ``CUSTOM`` background method to :ref:`csphagen` script
 * [`2273 <https://cta-redmine.irap.omp.eu/issues/2273>`_] -
   Support On/Off analysis in :ref:`cslightcrv` and :ref:`csphasecrv`
 * [`2267 <https://cta-redmine.irap.omp.eu/issues/2267>`_] -
@@ -13,8 +13,9 @@ Synopsis
 This script can be used to derive from an observation or a set of observations
 the count spectra in a source region and in background regions, as well as the
 detector response (effective area, energy redistribution matrix), that can be
-used to perform a classical 1D spectral analysis. The output files are saved in
-the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF).
+used to perform a classical 1D spectral analysis. Regions can be placed
+automatically or by hand. The output files are saved in the OGIP format normally
+used in X-ray astronomy (PHA, ARF, RMF).
 `See here <https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html>`__.
 The outputs are:
@@ -30,50 +31,48 @@ General parameters
 ------------------
 ``inobs [file]``
-    Input event list or observation definition XML file.
+    Input event list or observation definition XML file
 ``caldb [string]``
-    Calibration database.
+    Calibration database
 ``irf [string]``
-    Instrument response function.
+    Instrument response function
 ``(inexclusion = NONE) [file]``
     Optional FITS file containing a WCS map in the first hdu that defines sky
-    regions not to be used for background estimation (where map value != 0).
+    regions not to be used for background estimation (where map value != 0)
 ``outobs [string]``
-    Output observation definition XML file.
+    Output observation definition XML file
 ``(prefix = onoff) [string]``
-    Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files.
+    Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files
 ``emin [real]``
-    Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used.
+    Lower energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used
 ``emax [real]``
-    Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used.
+    Upper energy limit (TeV) if ``LIN`` or ``LOG`` binning algorithms are used
 ``enumbins [integer]``
     Number of energy bins. At least 30 bins per decade are recommended for
     proper evaluation of the instrument response.
-``ebinalg [string]``
-    Energy binning algorithm:
-    ``LOG`` : logarithmically spaced energy bins;
-    ``LIN`` : linearly spaced energy bins;
-    ``FILE``: energy bounds retrieved from file.
+``ebinalg <FILE|LIN|LOG> [string]``
+    Algorithm for defining energy bins (``FILE``: energy bounds retrieved from
+    file, see ``ebinfile`` parameter; ``LIN``: linearly spaced energy bins;
+    ``LOG``: logarithmically spaced energy bins)
 ``ebinfile [file]``
     Name of the file containing the energy bin definition if ``FILE`` algorithm
-    is used.
+    is used
-``(srcshape = CIRCLE) [string]``
-    Shape of the source region.
-    ``CIRCLE``: circular region around given position.
+``(srcshape = CIRCLE) <CIRCLE> [string]``
+    Shape of the source region (``CIRCLE``: circular region around given position)
-``coordsys [string]``
-    Coordinate system (``CEL`` - celestial, ``GAL`` - galactic)
+``coordsys <CEL|GAL> [string]``
+    Coordinate system (``CEL``: celestial; ``GAL``: galactic)
 ``ra [real]``
     Right Ascension of source region centre (deg)
@@ -90,10 +89,25 @@ General parameters
 ``rad [real]``
     Radius of source region circle (deg)
-``(bkgmethod = REFLECTED) [string]``
-    Method for background estimation.
-    ``REFLECTED:`` background evaluated in regions with the same shape as
-    source region reflected w.r.t. pointing direction for each observation.
+``srcregfile [file]``
+    Source region file (ds9 or FITS WCS map)
+
+``bkgmethod <REFLECTED|CUSTOM> [string]``
+    Method for background estimation:
+
+    - ``REFLECTED``: background is evaluated in regions with the same shape as
+      the source region reflected w.r.t. pointing direction for each observation
+
+    - ``CUSTOM``: background is evaluated in regions specified by user. For an
+      event list or a single observation in the observation definition XML file
+      a region file will be queried (see ``bkgregfile`` parameter). For multiple
+      observations specified in the observation definition XML file the name of
+      the region file will be extracted from the ``OffRegions`` parameter that
+      needs to be specified for each observation in the observation definition
+      XML file. Off region files can be either ds9 region files or FITS WCS maps.
+
+``bkgregfile [file]``
+    Background regions file (ds9 or FITS WCS map)
 ``(bkgregmin = 2) [integer]``
     Minimum number of background regions that are required for an observation.
@@ -102,20 +116,20 @@ General parameters
 ``(maxoffset = 4.0) [real]``
     Maximum offset in degrees of source from camera center to accept the
-    observation.
+    observation
 ``stack [boolean]``
     Specifies whether multiple observations should be stacked (``yes``) or
-    whether run-wise PHA, ARF and RMF files should be produced (``no``).
+    whether run-wise PHA, ARF and RMF files should be produced (``no``)
 ``(etruemin = 0.01) [real]``
-    Minimum true energy (TeV).
+    Minimum true energy (TeV)
 ``(etruemax = 0.01) [real]``
-    Maximum true energy (TeV).
+    Maximum true energy (TeV)
 ``(etruebins = 30) [integer]``
-    Number of bins per decade for true energy bins.
+    Number of bins per decade for true energy bins
 Standard parameters
+# Region file format: DS9 version 4.1
+global color=green dashlist=8 3 width=1
+fk5;circle(82.2412081162728,22.7975361111674,0.2)
+fk5;circle(82.5360705305998,23.1293202741537,0.2)
+fk5;circle(82.980058098929,23.260640435668,0.2)
+fk5;circle(83.4304728956696,23.1492236356079,0.2)
+fk5;circle(83.7424610541293,22.8309432073396,0.2)
+# Region file format: DS9 version 4.1
+global color=green dashlist=8 3 width=1
+fk5;circle(83.633,22.0145,0.2)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csphagen script.
 #
-# Copyright (C) 2017 Luigi Tibaldo
+# Copyright (C) 2017-2018 Luigi Tibaldo
 #
 # 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
@@ -48,10 +48,11 @@ class Test(test):
         self._exclusion      = self._datadir + '/crab_exclusion.fits'
         self._nreg_with_excl = 5
         self._nreg_wo_excl   = 8
+        self._nreg_bkg_reg   = 5
         self._nreg_mul       = [self._nreg_with_excl, 6]
         self._nbins          = 10
-        # number of expected background regions with/wo exclusion,
-        # and for two different runs
+        self._regfile_src    = self._datadir + '/crab_src_reg.reg'
+        self._regfile_bkg    = self._datadir + '/crab_bkg_reg.reg'
         # Return
         return
@@ -88,7 +89,8 @@ class Test(test):
                          ' coordsys="CEL" ra=83.633 dec=22.0145' + \
                          ' rad=0.2 stack="no" inexclusion="' + \
                          self._exclusion + \
-                         '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
+                         '" bkgmethod="REFLECTED" '+ \
+                         'outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
                          'logfile="csphagen_cmd1.log" chatter=2'
         # Check if execution of wrong command fails
@@ -112,7 +114,8 @@ class Test(test):
                          ' coordsys="CEL" ra=83.633 dec=22.0145' + \
                          ' rad=0.2 stack="no" inexclusion="' + \
                          self._exclusion + \
-                         '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
+                         '" bkgmethod="REFLECTED" '+ \
+                         'outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
                          'logfile="csphagen_cmd2.log" debug=yes chatter=1'
         # Check if execution failed
@@ -144,6 +147,7 @@ class Test(test):
         phagen['rad']         = 0.2
         phagen['stack']       = False
         phagen['inexclusion'] = self._exclusion
+        phagen['bkgmethod']   = 'REFLECTED'
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
@@ -152,7 +156,7 @@ class Test(test):
         phagen['logfile']     = 'csphagen_py1.log'
         phagen['chatter']     = 1
-        # Run script
+        # Execute script
         phagen.execute()
         # Check output
@@ -173,6 +177,7 @@ class Test(test):
         phagen['dec']       = 22.0145
         phagen['rad']       = 0.2
         phagen['stack']     = False
+        phagen['bkgmethod'] = 'REFLECTED'
         phagen['etruemin']  = 0.05
         phagen['etruemax']  = 150.0
         phagen['etruebins'] = 5
@@ -181,7 +186,7 @@ class Test(test):
         phagen['logfile']   = 'csphagen_py2.log'
         phagen['chatter']   = 2
-        # Run script
+        # Execute script
         phagen.execute()
         # Check output
@@ -209,6 +214,7 @@ class Test(test):
         phagen['rad']         = 0.2
         phagen['stack']       = False
         phagen['inexclusion'] = self._exclusion
+        phagen['bkgmethod']   = 'REFLECTED'
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
@@ -226,8 +232,7 @@ class Test(test):
                                self._nreg_mul[s])
         self._check_outobs('csphagen_py3', 2)
-        # Setup csphagen for test with multiple input observations and
-        # stacking
+        # Setup csphagen for test with multiple input observations and stacking
         phagen = cscripts.csphagen(obs)
         phagen['ebinalg']     = 'LOG'
         phagen['emin']        = 0.1
@@ -239,6 +244,7 @@ class Test(test):
         phagen['rad']         = 0.2
         phagen['stack']       = True
         phagen['inexclusion'] = self._exclusion
+        phagen['bkgmethod']   = 'REFLECTED'
         phagen['etruemin']    = 0.05
         phagen['etruemax']    = 150.0
         phagen['etruebins']   = 5
@@ -247,7 +253,7 @@ class Test(test):
         phagen['logfile']     = 'csphagen_py4.log'
         phagen['chatter']     = 4
-        # Run script
+        # Execute script
         phagen.execute()
         # Check output
@@ -256,6 +262,66 @@ class Test(test):
                                0, check_regions=False)
         self._check_outobs('csphagen_py4', 1)
+        # Setup csphagen for test with custom On and Off regions provided
+        phagen = cscripts.csphagen()
+        phagen['inobs']      = self._myevents1
+        phagen['caldb']      = self._caldb
+        phagen['irf']        = self._irf
+        phagen['ebinalg']    = 'LOG'
+        phagen['emin']       = 0.1
+        phagen['emax']       = 100.0
+        phagen['enumbins']   = self._nbins
+        phagen['bkgmethod']  = 'CUSTOM'
+        phagen['srcregfile'] = self._regfile_src
+        phagen['bkgregfile'] = self._regfile_bkg
+        phagen['etruemin']   = 0.05
+        phagen['etruemax']   = 150.0
+        phagen['etruebins']  = 5
+        phagen['stack']      = False
+        phagen['outobs']     = 'csphagen_py5.xml'
+        phagen['prefix']     = 'csphagen_py5'
+        phagen['logfile']    = 'csphagen_py5.log'
+        phagen['chatter']    = 2
+
+        # Execute script
+        phagen.execute()
+
+        # Check output
+        self._check_output('csphagen_py5', self._nbins, self._nreg_bkg_reg)
+        self._check_outobs('csphagen_py5', 1)
+
+        # Append off regions to observation container
+        for run in obs:
+            run.off_regions(gammalib.GSkyRegions(self._regfile_bkg))
+
+        # Setup csphagen for test with multiple input observations and stacking
+        phagen = cscripts.csphagen(obs)
+        phagen['ebinalg']     = 'LOG'
+        phagen['emin']        = 0.1
+        phagen['emax']        = 100.0
+        phagen['enumbins']    = self._nbins
+        phagen['coordsys']    = 'CEL'
+        phagen['stack']       = True
+        phagen['inexclusion'] = self._exclusion
+        phagen['bkgmethod']   = 'CUSTOM'
+        phagen['srcregfile']  = self._regfile_src
+        phagen['etruemin']    = 0.05
+        phagen['etruemax']    = 150.0
+        phagen['etruebins']   = 5
+        phagen['outobs']      = 'csphagen_py6.xml'
+        phagen['prefix']      = 'csphagen_py6'
+        phagen['logfile']     = 'csphagen_py6.log'
+        phagen['chatter']     = 4
+
+        # Execute script
+        phagen.execute()
+
+        # Check output
+        for s in range(2):
+            self._check_output('csphagen_py6_stacked', self._nbins,
+                               0, check_regions=False)
+        self._check_outobs('csphagen_py6', 1)
+
         # Return
         return
@@ -309,7 +375,7 @@ class Test(test):
             self.test_assert(table.contains(col),
                              'FITS file contains "' + col + '" column')
-        # Check SPECTRUM table
+        # Check EBOUNDS table
         table = fits['EBOUNDS']
         self._check_ebounds(table, bins)
@@ -433,3 +499,26 @@ class Test(test):
         # Return
         return
+
+    def _read_pha_counts(self, filename):
+        """
+        Read and integrate the counts in a pha file.
+        Pha file structure already tested in _check_pha()
+        """
+        # Open FITS file
+        fits = gammalib.GFits(filename)
+
+        # Get SPECTRUM table
+        table = fits['SPECTRUM']
+
+        # Integrate counts
+        counts_col = table['COUNTS']
+        counts     = 0
+        for channel in range( counts_col.nrows() ):
+            counts += counts_col[channel]
+
+        # Close FITS file
+        fits.close()
+
+        # Return
+        return counts