diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index f292d84..a21dc1a 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -37,6 +37,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
$(srcdir)/csphasecrv.py \
$(srcdir)/cspull.py \
$(srcdir)/csresmap.py \
+ $(srcdir)/csresspec.py \
$(srcdir)/csroot2caldb.py \
$(srcdir)/cssens.py \
$(srcdir)/csspec.py \
@@ -71,6 +72,7 @@ parfiles = $(srcdir)/cscaldb.par \
$(srcdir)/csphasecrv.par \
$(srcdir)/cspull.par \
$(srcdir)/csresmap.par \
+ $(srcdir)/csresspec.par \
$(srcdir)/csroot2caldb.par \
$(srcdir)/cssens.par \
$(srcdir)/csspec.par \
@@ -102,6 +104,7 @@ linkscripts = cscaldb \
csphasecrv \
cspull \
csresmap \
+ csresspec \
csroot2caldb \
cssens \
csspec \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index a007255..bc3ca41 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -1,7 +1,7 @@
# ==========================================================================
# cscripts Python module
#
-# Copyright (C) 2015-2017 Juergen Knoedlseder
+# Copyright (C) 2015-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
@@ -40,6 +40,7 @@ __all__ = [
"csphasecrv",
"cspull",
"csresmap",
+ "csresspec",
"csroot2caldb",
"cssens",
"csspec",
@@ -77,6 +78,7 @@ from cscripts.csphagen import csphagen
from cscripts.csphasecrv import csphasecrv
from cscripts.cspull import cspull
from cscripts.csresmap import csresmap
+from cscripts.csresspec import csresspec
from cscripts.csroot2caldb import csroot2caldb
from cscripts.cssens import cssens
from cscripts.csspec import csspec
diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py
index a5aab3b..5dabd68 100755
--- a/cscripts/csresmap.py
+++ b/cscripts/csresmap.py
@@ -19,9 +19,9 @@
#
# ==========================================================================
import sys
-import math
import gammalib
import ctools
+from cscripts import obsutils
# ============== #
@@ -209,67 +209,14 @@ class csresmap(ctools.csobservation):
# Get model map into GSkyMap object
modelmap = model.cube().counts().copy()
- # Store counts map as residual map. Note that we need a special
+ # Calculate residual maps
+ # Note that we need a special
# construct here to avoid memory leaks. This seems to be a SWIG feature
# as SWIG creates a new object when calling binning.cube()
- self._resmap = countmap.copy()
- self._resmap.stack_maps()
+ countmap1 = countmap.copy()
+ countmap1.stack_maps()
modelmap.stack_maps()
-
- # Get residual map algorithm type
- algorithm = self['algorithm'].string()
-
- # Subtract maps
- if algorithm == 'SUB':
- self._resmap -= modelmap
-
- # Subtract and divide by model map
- elif algorithm == 'SUBDIV':
- self._resmap -= modelmap
- self._resmap /= modelmap
-
- # Subtract and divide by sqrt of model map
- elif algorithm == 'SUBDIVSQRT':
- self._resmap -= modelmap
- self._resmap /= modelmap.sqrt()
-
- # Calculate significance from Li&Ma derivation
- elif algorithm == 'SIGNIFICANCE':
-
- # Compute sign map
- signmap = (self._resmap - modelmap).sign()
-
- # Loop over every bin in the map
- for i in range(self._resmap.npix()):
-
- # If the model value > 0.0 do the computation as normal ...
- model_val = modelmap[i]
- if model_val > 0.0:
-
- # If the data value is also > 0 then compute the
- # significance^2 and save it to the map ...
- data_val = self._resmap[i]
- if data_val > 0.0:
- log_val = math.log(data_val/model_val)
- self._resmap[i] = (data_val*log_val) + model_val - data_val
-
- # ... otherwise compute the reduced value of the above
- # expression. This is necessary to avoid computing log(0).
- else:
- self._resmap[i] = model_val
-
- # ... otherwise hard-code the significance to 0
- else:
- self._resmap[i] = 0.0
-
- # Compute significance map
- self._resmap *= 2.0
- self._resmap = self._resmap.sqrt()
- self._resmap = self._resmap * signmap
-
- # Raise exception if algorithm is unknown
- else:
- raise TypeError('Algorithm "'+algorithm+'" not known')
+ self._resmap = obsutils.residuals(self,countmap1,modelmap)
# Optionally publish map
if self['publish'].boolean():
diff --git a/cscripts/csresspec.par b/cscripts/csresspec.par
new file mode 100644
index 0000000..0e8b281
--- /dev/null
+++ b/cscripts/csresspec.par
@@ -0,0 +1,64 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: csresspec.par #
+# #
+# Author: Luigi Tibaldo #
+# luigi.tibaldo@irap.omp.eu #
+# IRAP #
+# #
+# Purpose: Parameter file for the csresspec script #
+# #
+#############################################################################
+#
+# Input and output parameters
+#============================
+inobs, f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
+inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
+modcube, f, a, NONE,,, "Input model cube file (generated with ctmodel)"
+expcube, f, a, NONE,,, "Input exposure cube file"
+psfcube, f, a, NONE,,, "Input PSF cube file"
+edispcube, f, a, NONE,,, "Input energy dispersion cube file"
+bkgcube, f, a, NONE,,, "Input background cube file"
+caldb, s, a, prod2,,, "Calibration database"
+irf, s, a, South_0.5h,,, "Instrument response function"
+edisp, b, h, no,,, "Apply energy dispersion?"
+outfile, f, a, resspec.fits,,, "Output residual spectrum file"
+
+#
+# Script parameters
+#==================
+statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
+components,b, h, no,,, "Calculate model for individual components"
+ebinalg, s, h, LOG,LIN|LOG|FILE,, "Binning algorithm"
+emin, r, a, 0.1,,, "Lower energy limit (TeV)"
+emax, r, a, 100.0,,, "Upper energy limit (TeV)"
+enumbins, i, a, 20,,, "Number of model cube energy bins"
+ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
+stack, b, a, no,,, "Stack observations"
+coordsys, s, a, CEL,CEL|GAL,, "Coordinate System"
+proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
+xref, r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
+yref, r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
+nxpix, i, a, 200,,, "Size of the X axis in pixels"
+nypix, i, a, 200,,, "Size of the Y axis in pixels"
+binsz, r, a, 0.02,,, "Pixel size (deg/pixel)"
+mask, b, a, no,,, "Mask data to calculate residuals in ROI"
+ra, r, a, 83.63,0,360, "RA for ROI centre (degrees)"
+dec, r, a, 22.01,-90,90, "Dec for ROI centre (degrees)"
+rad, r, a, 3.0,0,180, "Radius of ROI (degrees)"
+regfile, f, a, "NONE",,, "Input exclusion region file in ds9 format"
+algorithm, s, a, SIGNIFICANCE,SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE,, "Residuals computation algorithm"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish residual map on VO Hub?"
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug, b, h, no,,, "Debugging mode activated"
+mode, s, h, ql,,, "Mode of automatic parameters"
+logfile, f, h, csresspec.log,,, "Log filename"
diff --git a/cscripts/csresspec.py b/cscripts/csresspec.py
new file mode 100755
index 0000000..0c8ac79
--- /dev/null
+++ b/cscripts/csresspec.py
@@ -0,0 +1,739 @@
+#! /usr/bin/env python
+# ==========================================================================
+# Generates a residual spectrum.
+#
+# Copyright (C) 2017-2018 Luigi Tibaldo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import gammalib
+import ctools
+from cscripts import obsutils
+import sys
+
+
+# =============== #
+# csresspec class #
+# =============== #
+class csresspec(ctools.csobservation):
+ """
+ Generates a residual spectrum
+ """
+
+ # Constructor
+ def __init__(self, *argv):
+ """
+ Constructor
+ """
+ # Initialise application by calling the appropriate class constructor
+ self._init_csobservation(self.__class__.__name__, ctools.__version__,
+ argv)
+
+ # Initialise class members
+ self._use_maps = False
+ self._stack = False
+ self._mask = False
+ self._fits = gammalib.GFits()
+
+ # Return
+ return
+
+ # Private methods
+ def _get_parameters(self):
+ """
+ Get parameters from parfile and setup the observation
+ """
+ # Setup observations (require response and allow event list as well as
+ # counts cube)
+ self._setup_observations(self.obs(), True, True, True)
+
+ # Set observation statistic
+ self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
+
+ # Collect number of unbinned, binned and On/Off observations in
+ # observation container
+ n_unbinned = 0
+ n_binned = 0
+ n_onoff = 0
+ for obs in self.obs():
+ if obs.classname() == 'GCTAObservation':
+ if obs.eventtype() == 'CountsCube':
+ n_binned += 1
+ else:
+ n_unbinned += 1
+ elif obs.classname() == 'GCTAOnOffObservation':
+ n_onoff += 1
+ n_cta = n_unbinned + n_binned + n_onoff
+ n_other = self.obs().size() - n_cta
+
+ # Query whether to compute model for individual components
+ components = self['components'].boolean()
+
+ # If there is only one binned observation and no model for individual
+ # components is required, query for precomputed model file and set
+ # use_maps to True
+ if self.obs().size() == 1 and n_binned == 1 and not components:
+ modcube = self['modcube'].filename()
+ if modcube != 'NONE':
+ self._use_maps = True
+
+ # If there are unbinned observations query the energy binning parameters
+ if n_unbinned != 0:
+ self['ebinalg'].string()
+ if self['ebinalg'].string() == 'FILE':
+ self['ebinfile'].filename()
+ else:
+ self['emin'].real()
+ self['emax'].real()
+ self['enumbins'].integer()
+ if n_cta > n_unbinned:
+ n_notunbin = n_cta - n_unbinned
+
+ # If there is more than one observation, and observations are all
+ # unbinned or all onoff query user to know if they wish stacked results
+ if self.obs().size() > 1 and \
+ (n_unbinned == self.obs().size() or n_onoff == self.obs().size()):
+ self._stack = self['stack'].boolean()
+ # If we are to stack event lists query parameters for cube creation
+ if self._stack and n_unbinned == self.obs().size():
+ self['coordsys'].string()
+ self['proj'].string()
+ self['xref'].real()
+ self['yref'].real()
+ self['nxpix'].integer()
+ self['nypix'].integer()
+ self['binsz'].real()
+
+ # If we are not using a precomputed model and no models are available
+ # in the observation container query input XML model file
+ if not self._use_maps and self.obs().models().size() == 0:
+ self.obs().models(self['inmodel'].filename())
+
+ # Unless all observations are On/Off query for mask definition
+ if n_onoff == n_cta:
+ pass
+ else:
+ self._mask = self['mask'].boolean()
+ if self._mask:
+ self['ra'].real()
+ self['dec'].real()
+ self['rad'].real()
+ self['regfile'].filename()
+
+ # Unless all observations are On/Off, or we are using precomputed model
+ # maps query whether to use energy dispersion
+ if n_onoff == n_cta or self._use_maps:
+ pass
+ else:
+ self['edisp'].boolean()
+
+ # Query algorithm for residual computation
+ self['algorithm'].string()
+
+ # Read ahead output parameters
+ if self._read_ahead():
+ self['outfile'].filename()
+
+ # Write input parameters into logger
+ self._log_parameters(gammalib.TERSE)
+
+ # Write header for observation census
+ self._log_header1(gammalib.TERSE, 'Observation census')
+
+ # Log census of input observations
+ self._log_value(gammalib.NORMAL, 'Unbinned CTA observations', n_unbinned)
+ self._log_value(gammalib.NORMAL, 'Binned CTA observations', n_binned)
+ self._log_value(gammalib.NORMAL, 'On/off CTA observations', n_onoff)
+ self._log_value(gammalib.NORMAL, 'Other observations', n_other)
+ if n_other > 0:
+ msg = 'WARNING: Only CTA observation can be handled, all non-CTA ' \
+ + 'observations will be ignored.'
+ self._log_string(gammalib.TERSE, msg)
+
+ # Log for unbinned observations
+ if n_unbinned != 0:
+ msg = ' User defined energy binning will be used for %d unbinned '\
+ 'observations.' % (n_unbinned)
+ self._log_string(gammalib.TERSE, msg)
+ if n_cta > n_unbinned:
+ msg = ' The intrinsic binning will be used for the remaining '\
+ '%d CTA observations.' % (n_notunbin)
+ self._log_string(gammalib.TERSE, msg)
+
+ # Signal how energy dispersion is applied
+ if n_onoff == n_cta or self._use_maps:
+ msg = ' Energy dispersion is applied based on the input data/model ' \
+ + 'and not according to the edisp parameter'
+ self._log_string(gammalib.TERSE, msg)
+
+ # Return
+ return
+
+ def _stack_observations(self):
+ """
+ Stack multiple observations and replace observation container
+ with single stacked observation
+
+ Returns
+ -------
+ obs : `~gammalib.GObservations'
+ Observation container with a single stacked observation
+ """
+ # Use first observation to determine the type and apply the
+ # appropriate stacking method
+
+ # If On/Off use the GCTAOnOffObservation constructor
+ if self.obs()[0].classname() == 'GCTAOnOffObservation':
+ msg = 'Stacking %d On/Off observations.' % (self.obs().size())
+ self._log_string(gammalib.NORMAL, msg)
+ stacked = gammalib.GCTAOnOffObservation(self.obs())
+ stacked_obs = gammalib.GObservations()
+ stacked_obs.append(stacked)
+ stacked_obs.models(self.obs().models())
+
+ # If event list then bin observations
+ elif self.obs()[0].classname() == 'GCTAObservation':
+ msg = 'Stacking %d event lists.' % (self.obs().size())
+ self._log_string(gammalib.NORMAL, msg)
+ stacked_obs = obsutils.get_stacked_obs(self, self.obs())
+
+ # Return
+ return stacked_obs
+
+ def _bin_evlist(self, obs):
+ """
+ Turn single event list into counts cube
+
+ Parameters
+ ----------
+ obs : `~gammalib.GObservations()`
+ Observation container with single event list
+
+ Returns
+ -------
+ obs : `~gammalib.GObservations()`
+ Binned observation container dict with event list ROI and energy
+ range information
+ """
+ # Retrieve information about ROI in event list
+ roi = obs[0].roi()
+ ra = roi.centre().dir().ra_deg()
+ dec = roi.centre().dir().dec_deg()
+ rad = roi.radius()
+
+ # We will cover the whole ROI with 0.02 deg binning
+ npix = int(2.0 * rad / 0.02) + 1
+
+ # Log binning of events
+ self._log_string(gammalib.EXPLICIT, 'Binning events')
+
+ # Bin events
+ cntcube = ctools.ctbin(obs)
+ cntcube['xref'] = ra
+ cntcube['yref'] = dec
+ cntcube['binsz'] = 0.02
+ cntcube['nxpix'] = npix
+ cntcube['nypix'] = npix
+ cntcube['proj'] = 'TAN'
+ cntcube['coordsys'] = 'CEL'
+ cntcube['ebinalg'] = self['ebinalg'].string()
+ if self['ebinalg'].string() == 'FILE':
+ cntcube['ebinfile'] = self['ebinfile'].filename()
+ else:
+ cntcube['enumbins'] = self['enumbins'].integer()
+ cntcube['emin'] = self['emin'].real()
+ cntcube['emax'] = self['emax'].real()
+ cntcube.run()
+
+ # Retrieve the binned observation container
+ binned_obs = cntcube.obs().copy()
+
+ # Log response computation
+ self._log_string(gammalib.EXPLICIT, 'Computing binned response')
+
+ # Compute binned response
+ response = obsutils.get_stacked_response(obs, ra, dec,
+ binsz=0.02,
+ nxpix=npix,
+ nypix=npix,
+ coordsys='CEL',
+ proj='TAN',
+ emin=self['emin'].real(),
+ emax=self['emax'].real(),
+ edisp=self['edisp'].boolean())
+
+ if self['edisp'].boolean():
+ binned_obs[0].response(response['expcube'], response['psfcube'],
+ response['edispcube'], response['bkgcube'])
+ else:
+ binned_obs[0].response(response['expcube'], response['psfcube'],
+ response['bkgcube'])
+
+ # Get new models
+ models = response['models']
+
+ # Set models for observation container
+ binned_obs.models(models)
+
+ # Check if energy boundaries provided by user extend beyond
+ # the content of the event list
+ if self['emin'].real() > obs[0].events().emin().TeV():
+ emin = 'INDEF'
+ else:
+ emin = obs[0].events().emin().TeV()
+ if self['emax'].real() < obs[0].events().emax().TeV():
+ emax = 'INDEF'
+ else:
+ emax = obs[0].events().emax().TeV()
+
+ # Put ROI and E bound info in dictionary
+ info = {'was_list': True, 'roi_ra': ra, 'roi_dec': dec, 'roi_rad': rad,
+ 'emin': emin, 'emax': emax}
+
+ # Return new oberservation container
+ return binned_obs, info
+
+ def _masked_cube(self, cube, ra, dec, rad, emin='INDEF', emax='INDEF',
+ regfile='NONE'):
+ """
+ Mask an event cube and returns the masked cube
+
+ Parameters
+ ----------
+ cube : `~gammalib.GCTAEventCube()`
+ Event cube
+ ra : float (str 'INDEF' for no selection on direction)
+ Right Ascension (deg)
+ dec : float (str 'INDEF' for no selection on direction)
+ Declination (deg)
+ rad : float (str 'INDEF' for no selection on direction)
+ Radius (deg)
+ emin : float (str 'INDEF' for no selection on energy)
+ Minimum energy (TeV)
+ emax : float (str 'INDEF' for no selection on energy)
+ Maximum energy (TeV)
+
+ Returns
+ -------
+ cube : `~gammalib.GCTAEventCube()`
+ Event cube
+ """
+ # Turn cube into observation container to feed to ctcubemask
+ obs = gammalib.GCTAObservation()
+ obs.events(cube)
+ obs_cont = gammalib.GObservations()
+ obs_cont.append(obs)
+
+ # Use ctcubemask to mask event cube pixels
+ cubemask = ctools.ctcubemask(obs_cont)
+ cubemask['ra'] = ra
+ cubemask['dec'] = dec
+ cubemask['rad'] = rad
+ cubemask['emin'] = emin
+ cubemask['emax'] = emax
+ cubemask['regfile'] = regfile
+ cubemask.run()
+
+ # Extract copy of cube from observation container (copy is needed to
+ # avoid memory leaks in SWIG)
+ cube = cubemask.obs()[0].events().copy()
+
+ # Return cube
+ return cube
+
+ def _cube_to_spectrum(self, cube, evlist_info):
+ """
+ Derive from event cube a count spectrum
+
+ If data come from event list use only the ROI and energy range of
+ the original data. Apply user defined mask if requested.
+
+ Parameters
+ ----------
+ cube : `~gammalib.GCTAEventCube()`
+ Event cube
+ evlist_info : dict
+ Dictionary with information on original event list
+
+ Returns
+ -------
+ array : `~gammalib.GNdarray()'
+ Counts spectrum
+ """
+ # If we started from event list mask the ROI only
+ # for residual computation
+ if evlist_info['was_list']:
+ msg = 'Masking ROI from original event list'
+ self._log_string(gammalib.EXPLICIT, msg)
+ cube = self._masked_cube(cube, evlist_info['roi_ra'],
+ evlist_info['roi_dec'],
+ evlist_info['roi_rad'],
+ emin=evlist_info['emin'],
+ emax=evlist_info['emax'])
+
+ # Apply user mask
+ if self._mask:
+ msg = 'Masking ROI requested by user'
+ self._log_string(gammalib.EXPLICIT, msg)
+ cube = self._masked_cube(cube, self['ra'], self['dec'],
+ self['rad'],
+ regfile=self['regfile'])
+
+ # Extract skymap and clip at 0 to null masked areas
+ counts = cube.counts().copy()
+ counts = counts.clip(0.)
+
+ # Convert skymap into GNdarray count spectrum
+ counts = counts.counts()
+
+ # Return
+ return counts
+
+ def _residuals_table(self, obs_id, ebounds, counts, model, residuals):
+ """
+ Create a Fits Table and store counts, model, and residuals
+
+ Parameters
+ ----------
+ obs_id : str
+ Observation id
+ ebounds : `~gammalib.GEbounds()'
+ Energy boundaries
+ counts : `~gammalib.GNdarray()'
+ Counts spectrum
+ model : `~gammalib.GNdarray()'
+ Model spectrum
+ residuals : `~gammalib.GNdarray()'
+ Residual spectrum
+
+ Returns
+ -------
+ table : `~gammalib.GFitsBinTable()'
+ Residual spectrum as FITS binary table
+ """
+ # Create FITS table columns
+ nrows = ebounds.size()
+ energy_low = gammalib.GFitsTableDoubleCol('Emin', nrows)
+ energy_high = gammalib.GFitsTableDoubleCol('Emax', nrows)
+ counts_col = gammalib.GFitsTableDoubleCol('Counts', nrows)
+ model_col = gammalib.GFitsTableDoubleCol('Model', nrows)
+ resid_col = gammalib.GFitsTableDoubleCol('Residuals', nrows)
+ energy_low.unit('TeV')
+ energy_high.unit('TeV')
+
+ # Fill FITS table columns
+ for i in range(nrows):
+ energy_low[i] = ebounds.emin(i).TeV()
+ energy_high[i] = ebounds.emax(i).TeV()
+ counts_col[i] = counts[i]
+ model_col[i] = model[i]
+ resid_col[i] = residuals[i]
+
+ # Initialise FITS Table with extension set to obs id
+ table = gammalib.GFitsBinTable(nrows)
+ # If name not empty add leading blank
+ if obs_id != '':
+ obs_id = ' ' + obs_id
+ table.extname('RESIDUALS' + obs_id)
+
+ # Add Header card to specify algorithm used for residual computation
+ table.card('ALGORITHM', self['algorithm'].string(),
+ 'Algorithm used for computation of residuals')
+
+ # Append filled columns to fits table
+ table.append(energy_low)
+ table.append(energy_high)
+ table.append(counts_col)
+ table.append(model_col)
+ table.append(resid_col)
+
+ # Return binary table
+ return table
+
+ def _append_column(self, table, name, data):
+ """
+ Append optional column to residual table
+
+ Parameters
+ ----------
+ table : `~gammalib.GFitsBinTable()'
+ FITS binary table
+ name : str
+ Column name
+ data : `~gammalib.GEbounds()'
+ Data to be filled into new column
+
+ Returns
+ -------
+ table : `~gammalib.GFitsBinTable()'
+ FITS binary table
+ """
+ # Check size compatibility
+ if table.nrows() == data.size():
+ pass
+ # Otherwise throw error
+ else:
+ msg = 'csresspec._append_column: FITS table and data have ' \
+ 'incompatible size.'
+ raise RuntimeError(msg)
+
+ # Create column
+ column = gammalib.GFitsTableDoubleCol(name, table.nrows())
+
+ # Fill data
+ for i, value in enumerate(data):
+ column[i] = value
+
+ # Append new column to table
+ table.append(column)
+
+ # Return modified table
+ return table
+
+
+ # 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(), 'Observation')
+
+ # Stack observations if requested
+ if self._stack:
+ self._log_header1(gammalib.TERSE, 'Stacking Observations')
+ self.obs(self._stack_observations())
+
+ # Log processing header
+ self._log_header1(gammalib.TERSE, 'Processing Observations')
+
+ # Loop over observations and calculate residuals
+ for s, observation in enumerate(self.obs()):
+
+ # Retrieve and store obs id
+ obs_id = observation.id()
+
+ # If observation id is empty and there is more than one observation
+ # replace with incremental number
+ if obs_id == '' and self.obs().size() > 1:
+ obs_id = str(s)
+
+ # Log processing of observation
+ if self.obs().size() > 1:
+ self._log_header2(gammalib.NORMAL,
+ 'Processing observation %s' % obs_id)
+
+ # Turn into observation container and assign models
+ obs = gammalib.GObservations()
+ obs.append(observation)
+ obs.models(self.obs().models())
+
+ # if 3D observation
+ if obs[0].classname() == 'GCTAObservation':
+
+ # Prepare Observations
+
+ # If already binned set the evlist_info dictionary to have
+ # attribute was_list False
+ if obs[0].eventtype() == 'CountsCube':
+ evlist_info = {'was_list': False}
+ # Otherwise bin now
+ else:
+ # we remember if we binned an event list
+ # so that we can mask only the ROI for residual calculation
+ msg = 'Setting up binned observation'
+ self._log_string(gammalib.NORMAL, msg)
+ obs, evlist_info = self._bin_evlist(obs)
+
+ # Calculate Model and residuals
+ # If model cube is provided load it
+ if self._use_maps:
+ modcube = gammalib.GCTAEventCube(self['inmodel'].filename())
+ # Otherwise calculate it now
+ else:
+ msg = 'Computing model cube'
+ self._log_string(gammalib.NORMAL, msg)
+ modelcube = ctools.ctmodel(obs)
+ modelcube.run()
+ modcube = modelcube.cube().copy()
+
+ # Extract cntcube for residual computation
+ cntcube = obs[0].events().copy()
+
+ # Derive count spectra from cubes
+ msg = 'Computing counts, model, and residual spectra'
+ self._log_string(gammalib.NORMAL, msg)
+ counts = self._cube_to_spectrum(cntcube, evlist_info)
+ model = self._cube_to_spectrum(modcube, evlist_info)
+
+ # Calculate residuals
+ residuals = obsutils.residuals(self, counts, model)
+
+ # Extract energy bounds
+ ebounds = cntcube.ebounds()
+
+ # Fill results table
+ msg = 'Filling residual table'
+ self._log_string(gammalib.NORMAL, msg)
+ table = self._residuals_table(obs_id, ebounds, counts, model,
+ residuals)
+
+ # Calculate models of individual components if requested
+ if self['components'].boolean():
+ for component in obs.models():
+ self._log_value(gammalib.NORMAL,
+ 'Computing model for component',
+ component.name())
+ # Set model cube models to
+ # individual component
+ model_cont = gammalib.GModels()
+ model_cont.append(component)
+ modelcube.obs().models(model_cont)
+
+ # Run model cube
+ modelcube.run()
+
+ # Extract spectrum of individual component
+ modcube = modelcube.cube().copy()
+ model = self._cube_to_spectrum(modcube, evlist_info)
+
+ # append component to table
+ table = self._append_column(table, component.name(),
+ model)
+
+ # otherwise, if On/Off
+ elif obs[0].classname() == 'GCTAOnOffObservation':
+
+ onoff = obs[0]
+
+ # Calculate Counts, Model and residuals
+ msg = 'Computing counts, model, and residual spectra'
+ self._log_string(gammalib.NORMAL, msg)
+
+ # On spectrum
+ counts = onoff.on_spec().counts_spectrum()
+
+ # Model for On region
+ background = onoff.model_background(obs.models()).counts_spectrum()
+ alpha = onoff.on_spec().backscal_spectrum()
+ model = background.copy()
+ model *= alpha
+ model += onoff.model_gamma(obs.models()).counts_spectrum()
+
+ # On Residuals
+ residuals = obsutils.residuals(self, counts, model)
+
+ # Extract energy bounds
+ ebounds = onoff.on_spec().ebounds()
+
+ # Fill results table
+ msg = 'Filling residual table'
+ self._log_string(gammalib.NORMAL, msg)
+ table = self._residuals_table(obs_id, ebounds, counts, model,
+ residuals)
+
+ # Get Off spectrum and add to table
+ msg = 'Computing counts, model, and residual spectra for Off regions'
+ self._log_string(gammalib.NORMAL, msg)
+ counts_off = onoff.off_spec().counts_spectrum()
+ table = self._append_column(table, 'Counts_Off',
+ counts_off)
+
+ # Add background/Off model to table
+ table = self._append_column(table, 'Model_Off',
+ background)
+
+ # Calculate Off residuals and add to table
+ residuals_off = obsutils.residuals(self, counts_off, background)
+ table = self._append_column(table, 'Residuals_Off',
+ residuals_off)
+
+ # Calculate models of individual components if requested
+ if self['components'].boolean():
+ for component in obs.models():
+ # If background pass
+ # We always add the background at the end so that
+ # we accommodate WSTAT for which the background is not
+ # mandatory in the model
+ if component.classname() == 'GCTAModelIrfBackground':
+ pass
+ # Otherwise calculate gamma component and append to Table
+ else:
+ self._log_value(gammalib.NORMAL,
+ 'Computing model for component',
+ component.name())
+ # Create observation container for individual components
+ model_cont = gammalib.GModels()
+ model_cont.append(component)
+ # Calculate gamma model
+ model = onoff.model_gamma(model_cont)
+ model = model.counts_spectrum()
+ # Append to table
+ table = self._append_column(table, component.name(),
+ model)
+ # Add background already calculated
+ self._log_value(gammalib.NORMAL,
+ 'Computing model for component',
+ 'Background')
+ background *= onoff.on_spec().backscal_spectrum()
+ table = self._append_column(table, 'Background',
+ background)
+
+ # Append results table to output file
+ self._fits.append(table)
+
+ # Return
+ return
+
+ def save(self):
+ """
+ Save residuals
+ """
+ # Write header
+ self._log_header1(gammalib.TERSE, 'Save residuals')
+
+ # Continue only if FITS file is valid
+ if self._fits != None:
+
+ # Get outfile parameter
+ outfile = self['outfile'].filename()
+
+ # Log file name
+ self._log_value(gammalib.NORMAL, 'Residuals file', outfile.url())
+
+ # Save residuals
+ self._fits.saveto(outfile, self['clobber'].boolean())
+
+ # Return
+ return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Create instance of application
+ app = csresspec(sys.argv)
+
+ # Execute application
+ app.execute()
diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py
index f7018a2..5138794 100755
--- a/cscripts/obsutils.py
+++ b/cscripts/obsutils.py
@@ -1,7 +1,7 @@
# ==========================================================================
# Utility functions for observation handling
#
-# Copyright (C) 2011-2017 Juergen Knoedlseder
+# Copyright (C) 2011-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
@@ -829,3 +829,102 @@ def get_onoff_obs(cls, obs):
# Return On/Off oberservation container
return onoff_obs
+
+
+# ========================================= #
+# Calculate residuals from counts and model #
+# ========================================= #
+def residuals(cls, counts, model):
+ """
+ Calculate residuals given counts and models, according to algorithm
+ specified by user.
+
+ Can handle GSkyMap or GNdarray objects
+
+
+ Parameters
+ ----------
+ cls : `~ctools.cscript`
+ cscript class
+ counts : `~gammalib.GSkyMap/~gammalib.GNdarray'
+ Data counts
+ model : `~gammalib.GSkyMap/~gammalib.GNdarray'
+ Model counts
+
+ Returns
+ -------
+ residuals : `~gammalib.GSkyMap/~gammalib.GNdarray'
+ Residuals
+ """
+ # Find type of objects we are manipulating and set size to iterate later
+
+ # If GNdarray
+ if counts.classname() == 'GNdarray':
+ nelem = counts.size()
+ elif counts.classname() == 'GSkyMap':
+ nelem = counts.npix()
+ else:
+ msg = 'cscripts.obsutils.residuals only handles ' \
+ + 'gammalib.GNdarray or gammalib.GSkyMap objects.\n'
+ raise RuntimeError(msg)
+
+ # Copy counts to initialise residuals
+ residuals = counts.copy()
+
+ # Get residual map algorithm type
+ algorithm = cls['algorithm'].string()
+
+ # Subtract
+ if algorithm == 'SUB':
+ residuals -= model
+
+ # Subtract and divide by model
+ elif algorithm == 'SUBDIV':
+ residuals -= model
+ residuals /= model
+
+ # Subtract and divide by sqrt of model
+ elif algorithm == 'SUBDIVSQRT':
+ residuals -= model
+ residuals /= model.sqrt()
+
+ # Calculate significance from Li&Ma
+ elif algorithm == 'SIGNIFICANCE':
+
+ # Compute sign
+ sign = (residuals - model).sign()
+
+ # Loop over every bin
+ for i in range(nelem):
+
+ # If the model value > 0.0 do the computation as normal ...
+ model_val = model[i]
+ if model_val > 0.0:
+
+ # If the data value is also > 0 then compute the
+ # significance^2 and save it ...
+ data_val = residuals[i]
+ if data_val > 0.0:
+ log_val = math.log(data_val / model_val)
+ residuals[i] = (data_val * log_val) + model_val - data_val
+
+ # ... otherwise compute the reduced value of the above
+ # expression. This is necessary to avoid computing log(0).
+ else:
+ residuals[i] = model_val
+
+ # ... otherwise hard-code the significance to 0
+ else:
+ residuals[i] = 0.0
+
+ # Compute significance map
+ residuals *= 2.0
+ residuals = residuals.sqrt()
+ residuals *= sign
+
+ # Raise exception if algorithm is unknown
+ else:
+ raise TypeError('Algorithm "' + algorithm + '" not known')
+
+ # Return
+ return residuals
diff --git a/doc/source/users/reference_manual/csresspec.rst b/doc/source/users/reference_manual/csresspec.rst
new file mode 100644
index 0000000..8061786
--- /dev/null
+++ b/doc/source/users/reference_manual/csresspec.rst
@@ -0,0 +1,189 @@
+.. _csresspec:
+
+csresspec
+=========
+
+Generates residual spectrum.
+
+
+Synopsis
+--------
+
+This scripts generates a residual spectrum for a given model. It works for
+event lists, counts cubes, On/Off observations or observation definition files.
+For event lists parameters that define the spectral binning need to
+be provided so that the script can bin the data internally. The model
+is then convolved with the instrumental response function for the chosen
+binning (either specified by the user or intrinsic to the observations)
+and used for residual computation. It is possible to compute the model
+contributions from individual model components.
+
+If starting from event lists or count cubes,
+before residual computation the 3D counts and model cubes are collapsed into
+spectra by summing over all spatial direction. It is possible to specify a mask
+to compute counts, model, and residuals only for a region of interest.
+
+For On/Off observations, counts, model, and residual spectra for the Off regions
+are computed as well.
+
+For an observation definition file, when all observations are event lists or
+On/Off observations it is possible to stack them before residual computation. To
+do so spatial binning parameters for event lists are required. Otherwise the
+residuals are calculated separately for each observation.
+
+Four options exist then for residual computation:
+
+* the subtraction of the model from the counts (SUB)
+* the subtraction and division by the model (SUBDIV)
+* the subtraction and division by the square root of the model (SUBDIVSQRT)
+* the computation of the significance using Li & Ma (SIGNIFICANCE)
+
+The counts, model, and residual spectra are written into a FITS file.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+ Input event list, counts cube or observation definition XML file.
+
+``modcube [file]``
+ Input model cube file (generated with ctmodel).
+
+``inmodel [file]``
+ Input model XML file.
+
+``expcube [file]``
+ Input exposure cube file.
+
+``psfcube [file]``
+ Input PSF cube file.
+
+``edispcube [file]``
+ Input energy dispersion cube file.
+
+``bkgcube [file]``
+ Input background cube file.
+
+``caldb [string]``
+ Calibration database.
+
+``irf [string]``
+ Instrument response function.
+
+``(edisp = no) [boolean]``
+ Apply energy dispersion?
+
+``outfile [file]``
+ Output file.
+
+``(statistic = DEFAULT) [string]``
+ Optimization statistic. For On/Off observations determines how
+ the background counts are computed. ``DEFAULT``, i.e., ``CSTAT``,
+ will use the input background model. ``WSTAT`` uses the background
+ counts from likelihood profiling as defined for the W statistic.
+
+``ebinalg [string]``
+ Algorithm for defining energy bins.
+
+``emin [real]``
+ Lower energy value for first energy bin (in TeV).
+
+``emax [real]``
+ Upper energy value for last energy bin (in TeV).
+
+``enumbins [integer]``
+ Number of energy bins.
+
+``ebinfile [file]``
+ Name of the file containing the energy bin definition.
+
+``(stack = no) [boolean]``
+ Stack multiple observations for residual calculation.
+
+``coordsys [string]``
+ Coordinate system (CEL - celestial, GAL - galactic) for stacking
+ multiple event lists.
+
+``proj [string]``
+ Projection method for stacking multiple event lists.
+
+``xref [real]``
+ Right Ascension / Galactic longitude of cube centre (J2000, in degrees)
+ for stacking multiple event lists.
+
+``yref [real]``
+ Declination / Galactic latitude of cube centre (J2000, in degrees) for stacking
+ multiple event lists.
+
+``nxpix [integer]``
+ Number of cube bins in Right Ascension or Galactic longitude for stacking
+ multiple event lists.
+
+``nypix [integer]``
+ Number of cube bins in Declination or Galactic latitude for stacking
+ multiple event lists.
+
+``binsz [real]``
+ Cube bin size (in degrees/pixel) for stacking
+ multiple event lists.
+
+``ra [real]``
+ Right Ascension of circular selection region centre (J2000, in degrees).
+
+``dec [real]``
+ Declination of circular selection region centre (J2000, in degrees).
+
+``rad [real]``
+ Radius of circular selection region (in degrees).
+
+``regfile [file]``
+ Input exclusion region file in ds9 format.
+
+``algorithm [string]``
+ Algorithm used to generate the residual map:
+
+ ``SUB``: :math:`DATA - MODEL`
+
+ ``SUBDIV``: :math:`(DATA - MODEL) / MODEL`
+
+ ``SUBDIVSQRT``: :math:`(DATA - MODEL) / \sqrt{MODEL}`
+
+ ``SIGNIFICANCE``: :math:`{\rm sign}(DATA-MODEL) \times \sqrt{ 2 \times ( DATA \times \ln \left(\frac{DATA}{MODEL} \right) + MODEL - DATA ) }`
+
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+ Specifies whether the residual map should be published on VO Hub.
+
+``(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 residual map 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 = csresmap.log) [string]``
+ Log filename.
+
+
+Related tools or scripts
+------------------------
+
+None
diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst
index aab4ad0..8cf998e 100644
--- a/doc/source/users/reference_manual/index.rst
+++ b/doc/source/users/reference_manual/index.rst
@@ -57,6 +57,7 @@ cscripts
csphasecrv --- Computes phase curve
cspull --- Generates pull distribution
csresmap --- Generates residual map
+ csresspec --- Generates residual spectrum
csroot2caldb --- Creates a caldb entry from a ROOT file
cssens --- Computes CTA sensitivity
csspec --- Computes spectral points
diff --git a/doc/source/users/tutorials/beyond/beyond_residual_spectrum.rst b/doc/source/users/tutorials/beyond/beyond_residual_spectrum.rst
new file mode 100644
index 0000000..851681b
--- /dev/null
+++ b/doc/source/users/tutorials/beyond/beyond_residual_spectrum.rst
@@ -0,0 +1,72 @@
+.. _sec_residual_spectrum:
+
+Generating a residual spectrum
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In the previous tutorial you have used the :ref:`csresmap` script to
+inspect the spatial residuals. You can also inspect the spectral
+residuals using the :ref:`csresspec` script.
+
+The :ref:`csresspec` script can deal with any observation type (unbinned,
+binned, stacked, On/Off). For the example case with unbinned
+observations the script is then used as follows:
+
+.. code-block:: bash
+
+ $csresspec components=yes
+ Input event list, counts cube, or observation definition XML file [obs.xml]
+ Calibration database [prod2]
+ Instrument response function [South_0.5h]
+ Lower energy limit (TeV) [0.1]
+ Upper energy limit (TeV) [100.0]
+ Number of model cube energy bins [30]
+ Stack observations [yes]
+ Coordinate System (CEL|GAL) [CEL]
+ Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR]
+ First coordinate of image/source region center in degrees (RA or galactic l) (0-360) [83.63]
+ Second coordinate of image/source region center in degrees (DEC or galactic b) (-90-90) [22.01]
+ Size of the X axis in pixels [200] 51
+ Size of the Y axis in pixels [200] 51
+ Pixel size (deg/pixel) [0.02]
+ Input model definition XML file [crab_results.xml]
+ Mask data to calculate residuals in ROI [yes] no
+ Residuals computation algorithm (SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE) [SIGNIFICANCE]
+ Output residual spectrum file [residual_spectrum.fits]
+
+In our example, the residuals will cover energies from 100 GeV to 100
+TeV, in 30 energy bins. Since we have decided to stack multiple
+unbinned observations we have also entered parameters to perform a
+spatial binning. Note that we have defined the spatial grid such that
+we will obtain the residuals in a 1 degree square region around the
+Crab. Optionally, a mask can be applied to either unbinned or binned
+3D observations. We have used the hidden ``components`` parameter to
+include the individual model components in the results file.
+
+The :ref:`csresspec` script created the FITS file ``residual_spectrum.fits`` that contains
+the results of the residual computation. A graphical display of the fit results is
+shown below.
+
+.. figure:: residual_spectrum.png
+ :width: 600px
+ :align: center
+
+ *Data, model, and residuals from the spectral fit of the Crab nebula*
+
+..
+
+
+.. note::
+
+ The figure was created by typing:
+
+ .. code-block:: bash
+
+ $ $CTOOLS/share/examples/python/show_spectrum.py residual_spectrum.fits
+
+ Use the option ``-h`` to select the table to display if there are
+ multiple observations.
+
+
+
+
+
diff --git a/doc/source/users/tutorials/beyond/index.rst b/doc/source/users/tutorials/beyond/index.rst
index 4aa4ef7..bb5cb77 100644
--- a/doc/source/users/tutorials/beyond/index.rst
+++ b/doc/source/users/tutorials/beyond/index.rst
@@ -17,6 +17,7 @@ analyses.
stacked-analysis
beyond_tsmap
beyond_residual
+ beyond_residual_spectrum
advanced-python-fitting
beyond_irf
beyond_model
diff --git a/doc/source/users/tutorials/beyond/residual_spectrum.png b/doc/source/users/tutorials/beyond/residual_spectrum.png
new file mode 100644
index 0000000..f63f26b
--- /dev/null
+++ b/doc/source/users/tutorials/beyond/residual_spectrum.png
diff --git a/examples/Makefile.am b/examples/Makefile.am
index aebc2c3..7ec9cc2 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -1,7 +1,7 @@
#############################################################################
# Makefile for ctools examples #
# ------------------------------------------------------------------------- #
-# Copyright (C) 2016-2017 Juergen Knoedlseder #
+# Copyright (C) 2016-2018 Juergen Knoedlseder #
# ------------------------------------------------------------------------- #
# #
# This program is free software: you can redistribute it and/or modify #
@@ -28,6 +28,7 @@ scripts = README.md \
show_pha.py \
show_pull_evolution.py \
show_pull_histogram.py \
+ show_residuals.py \
show_response.py \
show_sensitivity.py \
show_spectrum.py \
diff --git a/examples/show_residuals.py b/examples/show_residuals.py
new file mode 100755
index 0000000..d4d3d93
--- /dev/null
+++ b/examples/show_residuals.py
@@ -0,0 +1,315 @@
+#! /usr/bin/env python
+# ==========================================================================
+# Display residuals generated by csresspec
+#
+# Copyright (C) 2017-2018 Luigi Tibaldo
+#
+# 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 cscripts
+try:
+ import matplotlib.pyplot as plt
+ from matplotlib import gridspec
+ plt.figure()
+ plt.close()
+except (ImportError, RuntimeError):
+ print('This script needs the "matplotlib" module')
+ sys.exit()
+
+
+# ================================= #
+# Fill counts/model/residual arrays #
+# ================================= #
+def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
+ c_resid, algorithm):
+ """
+ Helper function that fills in the counts, model, residuals lists and
+ calculate relative errors for given table row
+
+ Parameters
+ ----------
+ row : int
+ Table row index
+ counts : list of float
+ Counts
+ model : list of float
+ Model
+ resid : list of float
+ Residual
+ e_counts : list of float
+ Counts error
+ e_resid : list of float
+ Residual error
+ c_counts : `~gammalib.GFitsTableCol()'
+ Counts FITS table column
+ c_model : `~gammalib.GFitsTableCol()'
+ Model FITS table column
+ c_resid : `~gammalib.GFitsTableCol()'
+ Residual FITS table column
+ algorithm : str
+ Algorithm
+
+ Returns
+ -------
+ counts, model, resid, e_counts, e_resid : tuple of float
+ Values
+ """
+ # Extract values
+ counts.append(c_counts.real(row))
+ model.append(c_model.real(row))
+ resid.append(c_resid.real(row))
+
+ # Calculate count error
+ err = math.sqrt(c_counts.real(row))
+ if err == 1:
+ # this prevents visualization problem in matplotlib
+ err = 0.99
+ else:
+ pass
+ e_counts.append(err)
+
+ # Calculate residual error
+ if algorithm == 'SUB':
+ e_resid.append(err)
+ elif algorithm == 'SUBDIV':
+ e_resid.append(err / c_model.real(row))
+ elif algorithm == 'SUBDIVSQRT':
+ e_resid.append(err / math.sqrt(c_model.real(row)))
+ elif algorithm == 'SIGNIFICANCE':
+ e_resid.append(1.)
+
+ # Return tuple
+ return counts, model, resid, e_counts, e_resid
+
+
+# ============= #
+# Plot spectrum #
+# ============= #
+def plot_residuals(filename, plotfile, hdu):
+ """
+ Plot spectrum
+
+ Parameters
+ ----------
+ filename : str
+ Name of spectrum FITS file
+ plotfile : str
+ Plot file name
+ hdu : int
+ Number of observation to plot
+ """
+ # Read spectrum file
+ fits = gammalib.GFits(filename)
+ table = fits.table(1 + hdu)
+ c_emin = table['Emin']
+ c_emax = table['Emax']
+ c_counts = table['Counts']
+ c_model = table['Model']
+ c_resid = table['Residuals']
+ try:
+ c_counts_off = table['Counts_Off']
+ c_model_off = table['Model_Off']
+ c_resid_off = table['Residuals_Off']
+ is_onoff = True
+ except:
+ is_onoff = False
+
+ # Initialise arrays to be filled
+ ebounds = []
+ ed_engs = []
+ eu_engs = []
+ em_engs = []
+ counts = []
+ e_counts = []
+ model = []
+ resid = []
+ e_resid = []
+ counts_off = []
+ e_counts_off = []
+ model_off = []
+ resid_off = []
+ e_resid_off = []
+
+ # Residual algorithm
+ algorithm = table.header()['ALGORITHM'].string()
+
+ # Loop over rows of the file
+ nrows = table.nrows()
+
+ # Add first energy boundary
+ ebounds.append(c_emin.real(0))
+ for row in range(nrows):
+
+ # Boundaries
+ ebounds.append(c_emax.real(row))
+
+ # Geometrical mean energy and errors
+ emean = math.sqrt(c_emin.real(row) * c_emax.real(row))
+ em_engs.append(emean)
+ ed_engs.append(emean - c_emin.real(row))
+ eu_engs.append(c_emax.real(row) - emean)
+
+ # counts, model, residuals
+ counts, model, resid, e_counts, e_resid = fill_cmr(row, counts, model,
+ resid, e_counts,
+ e_resid, c_counts,
+ c_model, c_resid,
+ algorithm)
+ # On/off
+ if is_onoff:
+ counts_off, model_off, resid_off, e_counts_off, e_resid_off = fill_cmr(
+ row, counts_off,
+ model_off, resid_off,
+ e_counts_off,
+ e_resid_off,
+ c_counts_off,
+ c_model_off, c_resid_off,
+ algorithm)
+
+ # Add model value to be compatible with plt.step
+ model = [model[0]] + model
+ if is_onoff:
+ model_off = [model_off[0]] + model_off
+
+ # Plot
+ axarr = []
+ if is_onoff:
+ f = plt.figure(figsize=(8, 4))
+ gs = gridspec.GridSpec(2, 2, height_ratios=[2, 1])
+ gs.update(hspace=0, left=0.09, right=0.97)
+ axarr.append(f.add_subplot(gs[0]))
+ axarr.append(f.add_subplot(gs[2]))
+ axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
+ axarr.append(f.add_subplot(gs[3], sharex=axarr[1]))
+ else:
+ f = plt.figure()
+ gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
+ gs.update(hspace=0)
+ axarr.append(f.add_subplot(gs[0]))
+ axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
+ plt.setp(axarr[1].get_xticklabels(), visible=False)
+ axarr[0].set_xscale('log')
+ axarr[1].set_xscale('log')
+ axarr[0].set_yscale('log')
+
+ # Counts and model
+ axarr[0].errorbar(em_engs, counts, yerr=e_counts, xerr=[ed_engs, eu_engs],
+ fmt='ko', capsize=0, linewidth=2, zorder=2, label='Data')
+ axarr[0].step(ebounds, model, color='0.5', linewidth=2, zorder=1,
+ label='Model')
+ axarr[0].set_ylabel('Counts')
+
+ # Residuals
+ axarr[1].errorbar(em_engs, resid, yerr=e_resid, xerr=[ed_engs, eu_engs],
+ fmt='ko', capsize=0, linewidth=2, zorder=2)
+ axarr[1].axhline(0, color='0.5', linestyle='--')
+ axarr[1].set_xlabel('Energy (TeV)')
+ if algorithm == 'SUB':
+ axarr[1].set_ylabel('Residuals (counts)')
+ elif algorithm == 'SUBDIV':
+ axarr[1].set_ylabel('Residuals (fraction)')
+ elif algorithm == 'SUBDIVSQRT' or algorithm == 'SIGNIFICANCE':
+ axarr[1].set_ylabel(r'Residuals ($\sigma$)')
+
+ # On/Off
+ if is_onoff:
+ axarr[0].set_title('ON')
+ axarr[2].set_title('OFF')
+ plt.setp(axarr[3].get_xticklabels(), visible=False)
+ axarr[2].set_xscale('log')
+ axarr[3].set_xscale('log')
+ axarr[2].set_yscale('log')
+
+ # Counts and model
+ axarr[2].errorbar(em_engs, counts_off, yerr=e_counts_off,
+ xerr=[ed_engs, eu_engs],
+ fmt='ko', capsize=0, linewidth=2, zorder=2,
+ label='Data')
+ axarr[2].step(ebounds, model_off, color='0.5', linewidth=2, zorder=1,
+ label='Model')
+
+ # Residuals
+ axarr[3].errorbar(em_engs, resid_off, yerr=e_resid_off,
+ xerr=[ed_engs, eu_engs],
+ fmt='ko', capsize=0, linewidth=2, zorder=2)
+ axarr[3].axhline(0, color='0.5', linestyle='--')
+ axarr[3].set_xlabel('Energy (TeV)')
+
+ # Add spectra of individual components
+ skiplist = ['Counts', 'Model', 'Residuals',
+ 'Counts_Off', 'Model_Off', 'Residuals_Off',
+ 'Emin', 'Emax']
+ for s in range(table.ncols()):
+ if table[s].name() in skiplist:
+ pass
+ else:
+ component = []
+ for row in range(nrows):
+ component.append(table[s].real(row))
+ component = [component[0]] + component
+ axarr[0].step(ebounds, component, zorder=0,
+ label=table[s].name())
+
+ # Add legend
+ axarr[0].legend(loc='best')
+
+ # Show figure
+ if len(plotfile) > 0:
+ plt.savefig(plotfile)
+ else:
+ plt.show()
+
+ # Return
+ return
+
+
+# =============== #
+# Show residuals #
+# =============== #
+def show_residuals():
+ """
+ Show residuals
+ """
+ # Set usage string
+ usage = 'show_spectrum.py [-p plotfile] [-h hdu] [file]'
+
+ # Set default options
+ options = [{'option': '-p', 'value': ''},
+ {'option': '-h', 'value': 0}]
+
+ # Get arguments and options from command line arguments
+ args, options = cscripts.ioutils.get_args_options(options, usage)
+
+ # Extract script parameters from options
+ plotfile = options[0]['value']
+ hdu = options[1]['value']
+
+ # Show residuals
+ plot_residuals(args[0], plotfile, hdu)
+
+ # Return
+ return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Show residuals
+ show_residuals()
diff --git a/src/ctcubemask/ctcubemask.hpp b/src/ctcubemask/ctcubemask.hpp
index 48ebbd0..463ffc6 100644
--- a/src/ctcubemask/ctcubemask.hpp
+++ b/src/ctcubemask/ctcubemask.hpp
@@ -33,6 +33,7 @@
#include "GammaLib.hpp"
#include "GCTALib.hpp"
#include "ctobservation.hpp"
+#include "GCTAEventCube.hpp"
/* __Definitions _________________________________________________________ */
#define CTCUBEMASK_NAME "ctcubemask"
@@ -87,5 +88,4 @@ protected:
bool m_select_energy; //!< Perform energy selection
};
-
#endif /* CTCUBEMASK_HPP */
diff --git a/test/test_csresspec.py b/test/test_csresspec.py
new file mode 100755
index 0000000..88f1957
--- /dev/null
+++ b/test/test_csresspec.py
@@ -0,0 +1,277 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the csresspec script.
+#
+# Copyright (C) 2017-2018 Luigi Tibaldo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import gammalib
+import cscripts
+from testing import test
+
+
+# =============================== #
+# Test class for csresspec script #
+# =============================== #
+class Test(test):
+ """
+ Test class for csresspec script
+
+ This test class makes unit tests for the csresspec script by using it
+ from the command line and from Python.
+ """
+
+ # Constructor
+ def __init__(self):
+ """
+ Constructor
+ """
+ # Call base class constructor
+ test.__init__(self)
+
+ # Set members
+ self._inonoff = self._datadir + '/onoff_obs.xml'
+ self._onoff_model = self._datadir + '/onoff_model.xml'
+
+ # Return
+ return
+
+ # Set test functions
+ def set(self):
+ """
+ Set all test functions
+ """
+ # Set test name
+ self.name('csresspec')
+
+ # Append tests
+ self.append(self._test_cmd, 'Test csresspec on command line')
+ self.append(self._test_python, 'Test csresspec from Python')
+
+ # Return
+ return
+
+ # Test csresspec on command line
+ def _test_cmd(self):
+ """
+ Test csresspec on the command line.
+ """
+ # Set script name
+ csresspec = self._script('csresspec')
+
+ # Setup csresspec command
+ cmd = csresspec + ' inobs="' + self._events + '"' + \
+ ' outfile="csresspec_cmd1.fits"' + \
+ ' inmodel="' + self._model + '"' + \
+ ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
+ ' emin=0.1 emax=100.0 enumbins=20' + \
+ ' mask=no algorithm="SUB"' + \
+ ' logfile="csresspec_cmd1.log" chatter=4'
+
+ # Check if execution of wrong command fails
+ self.test_assert(self._execute('command_that_does_not_exist') != 0,
+ 'Self test of test script')
+
+ # Check if execution was successful
+ self.test_assert(self._execute(cmd) == 0,
+ 'Check successful execution from command line')
+
+ # Check result file
+ self._check_result_file('csresspec_cmd1.fits', 20, 5)
+
+ # Setup csresspec command
+ cmd = csresspec + ' inobs="event_file_that_does_not_exist.fits"' + \
+ ' outfile="csresspec_cmd2.fits"' + \
+ ' inmodel="' + self._model + '"' + \
+ ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
+ ' emin=0.1 emax=100.0 enumbins=20' + \
+ ' mask=no algorithm="SUBDIVSQRT"' + \
+ ' logfile="csresspec_cmd2.log" debug=yes chatter=2'
+
+ # Check if execution failed
+ self.test_assert(self._execute(cmd) != 0,
+ 'Check invalid input file when executed from command line')
+
+ # Check csresspec --help
+ self._check_help(csresspec)
+
+ # Return
+ return
+
+ # Test csresspec from Python
+ def _test_python(self):
+ """
+ Test csresspec from Python
+ """
+ # Set-up csresspec for event list
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._events
+ resspec['outfile'] = 'csresspec_py1.fits'
+ resspec['inmodel'] = self._model
+ resspec['caldb'] = self._caldb
+ resspec['irf'] = self._irf
+ resspec['emin'] = 0.1
+ resspec['emax'] = 100.0
+ resspec['enumbins'] = 20
+ resspec['mask'] = False
+ resspec['algorithm'] = 'SUB'
+ resspec['logfile'] = 'csresspec_py1.log'
+ resspec['chatter'] = 2
+
+ # Run csresspec script
+ resspec.logFileOpen() # Make sure we get a log file
+ resspec.run()
+ resspec.save()
+
+ # Check result file
+ self._check_result_file('csresspec_py1.fits', 20, 5)
+
+ # Set-up csresspec for event list
+ # with computation of individual model components
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._events
+ resspec['outfile'] = 'csresspec_py2.fits'
+ resspec['inmodel'] = self._model
+ resspec['caldb'] = self._caldb
+ resspec['irf'] = self._irf
+ resspec['emin'] = 0.1
+ resspec['emax'] = 100.0
+ resspec['enumbins'] = 20
+ resspec['mask'] = False
+ resspec['algorithm'] = 'SUBDIV'
+ resspec['components'] = True
+ resspec['logfile'] = 'csresspec_py2.log'
+ resspec['chatter'] = 2
+
+ # Run csresspec script
+ resspec.execute()
+
+ # Check result file
+ self._check_result_file('csresspec_py2.fits', 20, 7)
+
+ # Set-up csresspec for counts cube
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._cntcube
+ resspec['modcube'] = 'NONE'
+ resspec['expcube'] = 'NONE'
+ resspec['psfcube'] = 'NONE'
+ resspec['edispcube'] = 'NONE'
+ resspec['bkgcube'] = 'NONE'
+ resspec['caldb'] = self._caldb
+ resspec['irf'] = self._irf
+ resspec['inmodel'] = self._model
+ resspec['outfile'] = 'csresspec_py3.fits'
+ resspec['mask'] = False
+ resspec['algorithm'] = 'SUBDIVSQRT'
+ resspec['logfile'] = 'csresspec_py3.log'
+ resspec['chatter'] = 3
+
+ # Run csresspec script
+ resspec.execute()
+
+ # Check result file
+ self._check_result_file('csresspec_py3.fits', 20, 5)
+
+ # Set-up csresspec for counts cube
+ # with source mask
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._cntcube
+ resspec['modcube'] = 'NONE'
+ resspec['expcube'] = 'NONE'
+ resspec['psfcube'] = 'NONE'
+ resspec['edispcube'] = 'NONE'
+ resspec['bkgcube'] = 'NONE'
+ resspec['caldb'] = self._caldb
+ resspec['irf'] = self._irf
+ resspec['inmodel'] = self._model
+ resspec['outfile'] = 'csresspec_py4.fits'
+ resspec['mask'] = True
+ resspec['ra'] = 83.63
+ resspec['dec'] = 22.01
+ resspec['rad'] = 0.2
+ resspec['regfile'] = 'NONE'
+ resspec['algorithm'] = 'SUBDIVSQRT'
+ resspec['logfile'] = 'csresspec_py4.log'
+ resspec['chatter'] = 3
+
+ # Run csresspec script
+ resspec.execute()
+
+ # Check result file
+ self._check_result_file('csresspec_py4.fits', 20, 5)
+
+ # Set-up csresspec for On/Off observations
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._inonoff
+ resspec['inmodel'] = self._onoff_model
+ resspec['outfile'] = 'csresspec_py5.fits'
+ resspec['algorithm'] = 'SIGNIFICANCE'
+ resspec['stack'] = False
+ resspec['logfile'] = 'csresspec_py5.log'
+ resspec['chatter'] = 3
+
+ # Run csresspec script
+ resspec.execute()
+
+ # Check result file
+ self._check_result_file('csresspec_py5.fits', 20, 8, nobs=2)
+
+ # Set-up csresspec for On/Off observations
+ # in stacked mode and with WSTAT background
+ resspec = cscripts.csresspec()
+ resspec['inobs'] = self._inonoff
+ resspec['inmodel'] = self._onoff_model
+ resspec['outfile'] = 'csresspec_py6.fits'
+ resspec['algorithm'] = 'SIGNIFICANCE'
+ resspec['statistic'] = 'WSTAT'
+ resspec['stack'] = True
+ resspec['logfile'] = 'csresspec_py6.log'
+ resspec['chatter'] = 3
+
+ # Run csresspec script
+ resspec.execute()
+
+ # Check result file
+ self._check_result_file('csresspec_py6.fits', 20, 8)
+
+ # Return
+ return
+
+ # Check result file
+ def _check_result_file(self, filename, enumbins, ncols, nobs=1):
+ """
+ Check result file
+ """
+ # Load residual file
+ fits = gammalib.GFits(filename)
+
+ # Check number of HDUs = observations (+ first empty image)
+ self.test_value(fits.size(), 1 + nobs, 'Number of observations')
+
+ # Check first hdu with residuals
+ table = fits[1]
+
+ # Check number of rows = number of energy bins
+ self.test_value(table.nrows(), enumbins,
+ 'Check for %d rows in residuals' % enumbins)
+
+ # Check number of columns, 2 for energy bounds + 3 for 3D or 6 for
+ # On/Off + 1 for each model component if components = yes
+ self.test_value(table.ncols(), ncols,
+ 'Check for %d columns in spectrum' % ncols)
+
+ # Return
+ return
diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py
index 756deba..3840030 100755
--- a/test/test_python_cscripts.py
+++ b/test/test_python_cscripts.py
@@ -41,6 +41,7 @@ import test_csphagen
import test_csphasecrv
import test_cspull
import test_csresmap
+import test_csresspec
import test_csroot2caldb
import test_cssens
import test_csspec
@@ -133,6 +134,7 @@ def test(installed=False, debug=False):
test_csphasecrv.Test(),
test_cspull.Test(),
test_csresmap.Test(),
+ test_csresspec.Test(),
test_csroot2caldb.Test(),
test_cssens.Test(),
test_csspec.Test(),