diff --git a/AUTHORS b/AUTHORS
index daee14b..42793a0 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -15,4 +15,5 @@ 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)
+Luigi Tibaldo (luigi.tibaldo@irap.omp.eu)
Alexander Ziegler (alexander.ziegler@fau.de)
diff --git a/ChangeLog b/ChangeLog
index 2960fe1..4694694 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,9 @@
-2017-10-11
+2017-10-16
* Version 1.5.0 released
========================
+ Add csphagen script (#2208)
Add keywords to sky maps produced by ctskymap (#2194)
Enhance efficiency of computations in ctmodel (#2189)
Add csmodelsois script (#2180)
diff --git a/NEWS b/NEWS
index 49b1fc0..4238472 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
New Features and Important Changes in ctools 1.5.0
-11 October 2017
+16 October 2017
Introduction
@@ -132,6 +132,9 @@ a diffuse map cube (#2180).
Use t-type parameters for "tmin" and "tmax" to all specifying times in various
formats (#1864).
+The script csphagen was added to allow the generation of PHA, ARF and RMF datasets
+for classical On/Off IACT analysis (#2208).
+
Examples
--------
diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index d0bd9d1..381e166 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -33,6 +33,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
$(srcdir)/csobsdef.py \
$(srcdir)/csobsinfo.py \
$(srcdir)/csobsselect.py \
+ $(srcdir)/csphagen.py \
$(srcdir)/csphasecrv.py \
$(srcdir)/cspull.py \
$(srcdir)/csresmap.py \
@@ -66,6 +67,7 @@ parfiles = $(srcdir)/cscaldb.par \
$(srcdir)/csobsdef.par \
$(srcdir)/csobsinfo.par \
$(srcdir)/csobsselect.par \
+ $(srcdir)/csphagen.par \
$(srcdir)/csphasecrv.par \
$(srcdir)/cspull.par \
$(srcdir)/csresmap.par \
@@ -96,6 +98,7 @@ linkscripts = cscaldb \
csobsdef \
csobsinfo \
csobsselect \
+ csphagen \
csphasecrv \
cspull \
csresmap \
@@ -131,6 +134,7 @@ testscripts = $(srcdir)/tests/__init__.py \
$(top_srcdir)/test/test_csobsdef.py \
$(top_srcdir)/test/test_csobsinfo.py \
$(top_srcdir)/test/test_csobsselect.py \
+ $(top_srcdir)/test/test_csphagen.py \
$(top_srcdir)/test/test_csphasecrv.py \
$(top_srcdir)/test/test_cspull.py \
$(top_srcdir)/test/test_csresmap.py \
@@ -159,6 +163,9 @@ testdata = $(top_srcdir)/test/data/crab.xml \
$(top_srcdir)/test/data/crab_psfcube.fits.gz \
$(top_srcdir)/test/data/crab_edispcube.fits.gz \
$(top_srcdir)/test/data/crab_bkgcube.fits.gz \
+ $(top_srcdir)/test/data/crab_offaxis1.fits \
+ $(top_srcdir)/test/data/crab_offaxis2.fits \
+ $(top_srcdir)/test/data/crab_exclusion.fits \
$(top_srcdir)/test/data/tsmap_0.fits \
$(top_srcdir)/test/data/tsmap_1.fits \
$(top_srcdir)/test/data/exclusion.reg \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index 09f8982..a007255 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -36,6 +36,7 @@ __all__ = [
"csobsdef",
"csobsinfo",
"csobsselect",
+ "csphagen",
"csphasecrv",
"cspull",
"csresmap",
@@ -72,6 +73,7 @@ from cscripts.csobs2caldb import csobs2caldb
from cscripts.csobsdef import csobsdef
from cscripts.csobsinfo import csobsinfo
from cscripts.csobsselect import csobsselect
+from cscripts.csphagen import csphagen
from cscripts.csphasecrv import csphasecrv
from cscripts.cspull import cspull
from cscripts.csresmap import csresmap
diff --git a/cscripts/csobsselect.py b/cscripts/csobsselect.py
index 54f66e2..523e539 100755
--- a/cscripts/csobsselect.py
+++ b/cscripts/csobsselect.py
@@ -320,7 +320,7 @@ class csobsselect(ctools.cscript):
# Otherwise log filename and save file
else:
# Log filename
- self._log_value(gammalib.NORMAL, 'Observation definition XML file',
+ self._log_value(gammalib.NORMAL, 'Obs. definition XML file',
outobs.url())
# Save observations
diff --git a/cscripts/csphagen.par b/cscripts/csphagen.par
new file mode 100644
index 0000000..4cb633e
--- /dev/null
+++ b/cscripts/csphagen.par
@@ -0,0 +1,53 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: csphagen.par #
+# #
+# Author: Luigi Tibaldo #
+# luigi.tibaldo@irap.omp.eu #
+# IRAP #
+# #
+# Purpose: Parameter file for the csphagen script #
+# #
+#############################################################################
+#
+# Input and output parameters
+#============================
+inobs, f, a, obs.xml,,, "Input event list or observation definition XML file"
+caldb, s, a, prod2,,, "Calibration database"
+irf, s, a, South_0.5h,,, "Instrument response function"
+inexclusion, f, h, NONE,,, "Sky map of exclusion regions"
+outobs, s, a, onoff_obs.xml,,, "Output observation definition XML file"
+prefix, s, h, onoff,,, "Prefix for output files"
+
+#
+# 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?"
+
+#
+# Standard parameters
+#====================
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug, b, h, no,,, "Debugging mode activated"
+mode, s, h, ql,,, "Mode of automatic parameters"
+logfile, f, h, csphagen.log,,, "Log filename"
diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py
new file mode 100755
index 0000000..b35dc98
--- /dev/null
+++ b/cscripts/csphagen.py
@@ -0,0 +1,362 @@
+#! /usr/bin/env python
+# ==========================================================================
+# Computes the PHA spectra for source/background and ARF/RMF files using the
+# reflected region method
+#
+# Copyright (C) 2017 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import gammalib
+import ctools
+import math
+import sys
+
+
+# =============== #
+# csfindobs class #
+# =============== #
+class csphagen(ctools.cscript):
+ """
+ Generate PHA, ARF and RMF files for classical IACT spectral analysis
+ """
+
+ # Constructor
+ def __init__(self, *argv):
+ """
+ Constructor
+ """
+ # Set name, version
+ self._name = 'csphagen'
+ self._version = ctools.__version__
+
+ # Initialise observation container from constructor arguments
+ self._obs, argv = self._set_input_obs(argv)
+
+ # Initialise application by calling the appropriate class constructor
+ self._init_cscript(argv)
+
+ # Initialise other variables
+ self._outobs = gammalib.GObservations()
+ self._ebounds = gammalib.GEbounds()
+ self._src_dir = gammalib.GSkyDir()
+ self._src_reg = gammalib.GSkyRegions()
+ self._bkg_regs = []
+ self._excl_reg = None
+ self._has_exclusion = False
+
+ # Return
+ return
+
+ # Private methods
+ 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()
+
+ # 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()
+ dec = self['dec'].real()
+ self._src_dir.radec_deg(ra, dec)
+ elif coordsys == 'GAL':
+ glon = self['glon'].real()
+ glat = self['glat'].real()
+ self._src_dir.lb_deg(glon, glat)
+ self._srcshape = self['srcshape'].string()
+ if self._srcshape == 'CIRCLE':
+ self._rad = self['rad'].real()
+ self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
+
+ # Query background estimation method and parameters
+ bkgmethod = self['bkgmethod'].string()
+ if bkgmethod == 'REFLECTED':
+ self['bkgregmin'].integer()
+ self['maxoffset'].real()
+
+ # Exclusion map
+ if self['inexclusion'].filename().is_fits():
+ self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename())
+ self._has_exclusion = True
+ else:
+ self._has_exclusion = False
+
+ # Query remaining parametersStacking
+ self['stack'].boolean()
+
+ # Query ahead output parameters
+ if (self._read_ahead()):
+ self['outobs'].string()
+ self['prefix'].string()
+
+ # If there are no observations in container then get them from the
+ # parameter file
+ if self._obs.size() == 0:
+ self._obs = self._get_observations(False)
+
+ # Write input parameters into logger
+ self._log_parameters(gammalib.TERSE)
+
+ # Return
+ return
+
+ def _reflected_regions(self, obs):
+ """
+ Calculate list of reflected regions for a single observation (pointing)
+
+ Parameters
+ ----------
+ obs : `~gammalib.GObservation()`
+ CTA observation
+
+ Returns
+ -------
+ regions : `~gammalib.GSkyRegions`
+ List of reflected regions
+ """
+ # Initialise list of reflected regions
+ regions = gammalib.GSkyRegions()
+
+ # Get offset angle of source
+ pnt_dir = obs.pointing().dir()
+ offset = pnt_dir.dist_deg(self._src_dir)
+
+ # If ...
+ if offset <= self._rad or offset >= self['maxoffset'].real():
+ self._log_string(gammalib.EXPLICIT, 'Observation {} pointed at {} '
+ 'deg from source'.format(obs.id(), offset))
+
+ # ... otherwise
+ else:
+ posang = pnt_dir.posang_deg(self._src_dir)
+ if self._srcshape == 'CIRCLE':
+ # angular separation of reflected regions wrt camera center
+ # and number
+ alpha = 1.05 * 2 * self._rad / offset
+ # 1.05 ensures background regions do not overlap due to
+ # numerical precision issues
+ N = int(2.0 * math.pi / alpha)
+ if N < self['bkgregmin'].integer() + 3:
+ self._log_string(gammalib.EXPLICIT, 'Observation {}: '
+ 'insufficient regions for background '
+ 'estimation'.format(obs.id(), offset))
+
+ # Otherwise loop over position angle to create reflected
+ # regions
+ else:
+ alpha = 360.0 / N
+ for s in range(2, N - 1):
+ dphi = s * alpha
+ ctr_dir = pnt_dir.clone()
+ ctr_dir.rotate_deg(posang + dphi, offset)
+ region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
+ if self._has_exclusion:
+ if self._excl_reg.overlaps(region):
+ self._log_string(gammalib.VERBOSE, 'Observation '
+ '{}: reflected region overlaps '
+ 'with exclusion region'.
+ format(obs.id(), offset))
+ else:
+ regions.append(region)
+ else:
+ regions.append(region)
+
+ # Return reflected regions
+ return regions
+
+ # Public methods
+ def run(self):
+ """
+ Run the script
+ """
+ # Switch screen logging on in debug mode
+ if self._logDebug():
+ self._log.cout(True)
+
+ # Get parameters
+ self._get_parameters()
+
+ # Write observation into logger
+ self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
+
+ # Write header
+ self._log_header1(gammalib.NORMAL,
+ 'Generation of source and background spectra')
+
+ # Initialise run variables
+ self._outobs = gammalib.GObservations()
+ self._bkg_regs = []
+
+ # Loop through observations and generate pha, arf, rmf files
+ for obs in self._obs:
+
+ # Initialise 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)
+
+ # If there are reflected 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,
+ self._ebounds,
+ self._ebounds,
+ self._src_reg,
+ bkg_reg)
+ onoff.id(obs.id())
+ self._outobs.append(onoff)
+ self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()})
+ else:
+ self._log_string(gammalib.NORMAL, 'Observation {} not included '
+ 'in spectra generation'.format(obs.id()))
+
+ # Stack observations
+ if self._outobs.size() > 1 and self['stack'].boolean() == True:
+
+ # Write header
+ self._log_header1(gammalib.NORMAL,
+ 'Stacking {} observations'.format(self._outobs.size()))
+
+ # Stack observations
+ stacked_obs = gammalib.GCTAOnOffObservation(self._outobs)
+
+ # Put stacked observations in output container
+ self._outobs = gammalib.GObservations()
+ self._outobs.append(stacked_obs)
+
+ # Return
+ return
+
+ def save(self):
+ """
+ Save data
+ """
+ # Write header
+ self._log_header1(gammalib.TERSE, 'Save data')
+
+ # Get XMK output filename, prefix and clobber
+ outobs = self['outobs'].string()
+ prefix = self['prefix'].string()
+ clobber = self['clobber'].boolean()
+
+ # Loop over all observation in container
+ for obs in self._outobs:
+
+ # Set filenames
+ if self['stack'].boolean():
+ onname = prefix + '_stacked_pha_on.fits'
+ offname = prefix + '_stacked_pha_off.fits'
+ arfname = prefix + '_stacked_arf.fits'
+ rmfname = prefix + '_stacked_rmf.fits'
+ elif self._outobs.size() > 1:
+ onname = prefix + '_{}_pha_on.fits'.format(obs.id())
+ offname = prefix + '_{}_pha_off.fits'.format(obs.id())
+ arfname = prefix + '_{}_arf.fits'.format(obs.id())
+ rmfname = prefix + '_{}_rmf.fits'.format(obs.id())
+ else:
+ onname = prefix + '_pha_on.fits'
+ offname = prefix + '_pha_off.fits'
+ arfname = prefix + '_arf.fits'
+ rmfname = prefix + '_rmf.fits'
+
+ # Save files
+ obs.on_spec().save(onname, clobber)
+ obs.off_spec().save(offname, clobber)
+ obs.arf().save(arfname, clobber)
+ obs.rmf().save(rmfname, clobber)
+
+ # Log file names
+ self._log_value(gammalib.NORMAL, 'PHA on file', onname)
+ self._log_value(gammalib.NORMAL, 'PHA off file', offname)
+ self._log_value(gammalib.NORMAL, 'ARF file', arfname)
+ self._log_value(gammalib.NORMAL, 'RMF file', arfname)
+
+ # Save observation definition XML file
+ self._outobs.save(outobs)
+
+ # Log file name
+ self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs)
+
+ # Save ds9 On region file
+ regname = prefix + '_on.reg'
+ self._src_reg.save(regname)
+ self._log_value(gammalib.NORMAL, 'On region file', regname)
+
+ # Save ds9 Off region files
+ for bkg_reg in self._bkg_regs:
+
+ # Set filename
+ if len(self._bkg_regs) > 1:
+ regname = prefix + '_{}_off.reg'.format(bkg_reg['id'])
+ else:
+ regname = prefix + '_off.reg'
+
+ # Save ds9 region file
+ bkg_reg['regions'].save(regname)
+
+ # Log file name
+ self._log_value(gammalib.NORMAL, 'Off region file', regname)
+
+ # Return
+ return
+
+ def obs(self):
+ """
+ Return observation container
+ """
+ return self._outobs
+
+ def execute(self):
+ """
+ Execute the script
+ """
+ # Open logfile
+ self.logFileOpen()
+
+ # Read ahead output parameters
+ self._read_ahead(True)
+
+ # Run the script
+ self.run()
+
+ # Save result
+ self.save()
+
+ # Return
+ return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Create instance of application
+ app = csphagen(sys.argv)
+
+ # Execute application
+ app.execute()
diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst
index 1debf07..749cfe3 100644
--- a/doc/source/admin/release_history/1.5.rst
+++ b/doc/source/admin/release_history/1.5.rst
@@ -10,6 +10,7 @@ ctools 1.5 is a major release that adds significant functionality.
In particular, this release provides:
+* Support of classical IACT analysis
* Specifcation of times in various formats (UTC, MJD, JD, MET)
@@ -31,6 +32,8 @@ Improvements
New features
------------
+* [`2208 `_] -
+ Add :ref:`csphagen` script
* [`2180 `_] -
Add :ref:`csmodelsois` script
* [`1864 `_] -
diff --git a/doc/source/users/reference_manual/csphagen.rst b/doc/source/users/reference_manual/csphagen.rst
new file mode 100644
index 0000000..9e9fb94
--- /dev/null
+++ b/doc/source/users/reference_manual/csphagen.rst
@@ -0,0 +1,137 @@
+.. _csphagen:
+
+csphagen
+========
+
+Creates files necessary to perform a classical region-based spectral analysis
+on IACT data, i.e., source and background PHA files, ARF and RMF files.
+
+
+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).
+`See here `__.
+
+The outputs are:
+1) the PHA, ARF, RMF files, either separately for each observation, or stacked
+ for a set of observations,
+2) DS9 regions files listing the source and background regions for each
+ observation;
+3) a new observation definition XML file.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+ Input event list or observation definition XML file.
+
+``caldb [string]``
+ Calibration database.
+
+``irf [string]``
+ 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).
+
+``outobs [string]``
+ Output observation definition XML file.
+
+``(prefix = onoff) [string]``
+ 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.
+
+``emax [real]``
+ 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.
+
+``ebinfile [file]``
+ Name of the file containing the energy bin definition if FILE algorithm
+ is used.
+
+``(srcshape = CIRCLE) [string]``
+ Shape of the source region.
+ CIRCLE: circular region around given position.
+
+``coordsys [string]``
+ Coordinate system (CEL - celestial, GAL - galactic)
+
+``ra [real]``
+ Right Ascension of source region centre (deg)
+
+``dec [real]``
+ Declination of source region centre (deg)
+
+``glon [real]``
+ Galactic longitude of source region centre (deg)
+
+``glat [real]``
+ Galactic latitude of source region centre (deg)
+
+``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.
+
+``(bkgregmin = 2) [integer]``
+ Minimum number of background regions that are required for an observation.
+ If this number of background regions is not available the observation is
+ skipped.
+
+``(maxoffset = 4.0) [real]``
+ Maximum offset in degrees of source from camera center to accept the
+ observation.
+
+``stack [boolean]``
+ Specifies whether multiple observations should be stacked (``yes``) or
+ whether run-wise PHA, ARF and RMF files should be produced (``no``).
+
+
+Standard parameters
+-------------------
+
+``(chatter = 2) [integer]``
+ Verbosity of the executable:
+ ``chatter = 0``: no information will be logged
+
+ ``chatter = 1``: only errors will be logged
+
+ ``chatter = 2``: errors and actions will be logged
+
+ ``chatter = 3``: report about the task execution
+
+ ``chatter = 4``: detailed report about the task execution
+
+``(clobber = yes) [boolean]``
+ Specifies whether an existing output runlist should be overwritten.
+
+``(debug = no) [boolean]``
+ Enables debug mode. In debug mode the executable will dump any log file
+ output to the console.
+
+``(mode = ql) [string]``
+ Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = csphagen.log) [filename]``
+ Log filename.
+
+
diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst
index 9e6c17b..59cbbdc 100644
--- a/doc/source/users/reference_manual/index.rst
+++ b/doc/source/users/reference_manual/index.rst
@@ -50,6 +50,7 @@ cscripts
csobsdef --- Generates observation definition file
csobsinfo --- Shows observation container content
csobsselect --- Select observations from observation definition file
+ csphagen —-- Generates PHA, ARF, RMF files based on source/background regions
csphasecrv --- Computes phase curve
cspull --- Generates pull distribution
csresmap --- Generates residual map
diff --git a/test/data/crab_exclusion.fits b/test/data/crab_exclusion.fits
new file mode 100644
index 0000000..c6be7da
--- /dev/null
+++ b/test/data/crab_exclusion.fits
diff --git a/test/data/crab_offaxis1.fits b/test/data/crab_offaxis1.fits
new file mode 100644
index 0000000..584bab4
--- /dev/null
+++ b/test/data/crab_offaxis1.fits
@@ -0,0 +1,553 @@
+SIMPLE = T / file does conform to FITS standard BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 36 / width of table in bytes NAXIS2 = 3254 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 8 / number of fields in each row TTYPE1 = 'EVENT_ID' / label for field 1 TFORM1 = '1J ' / data format of field: 4-byte INTEGER TZERO1 = 2147483648 / offset for unsigned integers TSCAL1 = 1 / data are not scaled TTYPE2 = 'TIME ' / label for field 2 TFORM2 = '1D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field TTYPE3 = 'RA ' / label for field 3 TFORM3 = '1E ' / data format of field: 4-byte REAL TUNIT3 = 'deg ' / physical unit of field TTYPE4 = 'DEC ' / label for field 4 TFORM4 = '1E ' / data format of field: 4-byte REAL TUNIT4 = 'deg ' / physical unit of field TTYPE5 = 'ENERGY ' / label for field 5 TFORM5 = '1E ' / data format of field: 4-byte REAL TUNIT5 = 'TeV ' / physical unit of field TTYPE6 = 'DETX ' / label for field 6 TFORM6 = '1E ' / data format of field: 4-byte REAL TUNIT6 = 'deg ' / physical unit of field TTYPE7 = 'DETY ' / label for field 7 TFORM7 = '1E ' / data format of field: 4-byte REAL TUNIT7 = 'deg ' / physical unit of field TTYPE8 = 'MC_ID ' / label for field 8 TFORM8 = '1J ' / data format of field: 4-byte INTEGER EXTNAME = 'EVENTS ' / name of this extension DSTYP1 = 'TIME ' / Data sub-space type DSUNI1 = 's ' / Data sub-space unit DSVAL1 = 'TABLE ' / Data sub-space value DSREF1 = ':GTI ' / Data sub-space reference DSTYP2 = 'ENERGY ' / Data sub-space type DSUNI2 = 'TeV ' / Data sub-space unit DSVAL2 = '0.1:100 ' / Data sub-space value DSTYP3 = 'POS(RA,DEC)' / Data sub-space type DSUNI3 = 'deg ' / Data sub-space unit DSVAL3 = 'CIRCLE(83,22.5,5)' / Data sub-space value NDSKEYS = 3 / Number of data sub-space keys NMCIDS = 2 / Number of Monte Carlo identifiers MID00001= 1 / Monte Carlo identifier for model 1 MMN00001= 'Crab ' / Name of model 1 MID00002= 2 / Monte Carlo identifier for model 2 MMN00002= 'Background' / Name of model 2 CREATOR = 'GammaLib' / Program which created the file TELESCOP= 'CTA ' / Telescope OBS_ID = 0 / Observation identifier DATE_OBS= '2020-01-01' / Observation start date TIME_OBS= '00:00:00' / Observation start time DATE_END= '2020-01-01' / Observation end date TIME_END= '00:15:00' / Observation end time TSTART = 6.3110886918E+08 / [s] Mission time of start of observation TSTOP = 6.3110976918E+08 / [s] Mission time of end of observation MJDREFI = 51544 / [days] Integer part of time reference MJD MJDREFF = 5.0000000000E-01 / [days] Fractional part of time reference MJD TIMEUNIT= 's ' / Time unit TIMESYS = 'TT ' / Time system TIMEREF = 'LOCAL ' / Time reference TELAPSE = 8.9999999976E+02 / [s] Mission elapsed time ONTIME = 8.9999999976E+02 / [s] Total good time including deadtime LIVETIME= 8.8199999977E+02 / [s] Total livetime DEADC = 9.8000000000E-01 / Deadtime correction factor TIMEDEL = 1.0000000000E+00 / Time resolution OBJECT = ' ' / Observed object RA_OBJ = 0.0000000000E+00 / [deg] Target Right Ascension DEC_OBJ = 0.0000000000E+00 / [deg] Target Declination RA_PNT = 8.3000000000E+01 / [deg] Pointing Right Ascension DEC_PNT = 2.2500000000E+01 / [deg] Pointing Declination ALT_PNT = 9.0000000000E+01 / [deg] Average altitude of pointing AZ_PNT = 0.0000000000E+00 / [deg] Average azimuth of pointing RADECSYS= 'FK5 ' / Coordinate system EQUINOX = 2.0000000000E+03 / Epoch CONV_DEP= 0.0000000000E+00 / Convergence depth of telescopes CONV_RA = 0.0000000000E+00 / [deg] Convergence Right Ascension CONV_DEC= 0.0000000000E+00 / [deg] Convergence Declination OBSERVER= 'string ' / Observer N_TELS = 0 / Number of telescopes in event list TELLIST = 'string ' / Telescope IDs GEOLAT = 0.0000000000E+00 / [deg] Geographic latitude of array centre GEOLON = 0.0000000000E+00 / [deg] Geographic longitude of array centre ALTITUDE= 0.0000000000E+00 / [km] Altitude of array centre EUNIT = 'TeV ' / Energy unit EVTVER = 'draft1 ' / Event list version number END € AÂÎú‚¢!¹B§X‹A¯T>EYò¿!¥?ÒØ € AÂÎúƒˆÀ%B§p\A¯–•>:
+¿Å?*ÔÜ € AÂÎúƒß¾`B§?²A°
+G=λ ¾üÑ?4> € AÂÎú…€ªB§5¶A¯Â>r¢¿½þ?¥2 € AÂÎú†FPB§lA°yÐ=ÞQ|¾á&)> € AÂÎú†3dB§;A¯ft>±w¿å%ª € AÂÎú†…
›B§cA¯³š>P\¿ S½>üžê € AÂÎú‡R„B§IA°eÜ=èž’¾åâc?~ˆ € AÂÎú‡�WB§P5A¯ÕÂ>;O¿ðÒ?ã €
+AÂÎú‡âNB§�A±A =ïf!¾¯>?bL € AÂÎú‰÷ÞFB§.ÙA¯)?>9Ø¿‘Ÿ?� € AÂÎúŠí?B§WlA°X'>¾é@÷?) €
AÂÎú‹ªÅB§îA¯¼l>0C¿7•? ií € AÂÎú�‚^hB§"›A°aþ>E…Õ¾æþ°?T € AÂÎúŽ˜yB§nA¯æL>Dlã¿Ï“?)¶M € AÂÎúŽÔ?*B§X=A¯b³=æ]Ê¿N’?¬¸ € AÂÎú’d’B§BbA°¶ >4éÂ¾ÑØ*?\f € AÂÎú’œ‡HB§>QA¯x>¬¿±W?£ß € AÂÎú“‹SBB§ÖA°ùf=äE”¾Á'þ?Ü› € AÂÎú“É¥;B§MYA¯c,=Ù¥ä¿E)?ŸŽ € AÂÎú“âAÀB§zKA°¯Ò=êÒ=¾Ó/j?/D� € AÂÎú”yÛÖB§OÉA°¢'>$ö�¾ÖÈä?”¬ € AÂÎú”ᾄB§-jA°†š>:…6¾ÝͰ?«F € AÂÎú•ö8B§NPA¯¸Î>(™¿�`? € AÂÎú•E°SB§@(A±Î!=ö ¥¾‹ÙÏ?/+ € AÂÎú•uÇXB§�§A¯NX=Û2޿¼?2å¨ € AÂÎú–±înB¦ÿA±!—=òŠL¾·6>ì<0 € AÂÎú—ª!B§LÄA±4þ>@æ ¾²?g € AÂÎú˜OòB§VžA¯ˆ²=Ö§6¿�‡?ç € AÂÎú˜¤üRB§A°6[>³í¾ñð|?ÖN € AÂÎú˜ÖŠäB§>œA°ô‘>ÁU¾Â?¥?”z € AÂÎú™.B§1*A°Øa=ÏžÉXš?
^ € !AÂÎú™ÐÀB§.!A¯Óh=ß×½¿LÓ?4 € "AÂÎú™Þu*B¦õÓA°a^=Í<¾çK~>ãÙ} € #AÂÎúšyNpB§UˆA¯ÖN>0Ãï¿Ü‘?[<