diff --git a/ChangeLog b/ChangeLog
index baa09b2..1fa0f60 100755
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,9 @@
-2018-10-17
+2018-10-18
* Version 1.6.0 released
========================
+ Add csbkgmodel script (#2694)
Compute DETX and DETY in ctmodel if they are not available (#2694)
Do no longer set missing boundaries in ctool::get_observations and
ctool::setup_observations (#2686)
diff --git a/NEWS b/NEWS
index 51168a4..54343c5 100755
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,6 @@
New Features and Important Changes in ctools 1.6.0
-17 October 2018
+18 October 2018
Introduction
@@ -200,6 +200,13 @@ support - Support classes
None
+csbkgmodel - Generate background model for 3D analysis
+------------------------------------------------------
+Introduced the script that generates a background for a 3D analysis. This
+script is particularly useful to generate a background model for the analysis
+of the H.E.S.S. public data release (#2694).
+
+
cscaldb - Inspect calibration database
--------------------------------------
None
diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index 8c33c20..4a67755 100755
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -21,6 +21,7 @@
# Scripts to be distributed in $(pyexecdir)/cscripts
pythonscripts = $(top_builddir)/cscripts/__init__.py \
+ $(srcdir)/csbkgmodel.py \
$(srcdir)/cscaldb.py \
$(srcdir)/csebins.py \
$(srcdir)/csinfo.py \
@@ -58,7 +59,8 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
$(srcdir)/mputils.py
# Parfiles to be distributed in $(prefix)/syspfiles
-parfiles = $(srcdir)/cscaldb.par \
+parfiles = $(srcdir)/csbkgmodel.par \
+ $(srcdir)/cscaldb.par \
$(srcdir)/csebins.par \
$(srcdir)/cslightcrv.par \
$(srcdir)/csmodelinfo.par \
@@ -89,7 +91,8 @@ parfiles = $(srcdir)/cscaldb.par \
$(srcdir)/csiactobs.par
# Scripts for which symbolic links will be created in bindir
-linkscripts = cscaldb \
+linkscripts = csbkgmodel \
+ cscaldb \
csebins \
csinfo \
cslightcrv \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index 91aaa4c..5327f27 100755
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -24,6 +24,7 @@ cur_version = sys.version_info
# List all scripts in this module
__all__ = [
+ "csbkgmodel",
"cscaldb",
"csebins",
"csinfo",
@@ -63,9 +64,10 @@ if cur_version > req_version:
"csiactobs"])
# Import
+from cscripts.csbkgmodel import csbkgmodel
from cscripts.cscaldb import cscaldb
-from cscripts.csinfo import csinfo
from cscripts.csebins import csebins
+from cscripts.csinfo import csinfo
from cscripts.cslightcrv import cslightcrv
from cscripts.csmodelinfo import csmodelinfo
from cscripts.csmodelmerge import csmodelmerge
diff --git a/cscripts/csbkgmodel.par b/cscripts/csbkgmodel.par
new file mode 100755
index 0000000..5b2b91c
--- /dev/null
+++ b/cscripts/csbkgmodel.par
@@ -0,0 +1,39 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: csbkgmodel.par #
+# #
+# Author: Juergen Knoedlseder #
+# Jurgen.Knodlseder@irap.omp.eu #
+# IRAP #
+# #
+# Purpose: Parameter file for the csbkgmodel script #
+# #
+#############################################################################
+#
+# Input and output parameters
+#============================
+inobs, f, a, obs.xml,,, "Input observation definition XML file"
+outmodel, f, a, model.xml,,, "Output model definition XML file"
+
+#
+# Script parameters
+#==================
+instrument, s, a, CTA,,, "Instrument name (for multi-instrument analysis)"
+spatial, s, a, GAUSS,IRF|AEFF|GAUSS,, "Spatial model component"
+gradient, b, a, yes,,, "Allow for spatial gradient?"
+spectral, s, a, HESS,PLAW|HESS,, "Spectral model component"
+runwise, b, a, yes,,, "Generate runwise background model?"
+rad, r, h, 2.0,,, "Radius for event selection (degrees)"
+
+#
+# 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,,, "Activate debugging mode?"
+mode, s, h, ql,,, "Mode of automatic parameters"
+logfile, f, h, csbkgmodel.log,,, "Log filename"
diff --git a/cscripts/csbkgmodel.py b/cscripts/csbkgmodel.py
new file mode 100755
index 0000000..5acbc00
--- /dev/null
+++ b/cscripts/csbkgmodel.py
@@ -0,0 +1,412 @@
+#! /usr/bin/env python
+# ==========================================================================
+# Generates background model
+#
+# Copyright (C) 2018 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
+
+
+# ================ #
+# csbkgmodel class #
+# ================ #
+class csbkgmodel(ctools.csobservation):
+ """
+ Generates background model
+ """
+
+ # Constructor
+ def __init__(self, *argv):
+ """
+ Constructor
+ """
+ # Initialise application by calling the appropriate class constructor
+ self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
+
+ # Set members
+ self._models = gammalib.GModels()
+ self._instrument = ''
+
+ # Return
+ return
+
+
+ # Private methods
+ def _get_parameters(self):
+ """
+ Get parameters from parfile
+ """
+ # Setup observations (do not require response, allow event list and
+ # counts cube)
+ self._setup_observations(self.obs(), False, True, True)
+
+ # Set instrument name
+ self._instrument = self._get_instrument()
+
+ # Query input parameters
+ spatial = self['spatial'].string()
+ if self['spatial'].string() == 'GAUSS':
+ self['gradient'].boolean()
+ self['spectral'].string()
+ self['runwise'].boolean()
+ self['rad'].real()
+
+ # Query ahead output model filename
+ if self._read_ahead():
+ self['outmodel'].filename()
+
+ # Write input parameters into logger
+ self._log_parameters(gammalib.TERSE)
+
+ # Return
+ return
+
+ def _get_instrument(self):
+ """
+ Get instrument name
+
+ Extracts the instrument name from the observation container. If there
+ are multiple instruments in the container then use the "instrument"
+ parameter to set the instrument name.
+
+ Returns
+ -------
+ instrument : str
+ Instrument name
+ """
+ # Initialise instrument
+ instrument = ''
+
+ # Loop over all observations
+ for obs in self.obs():
+ if instrument == '':
+ instrument = obs.instrument()
+ elif instrument != obs.instrument():
+ instrument = self['instrument'].string()
+ break
+
+ # Return instrument
+ return instrument
+
+ def _generate_initial_model(self):
+ """
+ Generate initial background model
+
+ Returns
+ -------
+ model : `~gammalib.GModelData()`
+ Background model
+ """
+ # Handle IRF model
+ if self['spatial'].string() == 'IRF':
+ epivot = gammalib.GEnergy(1.0, 'TeV')
+ spectral = gammalib.GModelSpectralPlaw(1.0, 0.0, epivot)
+ model = gammalib.GCTAModelIrfBackground(spectral)
+
+ # Handle AEFF model
+ elif self['spatial'].string() == 'AEFF':
+ epivot = gammalib.GEnergy(1.0, 'TeV')
+ spectral = gammalib.GModelSpectralPlaw(1.0e-13, -2.5, epivot)
+ model = gammalib.GCTAModelAeffBackground(spectral)
+
+ # Handle GAUSS model
+ elif self['spatial'].string() == 'GAUSS':
+ if self['gradient'].boolean():
+ spatial = gammalib.GCTAModelSpatialMultiplicative()
+ factor1 = gammalib.GCTAModelRadialGauss(3.0)
+ factor2 = gammalib.GCTAModelSpatialGradient()
+ spatial.append(factor1)
+ spatial.append(factor2)
+ else:
+ spatial = gammalib.GCTAModelRadialGauss(3.0)
+ epivot = gammalib.GEnergy(1.0, 'TeV')
+ spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
+ spectral['Prefactor'].min(1.0e-8)
+ model = gammalib.GCTAModelBackground(spatial, spectral)
+
+ # Any other strings (should never occur)
+ else:
+ model = None
+
+ # Set background name and instrument
+ if model is not None:
+ model.name('Background')
+ model.instruments(self._instrument)
+
+ # Return model
+ return model
+
+ def _select_events(self, obs):
+ """
+ Select events within a given RoI radius
+
+ Parameters
+ ----------
+ obs : `~gammalib.GObservations()`
+ Observation container
+
+ Returns
+ -------
+ obs : `~gammalib.GObservations()`
+ Observation container
+ """
+ # Setup task parameters
+ select = ctools.ctselect(obs)
+ select['ra'] = 'UNDEF'
+ select['dec'] = 'UNDEF'
+ select['rad'] = self['rad'].real()
+ select['tmin'] = 'UNDEF'
+ select['tmax'] = 'UNDEF'
+ select['emin'] = 'UNDEF'
+ select['emax'] = 'UNDEF'
+ select['usethres'] = 'DEFAULT'
+
+ # Select events
+ select.run()
+
+ # Extract observation
+ obs = select.obs().copy()
+
+ # Return observation
+ return obs
+
+ def _create_nodes(self, model):
+ """
+ Replace spectral model by node function
+
+ Parameters
+ ----------
+ model : `~gammalib.GModelData()`
+ Input background model
+
+ Returns
+ -------
+ model : `~gammalib.GModelData()`
+ Background model with spectral node function
+ """
+ # Extract spectral model of background component
+ spectrum = model.spectral()
+
+ # Define 8 node energies
+ energies = gammalib.GEnergies()
+ energies.append(gammalib.GEnergy(0.1, 'TeV'))
+ energies.append(gammalib.GEnergy(0.2, 'TeV'))
+ energies.append(gammalib.GEnergy(0.4, 'TeV'))
+ energies.append(gammalib.GEnergy(0.8, 'TeV'))
+ energies.append(gammalib.GEnergy(1.6, 'TeV'))
+ energies.append(gammalib.GEnergy(3.2, 'TeV'))
+ energies.append(gammalib.GEnergy(6.4, 'TeV'))
+ energies.append(gammalib.GEnergy(12.8, 'TeV'))
+
+ # Create node spectrum
+ nodes = gammalib.GModelSpectralNodes(spectrum, energies)
+ nodes.autoscale()
+
+ # Set minimum node intensities (this avoids NaNs during model fitting)
+ for i in range(nodes.size()):
+ nodes[i].min(1.0e-6*nodes[i].value())
+
+ # Set node spectrum
+ model.spectral(nodes)
+
+ # Return
+ return model
+
+ def _generate_bkg(self, obs):
+ """
+ Generate background models
+
+ Parameters
+ ----------
+ obs : `~gammalib.GObservations()`
+ Observations container
+
+ Returns
+ -------
+ model : `~gammalib.GModelData()`
+ Background model component
+ """
+ # Write header for event selection
+ self._log_header3(gammalib.EXPLICIT, 'Select events from observation')
+
+ # Select events
+ obs = self._select_events(obs)
+
+ # Write header for initial background model generation
+ self._log_header3(gammalib.EXPLICIT, 'Generate initial background model')
+
+ # Generate initial background model
+ model = self._generate_initial_model()
+
+ # Attach initial background model
+ models = gammalib.GModels()
+ models.append(model)
+ obs.models(models)
+
+ # Write header for initial model fitting
+ self._log_header3(gammalib.EXPLICIT, 'Fit initial background model')
+
+ # Perform maximum likelihood fitting with initial model
+ like = ctools.ctlike(obs)
+ like.run()
+
+ # Extract fitted model
+ model = like.obs().models()[0].copy()
+
+ # If a HESS model is requested then refit a node spectrum
+ if self['spectral'].string() == 'HESS':
+
+ # Create nodes spectrum from fitted initial model
+ model = self._create_nodes(model)
+
+ # Attach node spectrum
+ models = gammalib.GModels()
+ models.append(model)
+ obs.models(models)
+
+ # Write header for node model fitting
+ self._log_header3(gammalib.EXPLICIT, 'Fit nodes background model')
+
+ # Perform maximum likelihood fitting with node model
+ like = ctools.ctlike(obs)
+ like.run()
+
+ # Extract fitted model
+ model = like.obs().models()[0].copy()
+
+ # Return model
+ return model
+
+ def _generate_runwise_bkg(self):
+ """
+ Generate background models
+
+ Parameters
+ ----------
+ responses : list of dict
+ List of response dictionaries
+ energy : float
+ Energy to insert (TeV)
+ comment : str
+ Reason for energy insertion
+ """
+ # Loop over observations
+ for run in self.obs():
+
+ # Write header
+ self._log_header2(gammalib.TERSE, self._get_obs_header(run))
+
+ # Build observation container with single run
+ obs = gammalib.GObservations()
+ obs.append(run)
+
+ # Generate background model for
+ model = self._generate_bkg(obs)
+
+ # Set model attributes
+ model.name(('Background_%s' % run.id()))
+ model.ids(run.id())
+
+ # Write model
+ self._log_string(gammalib.NORMAL, str(model))
+
+ # Append background model
+ self._models.append(model)
+
+ # 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, 'Generate background models')
+
+ # Clear models
+ self._models.clear()
+
+ # Generate background models
+ if self['runwise'].boolean():
+ self._generate_runwise_bkg()
+ else:
+ model = self._generate_bkg(self.obs())
+ self._models.append(model)
+
+ # Write models
+ self._log_models(gammalib.NORMAL, self._models, 'Output models')
+
+ # Return
+ return
+
+ def save(self):
+ """
+ Save models
+ """
+ # Write header
+ self._log_header1(gammalib.TERSE, 'Save models')
+
+ # Get models filename
+ outmodel = self['outmodel'].filename()
+
+ # Log file name
+ self._log_value(gammalib.NORMAL, 'Models file', outmodel.url())
+
+ # Save models
+ self._models.save(outmodel)
+
+ # Return
+ return
+
+ def models(self):
+ """
+ Return background models
+
+ Returns
+ -------
+ ebounds : `~gammalib.GModels()`
+ Background models
+ """
+ # Return
+ return self._models
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Create instance of application
+ app = csbkgmodel(sys.argv)
+
+ # Execute application
+ app.execute()
diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst
index ca63f7b..52d434c 100644
--- a/doc/source/admin/release_history/1.6.rst
+++ b/doc/source/admin/release_history/1.6.rst
@@ -10,6 +10,7 @@ ctools 1.6 is a major release that adds significant functionality.
In particular, this release provides:
+* Support for H.E.S.S. public data release analysis
* Support for Xspec spectral analysis
* Improved computation speed for :ref:`ctskymap` and :ref:`csphagen`
* Parallelisation of various cscripts
@@ -89,6 +90,8 @@ Improvements
New features
------------
+* [`2694 `_] -
+ Add :ref:`csbkgmodel` script
* [`2421 `_] -
Parallelise :ref:`cslightcrv`, :ref:`csphasecrv`, :ref:`csspec`, :ref:`cssens`,
:ref:`csphagen`, :ref:`cspull` and :ref:`cstsdist`
diff --git a/doc/source/users/reference_manual/csbkgmodel.rst b/doc/source/users/reference_manual/csbkgmodel.rst
new file mode 100755
index 0000000..a00fb8e
--- /dev/null
+++ b/doc/source/users/reference_manual/csbkgmodel.rst
@@ -0,0 +1,86 @@
+.. _csbkgmodel:
+
+csbkgmodel
+==========
+
+Generates background model for 3D analysis.
+
+
+Synopsis
+--------
+
+This script generates a background model for a 3D analysis.
+
+:ref:`csbkgmodel` writes an model definition XML file.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+ Event list, counts cube or input observation definition XML file.
+
+``outmodel [file]``
+ Output model XML file.
+
+``instrument [string]``
+ Instrument name for which background models should be generated. If the
+ input observation definition XML contains only observations for a given
+ instrument, the instrument name will be extracted from the observation
+ definition XML file.
+
+``spatial [string]``
+ Spatial model component. So far IRF template, effective area and radial
+ Gaussian models are supported.
+
+``gradient [boolean]``
+ Allow for a spatial gradient in the background event distribution?
+ This option only applies to the ``GAUSS`` spatial model component.
+
+``spectral [string]``
+ Spectral model component. ``PLAW`` specifies a simple power law model,
+ ``HESS`` specifies a piecewise broken power-law using 8 energy nodes
+ spread over the energy range from 100 GeV to 12.8 TeV.
+
+``runwise [boolean]``
+ Generate runwise background model? If ``yes`` is specified a background
+ model component for each observation in the observation definition XML
+ file will be added to the output model XML file.
+
+``(rad = 2.0) [real]``
+ Radius for event selection. Only events within this radius around the
+ pointing direction will be used for model fitting.
+
+
+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 = csbkgmodel.log) [filename]``
+ Log filename.
+
+
+Related tools or scripts
+------------------------
+
+:doc:`ctlike`
diff --git a/doc/source/users/reference_manual/csmodelinfo.rst b/doc/source/users/reference_manual/csmodelinfo.rst
index b4764a3..9a6c567 100755
--- a/doc/source/users/reference_manual/csmodelinfo.rst
+++ b/doc/source/users/reference_manual/csmodelinfo.rst
@@ -28,28 +28,28 @@ General parameters
------------------
``inmodel [file]``
- Event list, counts cube, or observation definition file
+ Input model XML file.
``outds9file [file]``
- Output DS9 region file containing soucre positions
+ Output DS9 region file containing soucre positions.
``(pnt_type=cross) [string]``
- Marker type for point sources
+ Marker type for point sources.
``(pnt_mark_size=12) [integer]``
- Marker size for point sources
+ Marker size for point sources.
``(show_labels=yes) [boolean]``
Add source labels?
``(width=2) [integer]``
- Line width for regions
+ Line width for regions.
``(fontfamily=helvetica) [string]``
- Font for source labels
+ Font for source labels.
``(fontsize=12) [integer]``
- Font size for source labels
+ Font size for source labels.
``(fontweight=normal) [string]``
Use normal or bold font?
@@ -61,10 +61,10 @@ General parameters
Show type of extended model in source name?
``(free_color=green) [string]``
- Color for sources with free parameters (any ds9 color or hex code)
+ Color for sources with free parameters (any ds9 color or hex code).
``(fixed_color=magenta) [string]``
- Color for sources without free parameters (any ds9 color or hex code)
+ Color for sources without free parameters (any ds9 color or hex code).
Standard parameters
diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst
index 1685791..dabb9ca 100755
--- a/doc/source/users/reference_manual/index.rst
+++ b/doc/source/users/reference_manual/index.rst
@@ -43,6 +43,7 @@ cscripts
.. toctree::
:maxdepth: 1
+ csbkgmodel --- Generates background model for 3D analysis
cscaldb --- Lists available instrument response functions
csebins --- Generates energy boundaries for stacked analysis
cslightcrv --- Computes light curve