From f6150a99c46b9e50cb37806238fe047c60b7a05d Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 18 Oct 2018 15:44:53 +0200 Subject: [PATCH] Add csbkgmodel script (#2694) --- ChangeLog | 3 ++- NEWS | 9 ++++++++- cscripts/Makefile.am | 7 +++++-- cscripts/__init__.py.in | 4 +++- cscripts/csbkgmodel.par | 39 +++++++++++++++++++++++++++++++++++++++ cscripts/csbkgmodel.py | 412 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/admin/release_history/1.6.rst | 3 +++ doc/source/users/reference_manual/csbkgmodel.rst | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/reference_manual/csmodelinfo.rst | 18 +++++++++--------- doc/source/users/reference_manual/index.rst | 1 + 10 files changed, 568 insertions(+), 14 deletions(-) create mode 100755 cscripts/csbkgmodel.par create mode 100755 cscripts/csbkgmodel.py create mode 100755 doc/source/users/reference_manual/csbkgmodel.rst 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 -- libgit2 0.21.4