Karim Sabri / ctools

Commit add26053feaa6b311245502b8ef9eef3351a5ac9

Authored by Jürgen Knödlseder
2 parents 275d27c8 b17eaf02

Merge branch '2208-create-csphagen' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	doc/source/admin/release_history/1.5.rst
Inline Side-by-side
Showing 16 changed files with 2078 additions and 3 deletions
... ... @@ -15,4 +15,5 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
15 15 Stefan Ohm (stefan.ohm81@gmail.com)
16 16 David Sanchez (david.sanchez@lapp.in2p3.fr)
17 17 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
  18 +Luigi Tibaldo (luigi.tibaldo@irap.omp.eu)
18 19 Alexander Ziegler (alexander.ziegler@fau.de)
... ...
ChangeLog
  View file @add2605
1   -2017-10-11
  1 +2017-10-16
2 2  
3 3 * Version 1.5.0 released
4 4 ========================
5 5  
  6 + Add csphagen script (#2208)
6 7 Add keywords to sky maps produced by ctskymap (#2194)
7 8 Enhance efficiency of computations in ctmodel (#2189)
8 9 Add csmodelsois script (#2180)
... ...
1 1 New Features and Important Changes in ctools 1.5.0
2 2  
3   -11 October 2017
  3 +16 October 2017
4 4  
5 5  
6 6 Introduction
... ... @@ -132,6 +132,9 @@ a diffuse map cube (#2180).
132 132 Use t-type parameters for "tmin" and "tmax" to all specifying times in various
133 133 formats (#1864).
134 134  
  135 +The script csphagen was added to allow the generation of PHA, ARF and RMF datasets
  136 +for classical On/Off IACT analysis (#2208).
  137 +
135 138  
136 139 Examples
137 140 --------
... ...
cscripts/Makefile.am
  View file @add2605
... ... @@ -33,6 +33,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
33 33 $(srcdir)/csobsdef.py \
34 34 $(srcdir)/csobsinfo.py \
35 35 $(srcdir)/csobsselect.py \
  36 + $(srcdir)/csphagen.py \
36 37 $(srcdir)/csphasecrv.py \
37 38 $(srcdir)/cspull.py \
38 39 $(srcdir)/csresmap.py \
... ... @@ -66,6 +67,7 @@ parfiles = $(srcdir)/cscaldb.par \
66 67 $(srcdir)/csobsdef.par \
67 68 $(srcdir)/csobsinfo.par \
68 69 $(srcdir)/csobsselect.par \
  70 + $(srcdir)/csphagen.par \
69 71 $(srcdir)/csphasecrv.par \
70 72 $(srcdir)/cspull.par \
71 73 $(srcdir)/csresmap.par \
... ... @@ -96,6 +98,7 @@ linkscripts = cscaldb \
96 98 csobsdef \
97 99 csobsinfo \
98 100 csobsselect \
  101 + csphagen \
99 102 csphasecrv \
100 103 cspull \
101 104 csresmap \
... ... @@ -131,6 +134,7 @@ testscripts = $(srcdir)/tests/__init__.py \
131 134 $(top_srcdir)/test/test_csobsdef.py \
132 135 $(top_srcdir)/test/test_csobsinfo.py \
133 136 $(top_srcdir)/test/test_csobsselect.py \
  137 + $(top_srcdir)/test/test_csphagen.py \
134 138 $(top_srcdir)/test/test_csphasecrv.py \
135 139 $(top_srcdir)/test/test_cspull.py \
136 140 $(top_srcdir)/test/test_csresmap.py \
... ... @@ -159,6 +163,9 @@ testdata = $(top_srcdir)/test/data/crab.xml \
159 163 $(top_srcdir)/test/data/crab_psfcube.fits.gz \
160 164 $(top_srcdir)/test/data/crab_edispcube.fits.gz \
161 165 $(top_srcdir)/test/data/crab_bkgcube.fits.gz \
  166 + $(top_srcdir)/test/data/crab_offaxis1.fits \
  167 + $(top_srcdir)/test/data/crab_offaxis2.fits \
  168 + $(top_srcdir)/test/data/crab_exclusion.fits \
162 169 $(top_srcdir)/test/data/tsmap_0.fits \
163 170 $(top_srcdir)/test/data/tsmap_1.fits \
164 171 $(top_srcdir)/test/data/exclusion.reg \
... ...
cscripts/__init__.py.in
  View file @add2605
... ... @@ -36,6 +36,7 @@ __all__ = [
36 36 "csobsdef",
37 37 "csobsinfo",
38 38 "csobsselect",
  39 + "csphagen",
39 40 "csphasecrv",
40 41 "cspull",
41 42 "csresmap",
... ... @@ -72,6 +73,7 @@ from cscripts.csobs2caldb import csobs2caldb
72 73 from cscripts.csobsdef import csobsdef
73 74 from cscripts.csobsinfo import csobsinfo
74 75 from cscripts.csobsselect import csobsselect
  76 +from cscripts.csphagen import csphagen
75 77 from cscripts.csphasecrv import csphasecrv
76 78 from cscripts.cspull import cspull
77 79 from cscripts.csresmap import csresmap
... ...
cscripts/csobsselect.py
  View file @add2605
... ... @@ -320,7 +320,7 @@ class csobsselect(ctools.cscript):
320 320 # Otherwise log filename and save file
321 321 else:
322 322 # Log filename
323   - self._log_value(gammalib.NORMAL, 'Observation definition XML file',
  323 + self._log_value(gammalib.NORMAL, 'Obs. definition XML file',
324 324 outobs.url())
325 325  
326 326 # Save observations
... ...
cscripts/csphagen.par 0 → 100644
  View file @add2605
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csphagen.par #
  8 +# #
  9 +# Author: Luigi Tibaldo #
  10 +# luigi.tibaldo@irap.omp.eu #
  11 +# IRAP #
  12 +# #
  13 +# Purpose: Parameter file for the csphagen script #
  14 +# #
  15 +#############################################################################
  16 +#
  17 +# Input and output parameters
  18 +#============================
  19 +inobs, f, a, obs.xml,,, "Input event list or observation definition XML file"
  20 +caldb, s, a, prod2,,, "Calibration database"
  21 +irf, s, a, South_0.5h,,, "Instrument response function"
  22 +inexclusion, f, h, NONE,,, "Sky map of exclusion regions"
  23 +outobs, s, a, onoff_obs.xml,,, "Output observation definition XML file"
  24 +prefix, s, h, onoff,,, "Prefix for output files"
  25 +
  26 +#
  27 +# Script parameters
  28 +#==================
  29 +emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  30 +emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  31 +enumbins, i, a, 120,,, "Number of energy bins"
  32 +ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm"
  33 +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
  34 +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region"
  35 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
  36 +ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)"
  37 +dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)"
  38 +glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
  39 +glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
  40 +rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
  41 +bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
  42 +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
  43 +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
  44 +stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?"
  45 +
  46 +#
  47 +# Standard parameters
  48 +#====================
  49 +chatter, i, h, 2,0,4, "Chattiness of output"
  50 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  51 +debug, b, h, no,,, "Debugging mode activated"
  52 +mode, s, h, ql,,, "Mode of automatic parameters"
  53 +logfile, f, h, csphagen.log,,, "Log filename"
... ...
cscripts/csphagen.py 0 → 100755
  View file @add2605
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Computes the PHA spectra for source/background and ARF/RMF files using the
  4 +# reflected region method
  5 +#
  6 +# Copyright (C) 2017 Luigi Tibaldo
  7 +#
  8 +# This program is free software: you can redistribute it and/or modify
  9 +# it under the terms of the GNU General Public License as published by
  10 +# the Free Software Foundation, either version 3 of the License, or
  11 +# (at your option) any later version.
  12 +#
  13 +# This program is distributed in the hope that it will be useful,
  14 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +# GNU General Public License for more details.
  17 +#
  18 +# You should have received a copy of the GNU General Public License
  19 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +#
  21 +# ==========================================================================
  22 +import gammalib
  23 +import ctools
  24 +import math
  25 +import sys
  26 +
  27 +
  28 +# =============== #
  29 +# csfindobs class #
  30 +# =============== #
  31 +class csphagen(ctools.cscript):
  32 + """
  33 + Generate PHA, ARF and RMF files for classical IACT spectral analysis
  34 + """
  35 +
  36 + # Constructor
  37 + def __init__(self, *argv):
  38 + """
  39 + Constructor
  40 + """
  41 + # Set name, version
  42 + self._name = 'csphagen'
  43 + self._version = ctools.__version__
  44 +
  45 + # Initialise observation container from constructor arguments
  46 + self._obs, argv = self._set_input_obs(argv)
  47 +
  48 + # Initialise application by calling the appropriate class constructor
  49 + self._init_cscript(argv)
  50 +
  51 + # Initialise other variables
  52 + self._outobs = gammalib.GObservations()
  53 + self._ebounds = gammalib.GEbounds()
  54 + self._src_dir = gammalib.GSkyDir()
  55 + self._src_reg = gammalib.GSkyRegions()
  56 + self._bkg_regs = []
  57 + self._excl_reg = None
  58 + self._has_exclusion = False
  59 +
  60 + # Return
  61 + return
  62 +
  63 + # Private methods
  64 + def _get_parameters(self):
  65 + """
  66 + Get parameters from parfile and setup observations
  67 + """
  68 + # Setup observations (require response and allow event list, don't
  69 + # allow counts cube)
  70 + self._setup_observations(self._obs, True, True, False)
  71 +
  72 + # Set energy bounds
  73 + self._ebounds = self._create_ebounds()
  74 +
  75 + # Initialise source position/region querying relevant parameters
  76 + self._src_dir = gammalib.GSkyDir()
  77 + self._src_reg = gammalib.GSkyRegions()
  78 + coordsys = self['coordsys'].string()
  79 + if coordsys == 'CEL':
  80 + ra = self['ra'].real()
  81 + dec = self['dec'].real()
  82 + self._src_dir.radec_deg(ra, dec)
  83 + elif coordsys == 'GAL':
  84 + glon = self['glon'].real()
  85 + glat = self['glat'].real()
  86 + self._src_dir.lb_deg(glon, glat)
  87 + self._srcshape = self['srcshape'].string()
  88 + if self._srcshape == 'CIRCLE':
  89 + self._rad = self['rad'].real()
  90 + self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
  91 +
  92 + # Query background estimation method and parameters
  93 + bkgmethod = self['bkgmethod'].string()
  94 + if bkgmethod == 'REFLECTED':
  95 + self['bkgregmin'].integer()
  96 + self['maxoffset'].real()
  97 +
  98 + # Exclusion map
  99 + if self['inexclusion'].filename().is_fits():
  100 + self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename())
  101 + self._has_exclusion = True
  102 + else:
  103 + self._has_exclusion = False
  104 +
  105 + # Query remaining parametersStacking
  106 + self['stack'].boolean()
  107 +
  108 + # Query ahead output parameters
  109 + if (self._read_ahead()):
  110 + self['outobs'].string()
  111 + self['prefix'].string()
  112 +
  113 + # If there are no observations in container then get them from the
  114 + # parameter file
  115 + if self._obs.size() == 0:
  116 + self._obs = self._get_observations(False)
  117 +
  118 + # Write input parameters into logger
  119 + self._log_parameters(gammalib.TERSE)
  120 +
  121 + # Return
  122 + return
  123 +
  124 + def _reflected_regions(self, obs):
  125 + """
  126 + Calculate list of reflected regions for a single observation (pointing)
  127 +
  128 + Parameters
  129 + ----------
  130 + obs : `~gammalib.GObservation()`
  131 + CTA observation
  132 +
  133 + Returns
  134 + -------
  135 + regions : `~gammalib.GSkyRegions`
  136 + List of reflected regions
  137 + """
  138 + # Initialise list of reflected regions
  139 + regions = gammalib.GSkyRegions()
  140 +
  141 + # Get offset angle of source
  142 + pnt_dir = obs.pointing().dir()
  143 + offset = pnt_dir.dist_deg(self._src_dir)
  144 +
  145 + # If ...
  146 + if offset <= self._rad or offset >= self['maxoffset'].real():
  147 + self._log_string(gammalib.EXPLICIT, 'Observation {} pointed at {} '
  148 + 'deg from source'.format(obs.id(), offset))
  149 +
  150 + # ... otherwise
  151 + else:
  152 + posang = pnt_dir.posang_deg(self._src_dir)
  153 + if self._srcshape == 'CIRCLE':
  154 + # angular separation of reflected regions wrt camera center
  155 + # and number
  156 + alpha = 1.05 * 2 * self._rad / offset
  157 + # 1.05 ensures background regions do not overlap due to
  158 + # numerical precision issues
  159 + N = int(2.0 * math.pi / alpha)
  160 + if N < self['bkgregmin'].integer() + 3:
  161 + self._log_string(gammalib.EXPLICIT, 'Observation {}: '
  162 + 'insufficient regions for background '
  163 + 'estimation'.format(obs.id(), offset))
  164 +
  165 + # Otherwise loop over position angle to create reflected
  166 + # regions
  167 + else:
  168 + alpha = 360.0 / N
  169 + for s in range(2, N - 1):
  170 + dphi = s * alpha
  171 + ctr_dir = pnt_dir.clone()
  172 + ctr_dir.rotate_deg(posang + dphi, offset)
  173 + region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
  174 + if self._has_exclusion:
  175 + if self._excl_reg.overlaps(region):
  176 + self._log_string(gammalib.VERBOSE, 'Observation '
  177 + '{}: reflected region overlaps '
  178 + 'with exclusion region'.
  179 + format(obs.id(), offset))
  180 + else:
  181 + regions.append(region)
  182 + else:
  183 + regions.append(region)
  184 +
  185 + # Return reflected regions
  186 + return regions
  187 +
  188 + # Public methods
  189 + def run(self):
  190 + """
  191 + Run the script
  192 + """
  193 + # Switch screen logging on in debug mode
  194 + if self._logDebug():
  195 + self._log.cout(True)
  196 +
  197 + # Get parameters
  198 + self._get_parameters()
  199 +
  200 + # Write observation into logger
  201 + self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
  202 +
  203 + # Write header
  204 + self._log_header1(gammalib.NORMAL,
  205 + 'Generation of source and background spectra')
  206 +
  207 + # Initialise run variables
  208 + self._outobs = gammalib.GObservations()
  209 + self._bkg_regs = []
  210 +
  211 + # Loop through observations and generate pha, arf, rmf files
  212 + for obs in self._obs:
  213 +
  214 + # Initialise background regions for this observation
  215 + bkg_reg = gammalib.GSkyRegions()
  216 +
  217 + # If reflected background is requested then created reflected
  218 + # background regions
  219 + if self['bkgmethod'].string() == 'REFLECTED':
  220 + bkg_reg = self._reflected_regions(obs)
  221 +
  222 + # If there are reflected regions then create On/Off observation
  223 + # and append it to the output container
  224 + if bkg_reg.size() >= self['bkgregmin'].integer():
  225 + onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir,
  226 + self._ebounds,
  227 + self._ebounds,
  228 + self._src_reg,
  229 + bkg_reg)
  230 + onoff.id(obs.id())
  231 + self._outobs.append(onoff)
  232 + self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()})
  233 + else:
  234 + self._log_string(gammalib.NORMAL, 'Observation {} not included '
  235 + 'in spectra generation'.format(obs.id()))
  236 +
  237 + # Stack observations
  238 + if self._outobs.size() > 1 and self['stack'].boolean() == True:
  239 +
  240 + # Write header
  241 + self._log_header1(gammalib.NORMAL,
  242 + 'Stacking {} observations'.format(self._outobs.size()))
  243 +
  244 + # Stack observations
  245 + stacked_obs = gammalib.GCTAOnOffObservation(self._outobs)
  246 +
  247 + # Put stacked observations in output container
  248 + self._outobs = gammalib.GObservations()
  249 + self._outobs.append(stacked_obs)
  250 +
  251 + # Return
  252 + return
  253 +
  254 + def save(self):
  255 + """
  256 + Save data
  257 + """
  258 + # Write header
  259 + self._log_header1(gammalib.TERSE, 'Save data')
  260 +
  261 + # Get XMK output filename, prefix and clobber
  262 + outobs = self['outobs'].string()
  263 + prefix = self['prefix'].string()
  264 + clobber = self['clobber'].boolean()
  265 +
  266 + # Loop over all observation in container
  267 + for obs in self._outobs:
  268 +
  269 + # Set filenames
  270 + if self['stack'].boolean():
  271 + onname = prefix + '_stacked_pha_on.fits'
  272 + offname = prefix + '_stacked_pha_off.fits'
  273 + arfname = prefix + '_stacked_arf.fits'
  274 + rmfname = prefix + '_stacked_rmf.fits'
  275 + elif self._outobs.size() > 1:
  276 + onname = prefix + '_{}_pha_on.fits'.format(obs.id())
  277 + offname = prefix + '_{}_pha_off.fits'.format(obs.id())
  278 + arfname = prefix + '_{}_arf.fits'.format(obs.id())
  279 + rmfname = prefix + '_{}_rmf.fits'.format(obs.id())
  280 + else:
  281 + onname = prefix + '_pha_on.fits'
  282 + offname = prefix + '_pha_off.fits'
  283 + arfname = prefix + '_arf.fits'
  284 + rmfname = prefix + '_rmf.fits'
  285 +
  286 + # Save files
  287 + obs.on_spec().save(onname, clobber)
  288 + obs.off_spec().save(offname, clobber)
  289 + obs.arf().save(arfname, clobber)
  290 + obs.rmf().save(rmfname, clobber)
  291 +
  292 + # Log file names
  293 + self._log_value(gammalib.NORMAL, 'PHA on file', onname)
  294 + self._log_value(gammalib.NORMAL, 'PHA off file', offname)
  295 + self._log_value(gammalib.NORMAL, 'ARF file', arfname)
  296 + self._log_value(gammalib.NORMAL, 'RMF file', arfname)
  297 +
  298 + # Save observation definition XML file
  299 + self._outobs.save(outobs)
  300 +
  301 + # Log file name
  302 + self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs)
  303 +
  304 + # Save ds9 On region file
  305 + regname = prefix + '_on.reg'
  306 + self._src_reg.save(regname)
  307 + self._log_value(gammalib.NORMAL, 'On region file', regname)
  308 +
  309 + # Save ds9 Off region files
  310 + for bkg_reg in self._bkg_regs:
  311 +
  312 + # Set filename
  313 + if len(self._bkg_regs) > 1:
  314 + regname = prefix + '_{}_off.reg'.format(bkg_reg['id'])
  315 + else:
  316 + regname = prefix + '_off.reg'
  317 +
  318 + # Save ds9 region file
  319 + bkg_reg['regions'].save(regname)
  320 +
  321 + # Log file name
  322 + self._log_value(gammalib.NORMAL, 'Off region file', regname)
  323 +
  324 + # Return
  325 + return
  326 +
  327 + def obs(self):
  328 + """
  329 + Return observation container
  330 + """
  331 + return self._outobs
  332 +
  333 + def execute(self):
  334 + """
  335 + Execute the script
  336 + """
  337 + # Open logfile
  338 + self.logFileOpen()
  339 +
  340 + # Read ahead output parameters
  341 + self._read_ahead(True)
  342 +
  343 + # Run the script
  344 + self.run()
  345 +
  346 + # Save result
  347 + self.save()
  348 +
  349 + # Return
  350 + return
  351 +
  352 +
  353 +# ======================== #
  354 +# Main routine entry point #
  355 +# ======================== #
  356 +if __name__ == '__main__':
  357 +
  358 + # Create instance of application
  359 + app = csphagen(sys.argv)
  360 +
  361 + # Execute application
  362 + app.execute()
... ...
doc/source/admin/release_history/1.5.rst
  View file @add2605
... ... @@ -10,6 +10,7 @@ ctools 1.5 is a major release that adds significant functionality.
10 10  
11 11 In particular, this release provides:
12 12  
  13 +* Support of classical IACT analysis
13 14 * Specifcation of times in various formats (UTC, MJD, JD, MET)
14 15  
15 16  
... ... @@ -31,6 +32,8 @@ Improvements
31 32 New features
32 33 ------------
33 34  
  35 +* [`2208 <https://cta-redmine.irap.omp.eu/issues/2208>`_] -
  36 + Add :ref:`csphagen` script
34 37 * [`2180 <https://cta-redmine.irap.omp.eu/issues/2180>`_] -
35 38 Add :ref:`csmodelsois` script
36 39 * [`1864 <https://cta-redmine.irap.omp.eu/issues/1864>`_] -
... ...
doc/source/users/reference_manual/csphagen.rst 0 → 100644
  View file @add2605
  1 +.. _csphagen:
  2 +
  3 +csphagen
  4 +========
  5 +
  6 +Creates files necessary to perform a classical region-based spectral analysis
  7 +on IACT data, i.e., source and background PHA files, ARF and RMF files.
  8 +
  9 +
  10 +Synopsis
  11 +--------
  12 +
  13 +This script can be used to derive from an observation or a set of observations
  14 +the count spectra in a source region and in background regions, as well as the
  15 +detector response (effective area, energy redistribution matrix), that can be
  16 +used to perform a classical 1D spectral analysis. The output files are saved in
  17 +the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF).
  18 +`See here <https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html>`__.
  19 +
  20 +The outputs are:
  21 +1) the PHA, ARF, RMF files, either separately for each observation, or stacked
  22 + for a set of observations,
  23 +2) DS9 regions files listing the source and background regions for each
  24 + observation;
  25 +3) a new observation definition XML file.
  26 +
  27 +
  28 +General parameters
  29 +------------------
  30 +
  31 +``inobs [file]``
  32 + Input event list or observation definition XML file.
  33 +
  34 +``caldb [string]``
  35 + Calibration database.
  36 +
  37 +``irf [string]``
  38 + Instrument response function.
  39 +
  40 +``(inexclusion = NONE) [file]``
  41 + Optional FITS file containing a WCS map in the first hdu that defines sky
  42 + regions not to be used for background estimation (where map value != 0).
  43 +
  44 +``outobs [string]``
  45 + Output observation definition XML file.
  46 +
  47 +``(prefix = onoff) [string]``
  48 + Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files.
  49 +
  50 +``emin [real]``
  51 + Lower energy limit (TeV) if LIN or LOG binning algorithms are used.
  52 +
  53 +``emax [real]``
  54 + Upper energy limit (TeV) if LIN or LOG binning algorithms are used.
  55 +
  56 +``enumbins [integer]``
  57 + Number of energy bins. At least 30 bins per decade are recommended for
  58 + proper evaluation of the instrument response.
  59 +
  60 +``ebinalg [string]’’
  61 + Energy binning algorithm: LOG - logarithmically spaced energy bins;
  62 + LIN - linearly spaced energy bins; FILE - energy bounds retrieved from file.
  63 +
  64 +``ebinfile [file]``
  65 + Name of the file containing the energy bin definition if FILE algorithm
  66 + is used.
  67 +
  68 +``(srcshape = CIRCLE) [string]``
  69 + Shape of the source region.
  70 + CIRCLE: circular region around given position.
  71 +
  72 +``coordsys [string]``
  73 + Coordinate system (CEL - celestial, GAL - galactic)
  74 +
  75 +``ra [real]``
  76 + Right Ascension of source region centre (deg)
  77 +
  78 +``dec [real]``
  79 + Declination of source region centre (deg)
  80 +
  81 +``glon [real]``
  82 + Galactic longitude of source region centre (deg)
  83 +
  84 +``glat [real]``
  85 + Galactic latitude of source region centre (deg)
  86 +
  87 +``rad [real]``
  88 + Radius of source region circle (deg)
  89 +
  90 +``(bkgmethod = REFLECTED) [string]’’
  91 + Method for background estimation.
  92 + REFLECTED: background evaluated in regions with the same shape as source
  93 + region reflected w.r.t. pointing direction for each observation.
  94 +
  95 +``(bkgregmin = 2) [integer]``
  96 + Minimum number of background regions that are required for an observation.
  97 + If this number of background regions is not available the observation is
  98 + skipped.
  99 +
  100 +``(maxoffset = 4.0) [real]``
  101 + Maximum offset in degrees of source from camera center to accept the
  102 + observation.
  103 +
  104 +``stack [boolean]``
  105 + Specifies whether multiple observations should be stacked (``yes``) or
  106 + whether run-wise PHA, ARF and RMF files should be produced (``no``).
  107 +
  108 +
  109 +Standard parameters
  110 +-------------------
  111 +
  112 +``(chatter = 2) [integer]``
  113 + Verbosity of the executable:
  114 + ``chatter = 0``: no information will be logged
  115 +
  116 + ``chatter = 1``: only errors will be logged
  117 +
  118 + ``chatter = 2``: errors and actions will be logged
  119 +
  120 + ``chatter = 3``: report about the task execution
  121 +
  122 + ``chatter = 4``: detailed report about the task execution
  123 +
  124 +``(clobber = yes) [boolean]``
  125 + Specifies whether an existing output runlist should be overwritten.
  126 +
  127 +``(debug = no) [boolean]``
  128 + Enables debug mode. In debug mode the executable will dump any log file
  129 + output to the console.
  130 +
  131 +``(mode = ql) [string]``
  132 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
  133 +
  134 +``(logfile = csphagen.log) [filename]``
  135 + Log filename.
  136 +
  137 +
... ...
doc/source/users/reference_manual/index.rst
  View file @add2605
... ... @@ -50,6 +50,7 @@ cscripts
50 50 csobsdef --- Generates observation definition file <csobsdef>
51 51 csobsinfo --- Shows observation container content <csobsinfo>
52 52 csobsselect --- Select observations from observation definition file <csobsselect>
  53 + csphagen —-- Generates PHA, ARF, RMF files based on source/background regions <csphagen>
53 54 csphasecrv --- Computes phase curve <csphasecrv>
54 55 cspull --- Generates pull distribution <cspull>
55 56 csresmap --- Generates residual map <csresmap>
... ...
test/data/crab_exclusion.fits 0 → 100644
View file @add2605
No preview for this file type
test/data/crab_offaxis1.fits 0 → 100644
View file @add2605
No preview for this file type
test/data/crab_offaxis2.fits 0 → 100644
View file @add2605
No preview for this file type
test/test_csphagen.py 0 → 100644
  View file @add2605
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the csphagen script.
  4 +#
  5 +# Copyright (C) 2017 Luigi Tibaldo
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import cscripts
  23 +from testing import test
  24 +
  25 +
  26 +# ================================ #
  27 +# Test class for csphagen script #
  28 +# ================================ #
  29 +class Test(test):
  30 + """
  31 + Test class for csphagen script
  32 +
  33 + This test class makes unit tests for the csphagen script by using it
  34 + from the command line and from Python.
  35 + """
  36 +
  37 + # Constructor
  38 + def __init__(self):
  39 + """
  40 + Constructor
  41 + """
  42 + # Call base class constructor
  43 + test.__init__(self)
  44 + self._myevents1 = self._datadir + '/crab_offaxis1.fits'
  45 + self._myevents2 = self._datadir + '/crab_offaxis2.fits'
  46 + self._exclusion = self._datadir + '/crab_exclusion.fits'
  47 + self._nreg_with_excl = 5
  48 + self._nreg_wo_excl = 8
  49 + self._nreg_mul = [self._nreg_with_excl, 6]
  50 + self._nbins = 20
  51 + # number of expected background regions with/wo exclusion,
  52 + # and for two different runs
  53 +
  54 + # Return
  55 + return
  56 +
  57 + # Set test functions
  58 + def set(self):
  59 + """
  60 + Set all test functions.
  61 + """
  62 + # Set test name
  63 + self.name('csphagen')
  64 +
  65 + # Append tests
  66 + self.append(self._test_cmd, 'Test csphagen on command line')
  67 + self.append(self._test_python, 'Test csphagen from Python')
  68 +
  69 + # Return
  70 + return
  71 +
  72 + # Test csphagen on command line
  73 + def _test_cmd(self):
  74 + """
  75 + Test csphagen on the command line.
  76 + """
  77 + # Set script name
  78 + csphagen = self._script('csphagen')
  79 +
  80 + # Setup csphagen command
  81 + cmd = csphagen + ' inobs="' + self._myevents1 + \
  82 + '" caldb="' + self._caldb + '" irf="' + self._irf + \
  83 + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
  84 + str(self._nbins) + \
  85 + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
  86 + ' rad=0.2 stack="no" inexclusion="' + \
  87 + self._exclusion + \
  88 + '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \
  89 + 'logfile="csphagen_cmd1.log" chatter=2'
  90 +
  91 + # Check if execution of wrong command fails
  92 + self.test_assert(self._execute('command_that_does_not_exist') != 0,
  93 + 'Self test of test script')
  94 +
  95 + # Check if execution was successful
  96 + self.test_assert(self._execute(cmd) == 0,
  97 + 'Check successful execution from command line')
  98 +
  99 + # Check output files
  100 + self._check_output('csphagen_cmd1', self._nbins, self._nreg_with_excl)
  101 + self._check_outobs('csphagen_cmd1', 1)
  102 +
  103 + # Setup csphagen command
  104 + cmd = csphagen + ' inobs="events_that_do_not_exist.fits" ' + \
  105 + 'caldb="' + self._caldb + '" irf="' + self._irf + \
  106 + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \
  107 + str(self._nbins) + \
  108 + ' coordsys="CEL" ra=83.633 dec=22.0145' + \
  109 + ' rad=0.2 stack="no" inexclusion="' + \
  110 + self._exclusion + \
  111 + '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \
  112 + 'logfile="csphagen_cmd2.log" chatter=1'
  113 +
  114 + # Check if execution failed
  115 + self.test_assert(self._execute(cmd) != 0,
  116 + 'Check invalid input file when executed from command line')
  117 +
  118 + # Check csphagen --help
  119 + self._check_help(csphagen)
  120 +
  121 + # Return
  122 + return
  123 +
  124 + def _test_python(self):
  125 + """
  126 + Test csphagen from Python
  127 + """
  128 + # Same test as from command line
  129 + phagen = cscripts.csphagen()
  130 + phagen['inobs'] = self._myevents1
  131 + phagen['caldb'] = self._caldb
  132 + phagen['irf'] = self._irf
  133 + phagen['ebinalg'] = 'LOG'
  134 + phagen['emin'] = 0.1
  135 + phagen['emax'] = 100.0
  136 + phagen['enumbins'] = self._nbins
  137 + phagen['coordsys'] = 'CEL'
  138 + phagen['ra'] = 83.633
  139 + phagen['dec'] = 22.0145
  140 + phagen['rad'] = 0.2
  141 + phagen['stack'] = False
  142 + phagen['inexclusion'] = self._exclusion
  143 + phagen['outobs'] = 'csphagen_py1.xml'
  144 + phagen['prefix'] = 'csphagen_py1'
  145 + phagen['logfile'] = 'csphagen_py1.log'
  146 + phagen['chatter'] = 1
  147 +
  148 + # Run script
  149 + phagen.execute()
  150 +
  151 + # Check output
  152 + self._check_output('csphagen_py1', self._nbins, self._nreg_with_excl)
  153 + self._check_outobs('csphagen_py1', 1)
  154 +
  155 + # Now test without exclusion region
  156 + phagen = cscripts.csphagen()
  157 + phagen['inobs'] = self._myevents1
  158 + phagen['caldb'] = self._caldb
  159 + phagen['irf'] = self._irf
  160 + phagen['ebinalg'] = 'LOG'
  161 + phagen['emin'] = 0.1
  162 + phagen['emax'] = 100.0
  163 + phagen['enumbins'] = self._nbins
  164 + phagen['coordsys'] = 'CEL'
  165 + phagen['ra'] = 83.633
  166 + phagen['dec'] = 22.0145
  167 + phagen['rad'] = 0.2
  168 + phagen['stack'] = False
  169 + phagen['outobs'] = 'csphagen_py2.xml'
  170 + phagen['prefix'] = 'csphagen_py2'
  171 + phagen['logfile'] = 'csphagen_py2.log'
  172 + phagen['chatter'] = 2
  173 +
  174 + # Run script
  175 + phagen.execute()
  176 +
  177 + # Check output
  178 + self._check_output('csphagen_py2', self._nbins, self._nreg_wo_excl)
  179 + self._check_outobs('csphagen_py2', 1)
  180 +
  181 + # Test with multiple input observations, no stacking
  182 + # Create observation container
  183 + obs = gammalib.GObservations()
  184 + for s, events in enumerate([self._myevents1, self._myevents2]):
  185 + run = gammalib.GCTAObservation(events)
  186 + run.id(str(s + 1))
  187 + run.response(self._irf, gammalib.GCaldb('cta', self._caldb))
  188 + obs.append(run)
  189 +
  190 + # Setup csphagen
  191 + phagen = cscripts.csphagen(obs)
  192 + phagen['ebinalg'] = 'LOG'
  193 + phagen['emin'] = 0.1
  194 + phagen['emax'] = 100.0
  195 + phagen['enumbins'] = self._nbins
  196 + phagen['coordsys'] = 'CEL'
  197 + phagen['ra'] = 83.633
  198 + phagen['dec'] = 22.0145
  199 + phagen['rad'] = 0.2
  200 + phagen['stack'] = False
  201 + phagen['inexclusion'] = self._exclusion
  202 + phagen['outobs'] = 'csphagen_py3.xml'
  203 + phagen['prefix'] = 'csphagen_py3'
  204 + phagen['logfile'] = 'csphagen_py3.log'
  205 + phagen['chatter'] = 3
  206 +
  207 + # Run script
  208 + phagen.execute()
  209 +
  210 + # Check output
  211 + for s in range(2):
  212 + self._check_output('csphagen_py3_' + str(s + 1), self._nbins,
  213 + self._nreg_mul[s])
  214 + self._check_outobs('csphagen_py3', 2)
  215 +
  216 + # Test with multiple input observations and stacking
  217 +
  218 + # Setup csphagen
  219 + phagen = cscripts.csphagen(obs)
  220 + phagen['ebinalg'] = 'LOG'
  221 + phagen['emin'] = 0.1
  222 + phagen['emax'] = 100.0
  223 + phagen['enumbins'] = self._nbins
  224 + phagen['coordsys'] = 'CEL'
  225 + phagen['ra'] = 83.633
  226 + phagen['dec'] = 22.0145
  227 + phagen['rad'] = 0.2
  228 + phagen['stack'] = True
  229 + phagen['inexclusion'] = self._exclusion
  230 + phagen['outobs'] = 'csphagen_py4.xml'
  231 + phagen['prefix'] = 'csphagen_py4'
  232 + phagen['logfile'] = 'csphagen_py4.log'
  233 + phagen['chatter'] = 4
  234 +
  235 + # Run script
  236 + phagen.execute()
  237 +
  238 + # Check output
  239 + for s in range(2):
  240 + self._check_output('csphagen_py4_stacked', self._nbins,
  241 + 0, check_regions=False)
  242 + self._check_outobs('csphagen_py4', 1)
  243 +
  244 + return
  245 +
  246 + def _check_ebounds(self, table, bins):
  247 + """
  248 + Check EBOUNDS table
  249 + """
  250 + cols = ['E_MIN', 'E_MAX']
  251 +
  252 + self.test_value(table.ncols(), len(cols),
  253 + 'Check for %d columns in energy bounds table' % len(
  254 + cols))
  255 + self.test_value(table.nrows(), bins,
  256 + 'Check for %d rows in energy bounds table' % bins)
  257 + for col in cols:
  258 + self.test_assert(table.contains(col),
  259 + 'FITS file contains "' + col + '" column')
  260 +
  261 + return
  262 +
  263 + def _check_pha(self, filename, bins):
  264 + """
  265 + Check PHA file
  266 + """
  267 + # Expected column names
  268 + cols = ['CHANNEL', 'COUNTS', 'STAT_ERR', 'SYS_ERR', 'QUALITY',
  269 + 'GROUPING', 'AREASCAL', 'BACKSCAL']
  270 +
  271 + # Open FITS file
  272 + fits = gammalib.GFits(filename)
  273 +
  274 + # Check FITS file structure
  275 + self.test_value(fits.size(), 3,
  276 + 'Check for 3 extensions in PHA file')
  277 + self.test_assert(fits.contains('SPECTRUM'),
  278 + 'FITS file contains "SPECTRUM" extension')
  279 + self.test_assert(fits.contains('EBOUNDS'),
  280 + 'FITS file contains "EBOUNDS" extension')
  281 +
  282 + # Get SPECTRUM table
  283 + table = fits['SPECTRUM']
  284 +
  285 + # Check FITS table structure
  286 + self.test_value(table.ncols(), len(cols),
  287 + 'Check for %d columns in PHA table' % len(cols))
  288 + self.test_value(table.nrows(), bins,
  289 + 'Check for %d rows in PHA table' % bins)
  290 + for col in cols:
  291 + self.test_assert(table.contains(col),
  292 + 'FITS file contains "' + col + '" column')
  293 +
  294 + # Check SPECTRUM table
  295 + table = fits['EBOUNDS']
  296 + self._check_ebounds(table, bins)
  297 +
  298 + # Close FITS file
  299 + fits.close()
  300 +
  301 + # Return
  302 + return
  303 +
  304 + def _check_arf(self, filename, bins):
  305 + """
  306 + Check ARF file
  307 + """
  308 + # Expected column names
  309 + cols = ['ENERG_LO', 'ENERG_HI', 'SPECRESP', 'BACKRESP']
  310 +
  311 + # Open FITS file
  312 + fits = gammalib.GFits(filename)
  313 +
  314 + # Check FITS file structure
  315 + self.test_value(fits.size(), 2,
  316 + 'Check for 2 extensions in ARF file')
  317 + self.test_assert(fits.contains('SPECRESP'),
  318 + 'FITS file contains "SPECRESP" extension')
  319 +
  320 + # Get SPECRESP table
  321 + table = fits['SPECRESP']
  322 +
  323 + # Check FITS table structure
  324 + self.test_value(table.ncols(), len(cols),
  325 + 'Check for %d columns in ARF table' % len(cols))
  326 + self.test_value(table.nrows(), bins,
  327 + 'Check for %d rows in ARF table' % bins)
  328 + for col in cols:
  329 + self.test_assert(table.contains(col),
  330 + 'FITS file contains "' + col + '" column')
  331 +
  332 + # Close FITS file
  333 + fits.close()
  334 +
  335 + # Return
  336 + return
  337 +
  338 + def _check_rmf(self, filename, bins):
  339 + """
  340 + Check RMF file
  341 + """
  342 + # Expected column names
  343 + cols = ['ENERG_LO', 'ENERG_HI', 'N_GRP', 'F_CHAN', 'N_CHAN', 'MATRIX']
  344 +
  345 + # Open FITS file
  346 + fits = gammalib.GFits(filename)
  347 +
  348 + # Check FITS file structure
  349 + self.test_value(fits.size(), 3,
  350 + 'Check for 3 extensions in RMF file')
  351 + self.test_assert(fits.contains('MATRIX'),
  352 + 'FITS file contains "MATRIX" extension')
  353 + self.test_assert(fits.contains('EBOUNDS'),
  354 + 'FITS file contains "EBOUNDS" extension')
  355 +
  356 + # Get MATRIX table
  357 + table = fits['MATRIX']
  358 +
  359 + # Check FITS table structure
  360 + self.test_value(table.ncols(), len(cols),
  361 + 'Check for %d columns in RMF table' % len(cols))
  362 + self.test_value(table.nrows(), bins,
  363 + 'Check for %d rows in RMF table' % bins)
  364 + for col in cols:
  365 + self.test_assert(table.contains(col),
  366 + 'FITS file contains "' + col + '" column')
  367 +
  368 + # Check EBOUNDS table
  369 + table = fits['EBOUNDS']
  370 + self._check_ebounds(table, bins)
  371 +
  372 + # Close FITS file
  373 + fits.close()
  374 +
  375 + # Return
  376 + return
  377 +
  378 + def _check_output(self, filenameroot, bins, nreg, check_regions=True):
  379 + """
  380 + Check the output from a csphagen run
  381 + """
  382 + # OGIP files
  383 + self._check_pha(filenameroot + '_pha_on.fits', bins)
  384 + self._check_pha(filenameroot + '_pha_off.fits', bins)
  385 + self._check_arf(filenameroot + '_arf.fits', bins)
  386 + self._check_rmf(filenameroot + '_rmf.fits', bins)
  387 +
  388 + # Optionally check for regions
  389 + if check_regions:
  390 +
  391 + # Check On region
  392 + onregion = (filenameroot+'_on.reg').replace('_1_', '_').replace('_2_', '_')
  393 + reg = gammalib.GSkyRegions(onregion)
  394 + self.test_value(reg.size(), 1, 'Check for 1 region in source region file')
  395 +
  396 + # Check for Off region
  397 + offregion = filenameroot + '_off.reg'
  398 + reg = gammalib.GSkyRegions(offregion)
  399 + self.test_value(reg.size(), nreg, 'Check for ' + str(nreg) +
  400 + ' region in background region file')
  401 +
  402 + # Return
  403 + return
  404 +
  405 + def _check_outobs(self, filenameroot, nout):
  406 + """
  407 + Check the output XML file containing ON/OFF observations
  408 + """
  409 + # Load observation container
  410 + obs = gammalib.GObservations(filenameroot + '.xml')
  411 +
  412 + # Check container size
  413 + self.test_value(obs.size(), nout, 'Check for ' + str(nout) +
  414 + ' observations in XML file')
  415 +
  416 + # Return
  417 + return
... ...
test/test_python_cscripts.py
  View file @add2605
... ... @@ -37,6 +37,7 @@ import test_csobs2caldb
37 37 import test_csobsdef
38 38 import test_csobsinfo
39 39 import test_csobsselect
  40 +import test_csphagen
40 41 import test_csphasecrv
41 42 import test_cspull
42 43 import test_csresmap
... ... @@ -127,6 +128,7 @@ def test(installed=False, debug=False):
127 128 test_csobsdef.Test(),
128 129 test_csobsinfo.Test(),
129 130 test_csobsselect.Test(),
  131 + test_csphagen.Test(),
130 132 test_csphasecrv.Test(),
131 133 test_cspull.Test(),
132 134 test_csresmap.Test(),
... ...