Add csbkgmodel script (#2694)

Jürgen Knödlseder
1 parent f4dea160
Showing 10 changed files with 568 additions and 14 deletions
ChangeLog
NEWS
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csbkgmodel.par
cscripts/csbkgmodel.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/csbkgmodel.rst
doc/source/users/reference_manual/csmodelinfo.rst
doc/source/users/reference_manual/index.rst
-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)
 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
@@ -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 \
@@ -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 &gt; 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
+#############################################################################
+#                                                                           #
+#                       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"
+#! /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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+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()
@@ -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 <https://cta-redmine.irap.omp.eu/issues/2694>`_] -
+  Add :ref:`csbkgmodel` script
 * [`2421 <https://cta-redmine.irap.omp.eu/issues/2421>`_] -
   Parallelise :ref:`cslightcrv`, :ref:`csphasecrv`, :ref:`csspec`, :ref:`cssens`,
   :ref:`csphagen`, :ref:`cspull` and :ref:`cstsdist`
+.. _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 <IRF|AEFF|GAUSS> [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 <PLAW|HESS> [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`
@@ -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) <circle|box|diamond|cross|x|arrow|boxcircle> [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) <helvetica|times|courier> [string]``
-    Font for source labels
+    Font for source labels.
 ``(fontsize=12) [integer]``
-    Font size for source labels
+    Font size for source labels.
 ``(fontweight=normal) <normal|bold> [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
@@ -43,6 +43,7 @@ cscripts
 .. toctree::
    :maxdepth: 1
+   csbkgmodel --- Generates background model for 3D analysis <csbkgmodel>
    cscaldb --- Lists available instrument response functions <cscaldb>
    csebins --- Generates energy boundaries for stacked analysis <csebins>
    cslightcrv --- Computes light curve <cslightcrv>