diff --git a/ChangeLog b/ChangeLog
index d9a724e..d11f822 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,9 @@
-2017-03-21
+2017-03-22
* Version 1.3.0 released
========================
+ Add csebins script (#1976)
Take into account RoI in ctskymap and use numerical integration for IRF (#1975)
Set local language to english (#1370)
Enhance logfile output of csroot2caldb script
diff --git a/NEWS b/NEWS
index 2148b77..7068abb 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
New Features and Important Changes in ctools 1.3.0.dev1
-21 March 2017
+22 March 2017
Introduction
@@ -112,7 +112,10 @@ None
cscripts
--------
-None
+The csebins script has been added that generates an energy boundaries
+extension that can be used for a stacked analysis. The energy boundaries
+are set as function of the variation in the effective area and background
+rates of the Instrument Response Function (#1976).
Examples
diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index d6a4ad1..4134487 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -22,6 +22,7 @@
# Scripts to be distributed in $(pyexecdir)/cscripts
pythonscripts = $(top_builddir)/cscripts/__init__.py \
$(srcdir)/cscaldb.py \
+ $(srcdir)/csebins.py \
$(srcdir)/csinfo.py \
$(srcdir)/cslightcrv.py \
$(srcdir)/csmodelinfo.py \
@@ -53,6 +54,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
# Parfiles to be distributed in $(prefix)/syspfiles
parfiles = $(srcdir)/cscaldb.par \
+ $(srcdir)/csebins.par \
$(srcdir)/cslightcrv.par \
$(srcdir)/csmodelinfo.par \
$(srcdir)/csmodelmerge.par \
@@ -79,6 +81,7 @@ parfiles = $(srcdir)/cscaldb.par \
# Scripts for which symbolic links will be created in bindir
linkscripts = cscaldb \
+ csebins \
csinfo \
cslightcrv \
csmodelinfo \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index 34c92d1..21554fa 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -25,6 +25,7 @@ cur_version = sys.version_info
# List all scripts in this module
__all__ = [
"cscaldb",
+ "csebins",
"csinfo",
"cslightcrv",
"csmodelinfo",
@@ -59,6 +60,7 @@ if cur_version > req_version:
# Import
from cscripts.cscaldb import cscaldb
from cscripts.csinfo import csinfo
+from cscripts.csebins import csebins
from cscripts.cslightcrv import cslightcrv
from cscripts.csmodelinfo import csmodelinfo
from cscripts.csmodelmerge import csmodelmerge
diff --git a/cscripts/csebins.par b/cscripts/csebins.par
new file mode 100644
index 0000000..bf11659
--- /dev/null
+++ b/cscripts/csebins.par
@@ -0,0 +1,39 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: csebins.par #
+# #
+# Author: Juergen Knoedlseder #
+# Jurgen.Knodlseder@irap.omp.eu #
+# IRAP #
+# #
+# Purpose: Parameter file for the csebins script #
+# #
+#############################################################################
+#
+# Input and output parameters
+#============================
+inobs, f, a, NONE,,, "Input observation definition XML file"
+caldb, s, a, prod2,,, "Calibration database"
+irf, s, a, South_0.5h,,, "Instrument response function"
+outfile, f, a, ebounds.fits,,, "Output energy boundaries file"
+
+#
+# Script parameters
+#==================
+emin, r, a, 0.02,,, "Minimum energy (TeV)"
+emax, r, a, 300.0,,, "Maximum energy (TeV)"
+aeffthres, r, h, 0.2,,, "Fractional effective area change threshold"
+bkgthres, r, h, 0.5,,, "Fractional background rate change threshold"
+
+#
+# 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, csebins.log,,, "Log filename"
diff --git a/cscripts/csebins.py b/cscripts/csebins.py
new file mode 100755
index 0000000..9375529
--- /dev/null
+++ b/cscripts/csebins.py
@@ -0,0 +1,335 @@
+#! /usr/bin/env python
+# ==========================================================================
+# Generates energy boundaries for stacked analysis
+#
+# Copyright (C) 2017 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
+# 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 sys
+import math
+import gammalib
+import ctools
+#from cscripts import obsutils
+
+
+# ============= #
+# csebins class #
+# ============= #
+class csebins(ctools.cscript):
+ """
+ Generates energy boundaries for stacked analysis
+ """
+
+ # Constructor
+ def __init__(self, *argv):
+ """
+ Constructor
+ """
+ # Set name and version
+ self._name = 'csebins'
+ self._version = '1.0.0'
+
+ # Set members
+ self._ebounds = gammalib.GEbounds()
+
+ # Initialise observation container from constructor arguments.
+ self._obs, argv = self._set_input_obs(argv)
+
+ # Initialise script by calling the appropriate class constructor.
+ self._init_cscript(argv)
+
+ # Return
+ return
+
+
+ # Private methods
+ def _get_parameters(self):
+ """
+ Get parameters from parfile
+ """
+ # Initialise observation container if it is currently empty
+ if self._obs.size() == 0:
+
+ # Query "inobs" parameter
+ filename = self['inobs'].filename()
+
+ # If an observation definition file was provided then load it,
+ # otherwise build a single observation with response information
+ if self._is_valid_filename(filename):
+ self._obs = gammalib.GObservations(filename)
+ else:
+ cta = gammalib.GCTAObservation()
+ caldb = gammalib.GCaldb('cta', self['caldb'].string())
+ rsp = gammalib.GCTAResponseIrf(self['irf'].string(), caldb)
+ cta.response(rsp)
+ self._obs.append(cta)
+
+ # Query input parameters
+ self['emin'].real()
+ self['emax'].real()
+ self['aeffthres'].real()
+ self['bkgthres'].real()
+
+ # Query ahead output model filename
+ if self._read_ahead():
+ self['outfile'].filename()
+
+ # Write input parameters into logger
+ self._log_parameters(gammalib.TERSE)
+
+ # Return
+ return
+
+ def _insert_energy(self, responses, energy, comment):
+ """
+ Set energy boundaries for one CTA observation
+
+ Parameters
+ ----------
+ responses : list of dict
+ List of response dictionaries
+ energy : float
+ Energy to insert (TeV)
+ comment : str
+ Reason for energy insertion
+ """
+ # Set thresholds of all response components
+ logE = math.log10(energy)
+ for rsp in responses:
+ rsp['thres'] = rsp['irf'](logE, 0.0, 0.0)
+
+ # Convert energy into gammalib energy
+ eng = gammalib.GEnergy(energy, 'TeV')
+
+ # Set insertion message
+ msg = str(energy)+' TeV'
+ if len(comment) > 0:
+ msg += ' ('+comment+')'
+
+ # If energy boundaries are empty then append an energy interval
+ # with zero size ...
+ if self._ebounds.size() == 0:
+
+ # Append energy interval with zero size
+ self._ebounds.append(eng, eng)
+
+ # Log first energy boundary
+ self._log_value(gammalib.NORMAL, 'First boundary', msg)
+
+ # ... otherwise append an interval
+ elif self._ebounds.size() > 0:
+
+ # If current boundaries have zero size then recover boundary energy,
+ # remove boundary, and append a boundary of non-zero size
+ if self._ebounds.emin() == self._ebounds.emax():
+
+ # Recover old energy, remove boundary and append new boundary
+ emin = self._ebounds.emin().copy() # Because emin becomes zero after remove
+ emax = eng
+ self._ebounds.remove(0)
+ self._ebounds.append(emin, emax)
+
+ # Log other energy boundary
+ self._log_value(gammalib.NORMAL, 'Boundary', msg)
+
+ else:
+
+ # Use old emax as new emin. Only append if interval is non zero
+ emin = self._ebounds.emax()
+ if emin < eng:
+ self._ebounds.append(emin, eng)
+
+ # Log other energy boundary
+ self._log_value(gammalib.NORMAL, 'Boundary', msg)
+
+ # Return
+ return
+
+ def _set_ebounds(self, obs):
+ """
+ Set energy boundaries for one CTA observation
+
+ Parameters
+ ----------
+ obs : `~gammalib.GObservations`
+ CTA observation
+ """
+ # Get parameters
+ emin = self['emin'].real()
+ emax = self['emax'].real()
+ aeffthres = self['aeffthres'].real()
+ bkgthres = self['bkgthres'].real()
+
+ # Build list of response dictionaries
+ responses = []
+
+ # Loop over all observations
+ for run in obs:
+
+ # Get response name
+ mission = run.response().caldb().mission()
+ instrument = run.response().caldb().instrument()
+ rspname = run.response().rspname()
+ name = '%s::%s::%s' % (mission, instrument, rspname)
+
+ # If response name exists already then examine next observation
+ exists = False
+ for rsp in responses:
+ if rsp['name'] == name:
+ exists = True
+ break
+ if exists:
+ continue
+
+ # Append effective area and background template to list
+ aeff = {'name': name, 'type': 'Aeff',
+ 'irf': run.response().aeff(), 'thres': 0.0}
+ bkg = {'name': name, 'type': 'Background',
+ 'irf': run.response().background(), 'thres': 0.0}
+ responses.append(aeff)
+ responses.append(bkg)
+
+ # Log response name
+ self._log_value(gammalib.NORMAL, 'Append response', '%s (%s) [%s]' %
+ (rspname, instrument, mission))
+
+ # Setup energy vector in log10(TeV)
+ logEs = [math.log10(1.0e-3*float(i)) for i in range(int(emin*1000),int(emax*1000))]
+
+ # Loop over all energies
+ for logE in logEs:
+
+ # Loop over all response components
+ for rsp in responses:
+
+ # Get IRF value
+ irf = rsp['irf'](logE, 0.0, 0.0)
+
+ # Get threshold
+ if rsp['type'] == 'Aeff':
+ threshold = aeffthres
+ else:
+ threshold = bkgthres
+
+ # If threshold is zero and response is positive then insert
+ # an energy ...
+ if rsp['thres'] == 0.0 and irf > 0.0:
+ self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
+ break
+
+ # ... otherwise if threshold is not zero then check whether
+ # the fractional change exceeds the threshold
+ elif rsp['thres'] != 0.0:
+ f = irf/rsp['thres']
+ if 1.0-f > threshold or f-1.0 > threshold:
+ self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
+ break
+
+ # Insert energy at end of range
+ self._insert_energy(responses, emax, 'End of range')
+
+ # Log number of energy boundaries
+ self._log_value(gammalib.NORMAL, 'Number of boundaries', self._ebounds.size())
+
+ # Return
+ return
+
+
+ # 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, 'Input observation')
+
+ # Write header
+ self._log_header1(gammalib.TERSE, 'Define energy boundaries')
+
+ # Define energy boundaries for observation container
+ self._set_ebounds(self._obs)
+
+ # Return
+ return
+
+ def execute(self):
+ """
+ Execute the script
+ """
+ # Open logfile
+ self.logFileOpen()
+
+ # Read ahead output parameters
+ self._read_ahead(True)
+
+ # Run the scrip
+ self.run()
+
+ # Save energy boundaries
+ self.save()
+
+ # Return
+ return
+
+ def save(self):
+ """
+ Save energy boundaries
+ """
+ # Write header
+ self._log_header1(gammalib.TERSE, 'Save energy boundaries')
+
+ # Get energy boundaries filename
+ outfile = self['outfile'].filename()
+
+ # Log file name
+ self._log_value(gammalib.NORMAL, 'Energy boundaries file', outfile.url())
+
+ # Save energy boundaries
+ self._ebounds.save(outfile, self._clobber())
+
+ # Return
+ return
+
+ def ebounds(self):
+ """
+ Return energy boundaries
+
+ Returns
+ -------
+ ebounds : `~gammalib.GEbounds()`
+ Energy boundaries
+ """
+ # Return
+ return self._ebounds
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Create instance of application
+ app = csebins(sys.argv)
+
+ # Execute application
+ app.execute()
diff --git a/doc/source/users/reference_manual/csebins.rst b/doc/source/users/reference_manual/csebins.rst
new file mode 100644
index 0000000..7a2ca11
--- /dev/null
+++ b/doc/source/users/reference_manual/csebins.rst
@@ -0,0 +1,85 @@
+.. _csebins:
+
+csebins
+=======
+
+Generates energy boundaries for stacked analysis.
+
+
+Synopsis
+--------
+
+This script generates the energy bounaries for a stacked analysis by inspecting
+the variation in the effective area and the background template of the response
+functions. The script takes either an observation definition XML file or a
+specific instrument response function on input.
+
+csebins writes an energy boundary extension into a FITS file.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+ Input observation definition XML file.
+
+``caldb [string]``
+ Calibration database.
+
+``irf [string]``
+ Instrumental response function.
+
+``outfile [file]``
+ Name of the energy boundary output file.
+
+``emin [real]``
+ Lower energy limit of energy boundaries (in TeV).
+
+``emax [real]``
+ Upper energy limit of energy boundaries (in TeV).
+
+``(aeffthres = 0.2) [real]``
+ Fractional change in effective area that leads to insertion of a new energy
+ boundary.
+
+``(bkgthres = 0.5) [real]``
+ Fractional change in background rate that leads to insertion of a new energy
+ boundary.
+
+
+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 energy boundaries output file 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 = csebins.log) [filename]``
+ Log filename.
+
+
+Related tools or scripts
+------------------------
+
+:doc:`ctbin`
+:doc:`ctexpcube`
+:doc:`ctpsfcube`
+:doc:`ctedispcube`
+:doc:`ctbkgcube`
diff --git a/doc/source/users/reference_manual/cslightcrv.rst b/doc/source/users/reference_manual/cslightcrv.rst
index 6ed5d78..d4c7509 100644
--- a/doc/source/users/reference_manual/cslightcrv.rst
+++ b/doc/source/users/reference_manual/cslightcrv.rst
@@ -46,7 +46,7 @@ General parameters
Applies energy dispersion to response computation.
``outfile [file]``
- Name of the source spectrum output file.
+ Name of the light curve output file.
``(tbinalg = GTI) [string]``
Algorithm for defining time bins.
@@ -125,7 +125,7 @@ Standard parameters
chatter = 4: detailed report about the task execution
``(clobber = yes) [boolean]``
- Specifies whether an existing source spectrum output file should be overwritten.
+ Specifies whether an existing light curve output file should be overwritten.
``(debug = no) [boolean]``
Enables debug mode. In debug mode the executable will dump any log file output to the console.
diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst
index 94b179e..332e84b 100644
--- a/doc/source/users/reference_manual/index.rst
+++ b/doc/source/users/reference_manual/index.rst
@@ -40,6 +40,7 @@ cscripts
:maxdepth: 1
cscaldb --- Lists available instrument response functions
+ csebins --- Generates energy boundaries for stacked analysis
cslightcrv --- Computes lightcurve
csmodelinfo --- Shows model container content
csmodelmerge --- Merges several model containers into one file
diff --git a/test/test_csebins.py b/test/test_csebins.py
new file mode 100755
index 0000000..743c10d
--- /dev/null
+++ b/test/test_csebins.py
@@ -0,0 +1,142 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cslightcrv script.
+#
+# Copyright (C) 2016-2017 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
+# 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 cscripts
+from testing import test
+
+
+# ============================= #
+# Test class for csebins script #
+# ============================= #
+class Test(test):
+ """
+ Test class for csebins script
+
+ This test class makes unit tests for the csebins script by using it from
+ the command line and from Python.
+ """
+
+ # Constructor
+ def __init__(self):
+ """
+ Constructor
+ """
+ # Call base class constructor
+ test.__init__(self)
+
+ # Return
+ return
+
+ # Set test functions
+ def set(self):
+ """
+ Set all test functions.
+ """
+ # Set test name
+ self.name('csebins')
+
+ # Append tests
+ self.append(self._test_cmd, 'Test csebins on command line')
+ self.append(self._test_python, 'Test csebins from Python')
+
+ # Return
+ return
+
+ # Test csebins on command line
+ def _test_cmd(self):
+ """
+ Test csebins on the command line.
+ """
+ # Set script name
+ cslightcrv = self._script('csebins')
+
+ # Setup cslightcrv command
+ cmd = cslightcrv+' inobs="NONE"'+ \
+ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+ ' emin=0.1 emax=100.0'+ \
+ ' outfile="csebins_cmd1.fits"'+ \
+ ' logfile="csebins_cmd1.log" chatter=1'
+
+ # Check if execution of wrong command fails
+ self.test_assert(self._execute('command_that_does_not_exist') != 0,
+ 'Self test of test script')
+
+ # Check if execution was successful
+ self.test_assert(self._execute(cmd) == 0,
+ 'Check successful execution from command line')
+
+ # Check energy boundaries
+ self._check_ebounds('csebins_cmd1.fits', 33)
+
+ # Setup cslightcrv command
+ cmd = cslightcrv+' inobs="NONE"'+ \
+ ' caldb="Database_that_does_not_exist" irf="'+self._irf+'"'+ \
+ ' emin=0.1 emax=100.0'+ \
+ ' outfile="csebins_cmd2.fits"'+ \
+ ' logfile="csebins_cmd2.log" chatter=1'
+
+ # Check if execution failed
+ self.test_assert(self._execute(cmd) != 0,
+ 'Check invalid input file when executed from command line')
+
+ # Return
+ return
+
+ # Test csebins from Python
+ def _test_python(self):
+ """
+ Test csebins from Python
+ """
+ # Set-up csebins
+ ebins = cscripts.csebins()
+ ebins['inobs'] = 'NONE'
+ ebins['caldb'] = self._caldb
+ ebins['irf'] = self._irf
+ ebins['emin'] = 0.1
+ ebins['emax'] = 100.0
+ ebins['outfile'] = 'csebins_py1.fits'
+ ebins['logfile'] = 'csebins_py1.log'
+ ebins['chatter'] = 2
+
+ # Run csebins script and energy boundaries
+ ebins.logFileOpen() # Make sure we get a log file
+ ebins.run()
+ ebins.save()
+
+ # Check energy boundaries
+ self._check_ebounds('csebins_py1.fits', 33)
+
+ # Return
+ return
+
+ # Check energy boundaries result file
+ def _check_ebounds(self, filename, bins):
+ """
+ Check energy boundaries file
+ """
+ # Load energy boundaries file
+ ebounds = gammalib.GEbounds(filename)
+
+ # Check number of energy boundaries
+ self.test_value(ebounds.size(), bins)
+
+ # Return
+ return
diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py
index 177dd47..944723a 100755
--- a/test/test_python_cscripts.py
+++ b/test/test_python_cscripts.py
@@ -2,7 +2,7 @@
# ==========================================================================
# This scripts performs unit tests for cscripts
#
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 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
@@ -25,6 +25,7 @@ import cscripts
sys.path.append(os.path.dirname(os.path.abspath(__file__)))
import test_cscript
import test_cscaldb
+import test_csebins
import test_csfindobs
import test_csinfo
import test_cslightcrv
@@ -113,6 +114,7 @@ def test(installed=False, debug=False):
# Define list of test suites
tests = [test_cscript.Test(),
test_cscaldb.Test(),
+ test_csebins.Test(),
test_csinfo.Test(),
test_cslightcrv.Test(),
test_csmodelinfo.Test(),