Merge branch '3436-off-runs' into integration

Jürgen Knödlseder
2 parents 59d3205f 1898bd89
Showing 6 changed files with 395 additions and 74 deletions
cscripts/Makefile.am
cscripts/csphagen.par
cscripts/csphagen.py
doc/source/users/reference_manual/csphagen.rst
test/data/background_events.fits.gz
test/test_csphagen.py
@@ -182,6 +182,7 @@ 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/background_events.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 \
@@ -17,6 +17,7 @@
 # Input and output parameters
 #============================
 inobs,       f, a, obs.xml,,, "Input event list or observation definition XML file"
+inobsoff,    f, a, obs_off.xml,,, "Input Off event list or observation definition XML file"
 inmodel,     f, a, NONE,,, "Input model definition XML file (if NONE, use point source)"
 srcname,     s, a, Crab,,, "Source name"
 caldb,       s, a, prod2,,, "Calibration database"
@@ -46,10 +47,10 @@ width,         r, a, 0.2,0,180, &quot;Width of source region rectangle (deg)&quot;
 height,        r, a, 0.2,0,180, "Height of source region rectangle (deg)"
 posang,        r, a, 0.0,-360,360, "Position angle of source region rectangle (deg)"
 srcregfile,    f, a, on.reg,,, "Source region file (ds9 or FITS WCS)"
-bkgmethod,     s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation"
+bkgmethod,     s, a, REFLECTED,REFLECTED|CUSTOM|OFF,, "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"
-bkgregskip,    i, h, 1,0,10, "Number of background regions to skip next to On region"
+bkgregmin,     i, h, 2,1,1000, "Minimum number of background regions for REFLECTED method"
+bkgregskip,    i, h, 1,0,10, "Number of background regions to skip next to On region for REFLECTED method"
 use_model_bkg, b, h, yes,,, "Use background model?"
 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?"
@@ -3,7 +3,7 @@
 # Computes the PHA spectra for source/background and ARF/RMF files using the
 # reflected region method
 #
-# Copyright (C) 2017-2020 Luigi Tibaldo
+# Copyright (C) 2017-2021 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
@@ -43,6 +43,7 @@ class csphagen(ctools.csobservation):
         self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
         # Initialise other variables
+        self._obs_off       = gammalib.GObservations()
         self._ebounds       = gammalib.GEbounds()
         self._etruebounds   = gammalib.GEbounds()
         self._src_dir       = gammalib.GSkyDir()
@@ -74,6 +75,7 @@ class csphagen(ctools.csobservation):
         """
         # Set pickled dictionary
         state = {'base'          : ctools.csobservation.__getstate__(self),
+                 'obs_off'       : self._obs_off,
                  'ebounds'       : self._ebounds,
                  'etruebounds'   : self._etruebounds,
                  'src_dir'       : self._src_dir,
@@ -103,6 +105,7 @@ class csphagen(ctools.csobservation):
             Pickled instance
         """
         ctools.csobservation.__setstate__(self, state['base'])
+        self._obs_off       = state['obs_off']
         self._ebounds       = state['ebounds']
         self._etruebounds   = state['etruebounds']
         self._src_dir       = state['src_dir']
@@ -205,21 +208,17 @@ class csphagen(ctools.csobservation):
         # Return region container
         return regs
-    def _get_parameters_bkgmethod_reflected(self):
+    def _get_source_parameters(self):
         """
-        Get parameters for REFLECTED background method
+        Get parameters to define source/On region
         """
+
         # Get source shape
         self._srcshape = self['srcshape'].string()
         # Query source direction
         self._query_src_direction()
-        # Query minimum number of background regions and number of
-        # background regions to skip next to On region
-        self['bkgregmin'].integer()
-        self['bkgregskip'].integer()
-
         # If source shape is a circle the append GSkyRegionCircle
         if self._srcshape == 'CIRCLE':
@@ -240,8 +239,21 @@ class csphagen(ctools.csobservation):
                                                               self._reg_height,
                                                               self._reg_posang))
-        # Query usage of background model
-        self['use_model_bkg'].boolean()
+        # Return
+        return
+
+    def _get_parameters_bkgmethod_reflected(self):
+        """
+        Get parameters for REFLECTED background method
+        """
+
+        # Query parameters for source/On region definition
+        self._get_source_parameters()
+
+        # Query minimum number of background regions and
+        # number of background regions to skip next to On region
+        self['bkgregmin'].integer()
+        self['bkgregskip'].integer()
         # Return
         return
@@ -284,6 +296,76 @@ class csphagen(ctools.csobservation):
         # Return
         return
+    def _get_parameters_bkgmethod_off(self):
+        """
+        Get parameters for OFF background method
+
+        Raises
+        ------
+        RuntimeError
+            On and Off observations must have same size
+        RuntimeError
+            Off observations must be event lists
+        """
+
+        # Set up Off observations. If there are no Off observations in the
+        # container then load them via user parameters
+        if self.obs_off().is_empty():
+
+            # Get Off observation file name
+            filename = self['inobsoff'].filename()
+
+            # If Off observation is a FITS file then load observation and
+            # append it to the Off observation container
+            if gammalib.GFilename(filename).is_fits():
+                self._obs_off.append(gammalib.GCTAObservation(filename))
+
+            # ... otherwise load XML file into Off observation container
+            else:
+                self._obs_off.load(filename)
+
+        # Check that size of On and Off observations are the same, otherwise
+        # raise error
+        if self.obs().size() != self.obs_off().size():
+            raise RuntimeError('On and Off observations must have the same size')
+
+        # Loop through observations
+        for obs in self.obs_off():
+
+            # Check that observation is event list, otherwise throw error
+            if obs.eventtype() != "EventList":
+                raise RuntimeError('Off observations must be event lists')
+
+            # Check that they have response, otherwise assign based on user parameter
+            if obs.has_response() == False:
+
+                # Get database and IRF
+                database = self["caldb"].string()
+                irf      = self["irf"].string()
+
+                # Create an XML element for response
+                parameter = "parameter name=\"Calibration\"" +\
+                               " database=\"" + database + "\"" +\
+                               " response=\"" + irf + "\""
+                xml = gammalib.GXmlElement()
+                xml.append(parameter)
+
+                # Create CTA response
+                response = gammalib.GCTAResponseIrf(xml)
+
+                # Attach response to observation
+                obs.response(response)
+
+        # Add models from Off observations to model container
+        for model in self.obs_off().models():
+            self._models.append(model)
+
+        # Query parameters for source/On region definition
+        self._get_source_parameters()
+
+        # Return
+        return
+
     def _get_parameters_bkgmethod(self):
         """
         Get background method parameters
@@ -296,9 +378,12 @@ class csphagen(ctools.csobservation):
             self._get_parameters_bkgmethod_reflected()
         elif bkgmethod == 'CUSTOM':
             self._get_parameters_bkgmethod_custom()
+        elif bkgmethod == 'OFF':
+            self._get_parameters_bkgmethod_off()
         # Query parameters that are needed for all background methods
         self['maxoffset'].real()
+        self['use_model_bkg'].boolean()
         # Return
         return
@@ -351,9 +436,13 @@ class csphagen(ctools.csobservation):
         else:
             self._has_exclusion = False
-        # Query remaining parameters
-        self['use_model_bkg'].boolean()
-        self['stack'].boolean()
+        # Get background method parameters (have to come after setting up of
+        # observations and models)
+        self._get_parameters_bkgmethod()
+
+        # If there are multiple observations query whether to stack them
+        if self.obs().size() > 1:
+            self['stack'].boolean()
         # Query ahead output parameters
         if (self._read_ahead()):
@@ -361,15 +450,6 @@ class csphagen(ctools.csobservation):
             self['outmodel'].filename()
             self['prefix'].string()
-        # If there are no observations in container then get them from the
-        # parameter file
-        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)
@@ -459,14 +539,6 @@ class csphagen(ctools.csobservation):
             if self._srcshape == 'CIRCLE':
                 msg += ' (circle rad=%.3f).' % (self._rad)
             self._log_string(gammalib.NORMAL, msg)
-
-        # Skip observation if it is pointed too far from the source
-        elif offset >= self['maxoffset'].real():
-            msg = ' Skip because observation is pointed at %.3f deg >= '\
-                  '"maxoffset=%.3f" from source.' \
-                  % (offset, self['maxoffset'].real())
-            self._log_string(gammalib.NORMAL, msg)
-
         # ... otherwise
         else:
             posang = pnt_dir.posang_deg(self._src_dir)
@@ -474,7 +546,6 @@ class csphagen(ctools.csobservation):
                 # Determine number of background regions to skip
                 N_skip  = self['bkgregskip'].integer()
-                N_start = 1 + N_skip
                 N_lim   = 1 + 2 * N_skip
                 # Compute the angular separation of reflected regions wrt
@@ -532,13 +603,82 @@ class csphagen(ctools.csobservation):
                             regions.append(region)
                             dphi += alpha
-                    # Log number of reflected regions
-                    msg = ' Use %d reflected regions.' % (regions.size())
-                    self._log_string(gammalib.NORMAL, msg)
+                    # Check again that we have enough background regions
+                    # now that we have checked for overlap with exclusion region
+                    if regions.size() >= self['bkgregmin'].integer():
+                        # Log number of reflected regions
+                        msg = ' Use %d reflected regions.' % (regions.size())
+                        self._log_string(gammalib.NORMAL, msg)
+                    # Otherwise log observation skipped and return empty region container
+                    else:
+                        msg = ' Skip because the number %d of regions' \
+                              'for background estimation not overlapping ' \
+                              'with the exclusion region is smaller than ' \
+                              '"bkgregmin"=%d.' % (bkg_reg.size(), self['bkgregmin'].integer())
+                        self._log_string(gammalib.NORMAL, msg)
+                        regions = gammalib.GSkyRegions()
         # Return reflected regions
         return regions
+    def _instrument_regions(self, obs, obs_off):
+        """
+        Compute background regions for Off observation
+        
+        Calculate background region in Off observation that corresponds to the
+        source region in the On observation in instrument coordinates
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCTAObservation()`
+            On CTA observation
+        obs_off : `~gammalib.GCTAObservation()`
+            Off CTA observation
+
+        Returns
+        -------
+        regions : `~gammalib.GSkyRegions`
+            Container with background region
+        """
+        # Initialise region container
+        regions = gammalib.GSkyRegions()
+
+        # Convert source position in On observation to instrument coordinates
+        instdir = obs.pointing().instdir(self._src_dir)
+
+        # Convert instrument position to sky direction for Off observation
+        off_dir = obs_off.pointing().skydir(instdir)
+
+        # Build region according to shape specified by user
+        # If circle
+        if self._srcshape == 'CIRCLE':
+            region = gammalib.GSkyRegionCircle(off_dir, self._rad)
+
+        # ... otherwise if rectangle
+        elif self._srcshape == 'RECT':
+            # Instrument coordinates take sky direction as reference
+            # so no need to change the position angle
+            region = gammalib.GSkyRegionRectangle(off_dir,
+                                                  self._reg_width,
+                                                  self._reg_height,
+                                                  self._reg_posang)
+
+        # Check if background region overlaps with exclusion region
+        is_valid = True
+        if self._has_exclusion:
+            if self._excl_reg.overlaps(region):
+                # Signal region overlap
+                msg = ' Background region overlaps with exclusion region.'
+                self._log_string(gammalib.EXPLICIT, msg)
+                is_valid = False
+
+        # If region is valid then append it to container
+        if is_valid:
+            regions.append(region)
+
+        # Return
+        return regions
+
     def _set_models(self, results):
         """
         Set models for On/Off fitting
@@ -587,7 +727,7 @@ class csphagen(ctools.csobservation):
                 # observations
                 for result in results:
                     if model.is_valid(result['instrument'], result['id']):
-                        if result['bkg_reg'].size() >= self['bkgregmin'].integer():
+                        if result['bkg_reg'].size() > 0:
                             use_model = True
                         break
@@ -699,7 +839,7 @@ class csphagen(ctools.csobservation):
         # Return energy boundaries
         return ebounds
-    def _set_background_regions(self, obs):
+    def _set_background_regions(self, obs, obs_off=None):
         """
         Set background regions for an observation
@@ -716,7 +856,7 @@ class csphagen(ctools.csobservation):
         # Initialise empty background regions for this observation
         bkg_reg = gammalib.GSkyRegions()
-        # If reflected background is requested then created reflected
+        # If reflected background is requested then create reflected
         # background regions
         if self['bkgmethod'].string() == 'REFLECTED':
             bkg_reg = self._reflected_regions(obs)
@@ -728,6 +868,11 @@ class csphagen(ctools.csobservation):
         elif self['bkgmethod'].string() == 'CUSTOM':
             bkg_reg = obs.off_regions().copy()
+        # ... otherwise if dedicated Off runs are use then
+        # use background region that correspond to the same instrument coordinates
+        if self['bkgmethod'].string() == 'OFF':
+            bkg_reg = self._instrument_regions(obs,obs_off)
+
         # Return background regions
         return bkg_reg
@@ -750,6 +895,13 @@ class csphagen(ctools.csobservation):
         bkg_reg = None
         obs     = self.obs()[i]
+        # Retrieve dedicated Off observation if it exists
+        if not self.obs_off().is_empty():
+            obs_off = self.obs_off()[i]
+        # Otherwise use the same as On
+        else:
+            obs_off = self.obs()[i]
+
         # Log header
         self._log_header3(gammalib.NORMAL,'%s observation "%s"' % \
                           (obs.instrument(), obs.id()))
@@ -768,32 +920,43 @@ class csphagen(ctools.csobservation):
                 msg = ' Background model not used, assume constant backround rate.'
             self._log_string(gammalib.NORMAL, msg)
-            # Set background regions for this observation
-            bkg_reg = self._set_background_regions(obs)
+            # Get offset angle of source
+            pnt_dir = obs.pointing().dir()
+            offset = pnt_dir.dist_deg(self._src_dir)
-            # If there are background regions then create On/Off observation
-            # and append it to the output container
-            if bkg_reg.size() >= self['bkgregmin'].integer():
+            # Skip observation if it is pointed too far from the source
+            if offset >= self['maxoffset'].real():
+                msg = ' Skip because observation is pointed at %.3f deg >= ' \
+                      '"maxoffset=%.3f" from source.' \
+                      % (offset, self['maxoffset'].real())
+                self._log_string(gammalib.NORMAL, msg)
+            # ... otherwise continue to process
+            else:
-                # Create On/Off observation
-                onoff = gammalib.GCTAOnOffObservation(obs,
-                                                      self._models,
-                                                      self._srcname,
-                                                      self._etruebounds,
-                                                      self._ebounds,
-                                                      self._src_reg,
-                                                      bkg_reg,
-                                                      use_model_bkg)
+                # Set background regions for this observation
+                bkg_reg = self._set_background_regions(obs,obs_off)
-                # Set On/Off observation ID
-                onoff.id(obs.id())
+                # If there are any background regions then create On/Off observation
+                # and append it to the output container
+                if bkg_reg.size() >= 0:
-            # Otherwise log observation skipped
-            else:
-                msg = ' Skip because the number %d of regions ' \
-                      'for background estimation is smaller than ' \
-                      '"bkgregmin"=%d.' % (bkg_reg.size(), self['bkgregmin'].integer())
-                self._log_string(gammalib.NORMAL, msg)
+                    # Create On/Off observation
+                    onoff = gammalib.GCTAOnOffObservation(obs, obs_off,
+                                                          self._models,
+                                                          self._srcname,
+                                                          self._etruebounds,
+                                                          self._ebounds,
+                                                          self._src_reg,
+                                                          bkg_reg,
+                                                          use_model_bkg)
+
+                    # Set On/Off observation ID
+                    onoff.id(obs.id())
+
+                # Otherwise log observation skipped
+                else:
+                    msg = ' Skip because no valid Off regions could be determined'
+                    self._log_string(gammalib.NORMAL, msg)
         # Construct dictionary with results
         result = {'onoff'     : onoff,
@@ -854,6 +1017,8 @@ class csphagen(ctools.csobservation):
         # Write observation into logger
         self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
+        if not self.obs_off().is_empty():
+            self._log_observations(gammalib.NORMAL, self._obs_off, 'Off Observation')
         # Set true energy bins
         self._etruebounds = self._etrue_ebounds()
@@ -1032,6 +1197,27 @@ class csphagen(ctools.csobservation):
         # Return
         return self._excl_reg
+    def obs_off(self, obs=None):
+        """
+        Return and optionally set the Off observations
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCTAObservations`
+            Off observations container
+
+        Returns
+        -------
+        observation container : `~gammalib.GCTAObservations`
+            Off observations container
+        """
+        # If an observation container is provided then set the Off observations ...
+        if obs is not None:
+            self._obs_off = obs
+
+        # Return
+        return self._obs_off
+
 # ======================== #
 # Main routine entry point #
@@ -17,17 +17,16 @@ the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF),
 If the ``bkgmethod`` parameter is ``REFLECTED``, the script will use an On region
 shape defined by the ``srcshape`` parameter and generate reflected regions by
 placing regions of the same shape that are reflected with respect to the
-pointing direction. So far circular and rectangular regions are supported.
-Circular regions are defined by a center that is specified by the ``ra`` and
-``dec`` parameters if ``coordsys=CEL`` or the ``glon`` and ``glat`` parameters if
+pointing direction within the same observation. So far circular and
+rectangular regions are supported. Circular regions are defined by a
+center that is specified by the ``ra`` and ``dec`` parameters if
+``coordsys=CEL`` or the ``glon`` and ``glat`` parameters if
 ``coordsys=GAL``, and a radius that is specified by the ``rad`` parameter.
 Rectangular regions are defined by a center, again specified by either the
 ``ra`` and ``dec`` or the ``glon`` and ``glat`` parameters, a ``width`` parameter,
 a ``height`` parameter and a ``posang`` parameter that specified the position
 angle measured counter clockwise from celestial North. Note that the ``height``
-axis is pointing in the direction of the position angle. Reflected regions are
-only generated for offset angles with respect to the pointing direction that are
-smaller than the value specified by the ``maxoffset`` parameter. If less than
+axis is pointing in the direction of the position angle. If less than
 ``bkgregmin`` reflected regions can be found because the offset angle between
 source and pointing direction is too small, the observation will be skipped.
 Regions to be excluded for background determination can be specified using a
@@ -44,7 +43,7 @@ respectively. Both parameters accept either a ds9 region file or a sky map where
 all non-zero pixel values define a region. In case that an observation definition
 XML file is specified as input using the ``inobs``, background region files can
 be specified separately for each observation using the ``OffRegions`` parameter
-in the format
+in the format.
 .. code-block:: xml
@@ -53,6 +52,16 @@ in the format
      ...
    </observation>
+Finally,  if ``bkgmethod=OFF`` is specified, the background regions
+will be taken from independent Off observations given by the
+``inobsoff`` parameter. The On region can be defined in the same way
+as for ``bkgmethod=REFLECTED``, and the Off regions are the
+corresponding regions in instrument coordinates from the Off observations.
+
+For all background estimation methods only observations with a
+pointing direction at an offset smaller than ``maxoffset`` from the
+source are used. 
+
 :ref:`csphagen` supports multiprocessing. By default each event list will be treated
 in parallel over as many processes as the number of CPUs available on your
 machine. The maximum number of parallel processes can be set by the user through the
@@ -78,6 +87,11 @@ General parameters
 ``inobs [file]``
     Input event list or observation definition XML file.
+``inobsoff [file]``
+    Input event list or observation definition XML file for dedicated
+    Off observations. The number of observations must be the same as
+    for ``inobs``.
+
 ``inmodel [file]``
     Input model XML file (if ``NONE`` a point source at the centre of the
     source region is used).
@@ -172,7 +186,7 @@ General parameters
 ``srcregfile [file]``
     Source region file (ds9 or FITS WCS map).
-``bkgmethod <REFLECTED|CUSTOM> [string]``
+``bkgmethod <REFLECTED|CUSTOM|OFF> [string]``
     Method for background estimation:
     - ``REFLECTED``: background is evaluated in regions with the same shape as
@@ -184,18 +198,25 @@ General parameters
       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.
+      XML file. Off region files can be either ds9 region files or
+      FITS WCS maps.
+
+      - ``OFF``: background is evaluated from dedicated Off
+	observations in regions that correspond to the source region
+	in the On observations in instrument coordinates.
 ``bkgregfile [file]``
     Background regions file (ds9 or FITS WCS map).
 ``(bkgregmin = 2) [integer]``
-    Minimum number of background regions that are required for an observation.
+    Minimum number of background regions that are required for an
+    observation when the ``REFLECTED`` background estimation is used.
     If this number of background regions is not available the observation is
     skipped.
 ``(bkgregskip = 1) [integer]``
-    Number of background regions that should be skipped next to the On regions.
+    Number of background regions that should be skipped next to the On
+    region when the ``REFLECTED`` background estimation is used.
     Typically, one region is skipped so that the Off regions are taken sufficiently
     distant from the On region, but in some cases it may be useful to keep the
     background regions next to the On region.
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csphagen script.
 #
-# Copyright (C) 2017-2020 Luigi Tibaldo
+# Copyright (C) 2017-2021 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
@@ -46,6 +46,7 @@ class Test(test):
         self._myevents1      = self._datadir + '/crab_offaxis1.fits'
         self._myevents2      = self._datadir + '/crab_offaxis2.fits'
         self._exclusion      = self._datadir + '/crab_exclusion.fits'
+        self._offevents      = self._datadir + '/background_events.fits'
         self._nreg_with_excl = 5
         self._nreg_wo_excl   = 8
         self._nreg_bkg_reg   = 5
@@ -414,6 +415,117 @@ class Test(test):
         self._check_output('csphagen_py8', self._nbins, 9)
         self._check_outobs('csphagen_py8_obs.xml', 1)
+        # Setup csphagen for test with dedicated Off observations
+        phagen = cscripts.csphagen()
+        phagen['inobs']      = self._myevents1
+        phagen['inobsoff']   = self._offevents
+        phagen['inmodel']    = 'NONE'
+        phagen['caldb']      = self._caldb
+        phagen['irf']        = self._irf
+        phagen['srcshape']   = 'CIRCLE'
+        phagen['coordsys']   = 'CEL'
+        phagen['ra']         = 83.633
+        phagen['dec']        = 22.0145
+        phagen['rad']        = 0.2
+        phagen['ebinalg']    = 'LOG'
+        phagen['emin']       = 0.1
+        phagen['emax']       = 100.0
+        phagen['enumbins']   = self._nbins
+        phagen['bkgmethod']  = 'OFF'
+        phagen['etruemin']   = 0.05
+        phagen['etruemax']   = 150.0
+        phagen['etruebins']  = 5
+        phagen['stack']      = False
+        phagen['outobs']     = 'csphagen_py9_obs.xml'
+        phagen['outmodel']   = 'csphagen_py9_model.xml'
+        phagen['prefix']     = 'csphagen_py9'
+        phagen['logfile']    = 'csphagen_py9.log'
+        phagen['chatter']    = 2
+
+        # Execute script
+        phagen.execute()
+
+        # Check output
+        self._check_output('csphagen_py9', self._nbins, 1)
+        self._check_outobs('csphagen_py9_obs.xml', 1)
+
+        # Setup csphagen for test with dedicated Off observations and
+        # exclusion region
+        phagen = cscripts.csphagen()
+        phagen['inobs']       = self._myevents1
+        phagen['inobsoff']    = self._offevents
+        phagen['inmodel']     = 'NONE'
+        phagen['caldb']       = self._caldb
+        phagen['irf']         = self._irf
+        phagen['srcshape']    = 'CIRCLE'
+        phagen['coordsys']    = 'CEL'
+        phagen['ra']          = 83.633
+        phagen['dec']         = 22.0145
+        phagen['rad']         = 0.2
+        phagen['ebinalg']     = 'LOG'
+        phagen['emin']        = 0.1
+        phagen['emax']        = 100.0
+        phagen['enumbins']    = self._nbins
+        phagen['inexclusion'] = self._exclusion
+        phagen['bkgmethod']   = 'OFF'
+        phagen['etruemin']    = 0.05
+        phagen['etruemax']    = 150.0
+        phagen['etruebins']   = 5
+        phagen['stack']       = False
+        phagen['outobs']      = 'csphagen_py10_obs.xml'
+        phagen['outmodel']    = 'csphagen_py10_model.xml'
+        phagen['prefix']      = 'csphagen_py10'
+        phagen['logfile']     = 'csphagen_py10.log'
+        phagen['chatter']     = 2
+
+        # Execute script
+        phagen.execute()
+
+        # Check output
+        self._check_output('csphagen_py10', self._nbins, 1)
+        self._check_outobs('csphagen_py10_obs.xml', 1)
+
+        # Test with multiple input observations and dedicated Off runs
+        # Create Off observation container
+        obs_off = gammalib.GObservations()
+        for s, events in enumerate([self._offevents, self._offevents]):
+            run = gammalib.GCTAObservation(events)
+            run.id(str(s + 1))
+            run.response(self._irf, gammalib.GCaldb('cta', self._caldb))
+            obs_off.append(run)
+
+        # Setup csphagen
+        phagen = cscripts.csphagen(obs)
+        phagen.obs_off(obs_off)
+        phagen['inmodel']    = 'NONE'
+        phagen['srcshape']   = 'CIRCLE'
+        phagen['coordsys']   = 'CEL'
+        phagen['ra']         = 83.633
+        phagen['dec']        = 22.0145
+        phagen['rad']        = 0.2
+        phagen['ebinalg']    = 'LOG'
+        phagen['emin']       = 0.1
+        phagen['emax']       = 100.0
+        phagen['enumbins']   = self._nbins
+        phagen['bkgmethod']  = 'OFF'
+        phagen['etruemin']   = 0.05
+        phagen['etruemax']   = 150.0
+        phagen['etruebins']  = 5
+        phagen['stack']      = False
+        phagen['outobs']     = 'csphagen_py11_obs.xml'
+        phagen['outmodel']   = 'csphagen_py11_model.xml'
+        phagen['prefix']     = 'csphagen_py11'
+        phagen['logfile']    = 'csphagen_py11.log'
+        phagen['chatter']    = 2
+
+        # Execute script
+        phagen.execute()
+
+        # Check output
+        self._check_output('csphagen_py11_1', self._nbins, 1)
+        self._check_output('csphagen_py11_2', self._nbins, 1)
+        self._check_outobs('csphagen_py11_obs.xml', 2)
+
         # Return
         return