Simon Bonnefoy / ctools

Commit 47cd6656e630a3e062e68b2bd1518587734f00d6

Authored by Jürgen Knödlseder
2 parents edd8ccb7 c6f557e5

Merge branch '2289-add-csresspec' into integration 

# Conflicts:
#	doc/source/users/tutorials/beyond/index.rst
Inline Side-by-side
Showing 16 changed files with 1775 additions and 63 deletions
cscripts/Makefile.am
  View file @47cd665
... ... @@ -37,6 +37,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
37 37 $(srcdir)/csphasecrv.py \
38 38 $(srcdir)/cspull.py \
39 39 $(srcdir)/csresmap.py \
  40 + $(srcdir)/csresspec.py \
40 41 $(srcdir)/csroot2caldb.py \
41 42 $(srcdir)/cssens.py \
42 43 $(srcdir)/csspec.py \
... ... @@ -71,6 +72,7 @@ parfiles = $(srcdir)/cscaldb.par \
71 72 $(srcdir)/csphasecrv.par \
72 73 $(srcdir)/cspull.par \
73 74 $(srcdir)/csresmap.par \
  75 + $(srcdir)/csresspec.par \
74 76 $(srcdir)/csroot2caldb.par \
75 77 $(srcdir)/cssens.par \
76 78 $(srcdir)/csspec.par \
... ... @@ -102,6 +104,7 @@ linkscripts = cscaldb \
102 104 csphasecrv \
103 105 cspull \
104 106 csresmap \
  107 + csresspec \
105 108 csroot2caldb \
106 109 cssens \
107 110 csspec \
... ...
cscripts/__init__.py.in
  View file @47cd665
1 1 # ==========================================================================
2 2 # cscripts Python module
3 3 #
4   -# Copyright (C) 2015-2017 Juergen Knoedlseder
  4 +# Copyright (C) 2015-2018 Juergen Knoedlseder
5 5 #
6 6 # This program is free software: you can redistribute it and/or modify
7 7 # it under the terms of the GNU General Public License as published by
... ... @@ -40,6 +40,7 @@ __all__ = [
40 40 "csphasecrv",
41 41 "cspull",
42 42 "csresmap",
  43 + "csresspec",
43 44 "csroot2caldb",
44 45 "cssens",
45 46 "csspec",
... ... @@ -77,6 +78,7 @@ from cscripts.csphagen import csphagen
77 78 from cscripts.csphasecrv import csphasecrv
78 79 from cscripts.cspull import cspull
79 80 from cscripts.csresmap import csresmap
  81 +from cscripts.csresspec import csresspec
80 82 from cscripts.csroot2caldb import csroot2caldb
81 83 from cscripts.cssens import cssens
82 84 from cscripts.csspec import csspec
... ...
cscripts/csresmap.py
  View file @47cd665
... ... @@ -19,9 +19,9 @@
19 19 #
20 20 # ==========================================================================
21 21 import sys
22   -import math
23 22 import gammalib
24 23 import ctools
  24 +from cscripts import obsutils
25 25  
26 26  
27 27 # ============== #
... ... @@ -209,67 +209,14 @@ class csresmap(ctools.csobservation):
209 209 # Get model map into GSkyMap object
210 210 modelmap = model.cube().counts().copy()
211 211  
212   - # Store counts map as residual map. Note that we need a special
  212 + # Calculate residual maps
  213 + # Note that we need a special
213 214 # construct here to avoid memory leaks. This seems to be a SWIG feature
214 215 # as SWIG creates a new object when calling binning.cube()
215   - self._resmap = countmap.copy()
216   - self._resmap.stack_maps()
  216 + countmap1 = countmap.copy()
  217 + countmap1.stack_maps()
217 218 modelmap.stack_maps()
218   -
219   - # Get residual map algorithm type
220   - algorithm = self['algorithm'].string()
221   -
222   - # Subtract maps
223   - if algorithm == 'SUB':
224   - self._resmap -= modelmap
225   -
226   - # Subtract and divide by model map
227   - elif algorithm == 'SUBDIV':
228   - self._resmap -= modelmap
229   - self._resmap /= modelmap
230   -
231   - # Subtract and divide by sqrt of model map
232   - elif algorithm == 'SUBDIVSQRT':
233   - self._resmap -= modelmap
234   - self._resmap /= modelmap.sqrt()
235   -
236   - # Calculate significance from Li&Ma derivation
237   - elif algorithm == 'SIGNIFICANCE':
238   -
239   - # Compute sign map
240   - signmap = (self._resmap - modelmap).sign()
241   -
242   - # Loop over every bin in the map
243   - for i in range(self._resmap.npix()):
244   -
245   - # If the model value > 0.0 do the computation as normal ...
246   - model_val = modelmap[i]
247   - if model_val > 0.0:
248   -
249   - # If the data value is also > 0 then compute the
250   - # significance^2 and save it to the map ...
251   - data_val = self._resmap[i]
252   - if data_val > 0.0:
253   - log_val = math.log(data_val/model_val)
254   - self._resmap[i] = (data_val*log_val) + model_val - data_val
255   -
256   - # ... otherwise compute the reduced value of the above
257   - # expression. This is necessary to avoid computing log(0).
258   - else:
259   - self._resmap[i] = model_val
260   -
261   - # ... otherwise hard-code the significance to 0
262   - else:
263   - self._resmap[i] = 0.0
264   -
265   - # Compute significance map
266   - self._resmap *= 2.0
267   - self._resmap = self._resmap.sqrt()
268   - self._resmap = self._resmap * signmap
269   -
270   - # Raise exception if algorithm is unknown
271   - else:
272   - raise TypeError('Algorithm "'+algorithm+'" not known')
  219 + self._resmap = obsutils.residuals(self,countmap1,modelmap)
273 220  
274 221 # Optionally publish map
275 222 if self['publish'].boolean():
... ...
cscripts/csresspec.par 0 → 100644
  View file @47cd665
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: csresspec.par #
  8 +# #
  9 +# Author: Luigi Tibaldo #
  10 +# luigi.tibaldo@irap.omp.eu #
  11 +# IRAP #
  12 +# #
  13 +# Purpose: Parameter file for the csresspec script #
  14 +# #
  15 +#############################################################################
  16 +#
  17 +# Input and output parameters
  18 +#============================
  19 +inobs, f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
  20 +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
  21 +modcube, f, a, NONE,,, "Input model cube file (generated with ctmodel)"
  22 +expcube, f, a, NONE,,, "Input exposure cube file"
  23 +psfcube, f, a, NONE,,, "Input PSF cube file"
  24 +edispcube, f, a, NONE,,, "Input energy dispersion cube file"
  25 +bkgcube, f, a, NONE,,, "Input background cube file"
  26 +caldb, s, a, prod2,,, "Calibration database"
  27 +irf, s, a, South_0.5h,,, "Instrument response function"
  28 +edisp, b, h, no,,, "Apply energy dispersion?"
  29 +outfile, f, a, resspec.fits,,, "Output residual spectrum file"
  30 +
  31 +#
  32 +# Script parameters
  33 +#==================
  34 +statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
  35 +components,b, h, no,,, "Calculate model for individual components"
  36 +ebinalg, s, h, LOG,LIN|LOG|FILE,, "Binning algorithm"
  37 +emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  38 +emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  39 +enumbins, i, a, 20,,, "Number of model cube energy bins"
  40 +ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
  41 +stack, b, a, no,,, "Stack observations"
  42 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System"
  43 +proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
  44 +xref, r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
  45 +yref, r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
  46 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  47 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  48 +binsz, r, a, 0.02,,, "Pixel size (deg/pixel)"
  49 +mask, b, a, no,,, "Mask data to calculate residuals in ROI"
  50 +ra, r, a, 83.63,0,360, "RA for ROI centre (degrees)"
  51 +dec, r, a, 22.01,-90,90, "Dec for ROI centre (degrees)"
  52 +rad, r, a, 3.0,0,180, "Radius of ROI (degrees)"
  53 +regfile, f, a, "NONE",,, "Input exclusion region file in ds9 format"
  54 +algorithm, s, a, SIGNIFICANCE,SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE,, "Residuals computation algorithm"
  55 +
  56 +#
  57 +# Standard parameters
  58 +#====================
  59 +publish, b, h, no,,, "Publish residual map on VO Hub?"
  60 +chatter, i, h, 2,0,4, "Chattiness of output"
  61 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  62 +debug, b, h, no,,, "Debugging mode activated"
  63 +mode, s, h, ql,,, "Mode of automatic parameters"
  64 +logfile, f, h, csresspec.log,,, "Log filename"
... ...
cscripts/csresspec.py 0 → 100755
  View file @47cd665
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Generates a residual spectrum.
  4 +#
  5 +# Copyright (C) 2017-2018 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 ctools
  23 +from cscripts import obsutils
  24 +import sys
  25 +
  26 +
  27 +# =============== #
  28 +# csresspec class #
  29 +# =============== #
  30 +class csresspec(ctools.csobservation):
  31 + """
  32 + Generates a residual spectrum
  33 + """
  34 +
  35 + # Constructor
  36 + def __init__(self, *argv):
  37 + """
  38 + Constructor
  39 + """
  40 + # Initialise application by calling the appropriate class constructor
  41 + self._init_csobservation(self.__class__.__name__, ctools.__version__,
  42 + argv)
  43 +
  44 + # Initialise class members
  45 + self._use_maps = False
  46 + self._stack = False
  47 + self._mask = False
  48 + self._fits = gammalib.GFits()
  49 +
  50 + # Return
  51 + return
  52 +
  53 + # Private methods
  54 + def _get_parameters(self):
  55 + """
  56 + Get parameters from parfile and setup the observation
  57 + """
  58 + # Setup observations (require response and allow event list as well as
  59 + # counts cube)
  60 + self._setup_observations(self.obs(), True, True, True)
  61 +
  62 + # Set observation statistic
  63 + self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
  64 +
  65 + # Collect number of unbinned, binned and On/Off observations in
  66 + # observation container
  67 + n_unbinned = 0
  68 + n_binned = 0
  69 + n_onoff = 0
  70 + for obs in self.obs():
  71 + if obs.classname() == 'GCTAObservation':
  72 + if obs.eventtype() == 'CountsCube':
  73 + n_binned += 1
  74 + else:
  75 + n_unbinned += 1
  76 + elif obs.classname() == 'GCTAOnOffObservation':
  77 + n_onoff += 1
  78 + n_cta = n_unbinned + n_binned + n_onoff
  79 + n_other = self.obs().size() - n_cta
  80 +
  81 + # Query whether to compute model for individual components
  82 + components = self['components'].boolean()
  83 +
  84 + # If there is only one binned observation and no model for individual
  85 + # components is required, query for precomputed model file and set
  86 + # use_maps to True
  87 + if self.obs().size() == 1 and n_binned == 1 and not components:
  88 + modcube = self['modcube'].filename()
  89 + if modcube != 'NONE':
  90 + self._use_maps = True
  91 +
  92 + # If there are unbinned observations query the energy binning parameters
  93 + if n_unbinned != 0:
  94 + self['ebinalg'].string()
  95 + if self['ebinalg'].string() == 'FILE':
  96 + self['ebinfile'].filename()
  97 + else:
  98 + self['emin'].real()
  99 + self['emax'].real()
  100 + self['enumbins'].integer()
  101 + if n_cta > n_unbinned:
  102 + n_notunbin = n_cta - n_unbinned
  103 +
  104 + # If there is more than one observation, and observations are all
  105 + # unbinned or all onoff query user to know if they wish stacked results
  106 + if self.obs().size() > 1 and \
  107 + (n_unbinned == self.obs().size() or n_onoff == self.obs().size()):
  108 + self._stack = self['stack'].boolean()
  109 + # If we are to stack event lists query parameters for cube creation
  110 + if self._stack and n_unbinned == self.obs().size():
  111 + self['coordsys'].string()
  112 + self['proj'].string()
  113 + self['xref'].real()
  114 + self['yref'].real()
  115 + self['nxpix'].integer()
  116 + self['nypix'].integer()
  117 + self['binsz'].real()
  118 +
  119 + # If we are not using a precomputed model and no models are available
  120 + # in the observation container query input XML model file
  121 + if not self._use_maps and self.obs().models().size() == 0:
  122 + self.obs().models(self['inmodel'].filename())
  123 +
  124 + # Unless all observations are On/Off query for mask definition
  125 + if n_onoff == n_cta:
  126 + pass
  127 + else:
  128 + self._mask = self['mask'].boolean()
  129 + if self._mask:
  130 + self['ra'].real()
  131 + self['dec'].real()
  132 + self['rad'].real()
  133 + self['regfile'].filename()
  134 +
  135 + # Unless all observations are On/Off, or we are using precomputed model
  136 + # maps query whether to use energy dispersion
  137 + if n_onoff == n_cta or self._use_maps:
  138 + pass
  139 + else:
  140 + self['edisp'].boolean()
  141 +
  142 + # Query algorithm for residual computation
  143 + self['algorithm'].string()
  144 +
  145 + # Read ahead output parameters
  146 + if self._read_ahead():
  147 + self['outfile'].filename()
  148 +
  149 + # Write input parameters into logger
  150 + self._log_parameters(gammalib.TERSE)
  151 +
  152 + # Write header for observation census
  153 + self._log_header1(gammalib.TERSE, 'Observation census')
  154 +
  155 + # Log census of input observations
  156 + self._log_value(gammalib.NORMAL, 'Unbinned CTA observations', n_unbinned)
  157 + self._log_value(gammalib.NORMAL, 'Binned CTA observations', n_binned)
  158 + self._log_value(gammalib.NORMAL, 'On/off CTA observations', n_onoff)
  159 + self._log_value(gammalib.NORMAL, 'Other observations', n_other)
  160 + if n_other > 0:
  161 + msg = 'WARNING: Only CTA observation can be handled, all non-CTA ' \
  162 + + 'observations will be ignored.'
  163 + self._log_string(gammalib.TERSE, msg)
  164 +
  165 + # Log for unbinned observations
  166 + if n_unbinned != 0:
  167 + msg = ' User defined energy binning will be used for %d unbinned '\
  168 + 'observations.' % (n_unbinned)
  169 + self._log_string(gammalib.TERSE, msg)
  170 + if n_cta > n_unbinned:
  171 + msg = ' The intrinsic binning will be used for the remaining '\
  172 + '%d CTA observations.' % (n_notunbin)
  173 + self._log_string(gammalib.TERSE, msg)
  174 +
  175 + # Signal how energy dispersion is applied
  176 + if n_onoff == n_cta or self._use_maps:
  177 + msg = ' Energy dispersion is applied based on the input data/model ' \
  178 + + 'and not according to the edisp parameter'
  179 + self._log_string(gammalib.TERSE, msg)
  180 +
  181 + # Return
  182 + return
  183 +
  184 + def _stack_observations(self):
  185 + """
  186 + Stack multiple observations and replace observation container
  187 + with single stacked observation
  188 +
  189 + Returns
  190 + -------
  191 + obs : `~gammalib.GObservations'
  192 + Observation container with a single stacked observation
  193 + """
  194 + # Use first observation to determine the type and apply the
  195 + # appropriate stacking method
  196 +
  197 + # If On/Off use the GCTAOnOffObservation constructor
  198 + if self.obs()[0].classname() == 'GCTAOnOffObservation':
  199 + msg = 'Stacking %d On/Off observations.' % (self.obs().size())
  200 + self._log_string(gammalib.NORMAL, msg)
  201 + stacked = gammalib.GCTAOnOffObservation(self.obs())
  202 + stacked_obs = gammalib.GObservations()
  203 + stacked_obs.append(stacked)
  204 + stacked_obs.models(self.obs().models())
  205 +
  206 + # If event list then bin observations
  207 + elif self.obs()[0].classname() == 'GCTAObservation':
  208 + msg = 'Stacking %d event lists.' % (self.obs().size())
  209 + self._log_string(gammalib.NORMAL, msg)
  210 + stacked_obs = obsutils.get_stacked_obs(self, self.obs())
  211 +
  212 + # Return
  213 + return stacked_obs
  214 +
  215 + def _bin_evlist(self, obs):
  216 + """
  217 + Turn single event list into counts cube
  218 +
  219 + Parameters
  220 + ----------
  221 + obs : `~gammalib.GObservations()`
  222 + Observation container with single event list
  223 +
  224 + Returns
  225 + -------
  226 + obs : `~gammalib.GObservations()`
  227 + Binned observation container dict with event list ROI and energy
  228 + range information
  229 + """
  230 + # Retrieve information about ROI in event list
  231 + roi = obs[0].roi()
  232 + ra = roi.centre().dir().ra_deg()
  233 + dec = roi.centre().dir().dec_deg()
  234 + rad = roi.radius()
  235 +
  236 + # We will cover the whole ROI with 0.02 deg binning
  237 + npix = int(2.0 * rad / 0.02) + 1
  238 +
  239 + # Log binning of events
  240 + self._log_string(gammalib.EXPLICIT, 'Binning events')
  241 +
  242 + # Bin events
  243 + cntcube = ctools.ctbin(obs)
  244 + cntcube['xref'] = ra
  245 + cntcube['yref'] = dec
  246 + cntcube['binsz'] = 0.02
  247 + cntcube['nxpix'] = npix
  248 + cntcube['nypix'] = npix
  249 + cntcube['proj'] = 'TAN'
  250 + cntcube['coordsys'] = 'CEL'
  251 + cntcube['ebinalg'] = self['ebinalg'].string()
  252 + if self['ebinalg'].string() == 'FILE':
  253 + cntcube['ebinfile'] = self['ebinfile'].filename()
  254 + else:
  255 + cntcube['enumbins'] = self['enumbins'].integer()
  256 + cntcube['emin'] = self['emin'].real()
  257 + cntcube['emax'] = self['emax'].real()
  258 + cntcube.run()
  259 +
  260 + # Retrieve the binned observation container
  261 + binned_obs = cntcube.obs().copy()
  262 +
  263 + # Log response computation
  264 + self._log_string(gammalib.EXPLICIT, 'Computing binned response')
  265 +
  266 + # Compute binned response
  267 + response = obsutils.get_stacked_response(obs, ra, dec,
  268 + binsz=0.02,
  269 + nxpix=npix,
  270 + nypix=npix,
  271 + coordsys='CEL',
  272 + proj='TAN',
  273 + emin=self['emin'].real(),
  274 + emax=self['emax'].real(),
  275 + edisp=self['edisp'].boolean())
  276 +
  277 + if self['edisp'].boolean():
  278 + binned_obs[0].response(response['expcube'], response['psfcube'],
  279 + response['edispcube'], response['bkgcube'])
  280 + else:
  281 + binned_obs[0].response(response['expcube'], response['psfcube'],
  282 + response['bkgcube'])
  283 +
  284 + # Get new models
  285 + models = response['models']
  286 +
  287 + # Set models for observation container
  288 + binned_obs.models(models)
  289 +
  290 + # Check if energy boundaries provided by user extend beyond
  291 + # the content of the event list
  292 + if self['emin'].real() > obs[0].events().emin().TeV():
  293 + emin = 'INDEF'
  294 + else:
  295 + emin = obs[0].events().emin().TeV()
  296 + if self['emax'].real() < obs[0].events().emax().TeV():
  297 + emax = 'INDEF'
  298 + else:
  299 + emax = obs[0].events().emax().TeV()
  300 +
  301 + # Put ROI and E bound info in dictionary
  302 + info = {'was_list': True, 'roi_ra': ra, 'roi_dec': dec, 'roi_rad': rad,
  303 + 'emin': emin, 'emax': emax}
  304 +
  305 + # Return new oberservation container
  306 + return binned_obs, info
  307 +
  308 + def _masked_cube(self, cube, ra, dec, rad, emin='INDEF', emax='INDEF',
  309 + regfile='NONE'):
  310 + """
  311 + Mask an event cube and returns the masked cube
  312 +
  313 + Parameters
  314 + ----------
  315 + cube : `~gammalib.GCTAEventCube()`
  316 + Event cube
  317 + ra : float (str 'INDEF' for no selection on direction)
  318 + Right Ascension (deg)
  319 + dec : float (str 'INDEF' for no selection on direction)
  320 + Declination (deg)
  321 + rad : float (str 'INDEF' for no selection on direction)
  322 + Radius (deg)
  323 + emin : float (str 'INDEF' for no selection on energy)
  324 + Minimum energy (TeV)
  325 + emax : float (str 'INDEF' for no selection on energy)
  326 + Maximum energy (TeV)
  327 +
  328 + Returns
  329 + -------
  330 + cube : `~gammalib.GCTAEventCube()`
  331 + Event cube
  332 + """
  333 + # Turn cube into observation container to feed to ctcubemask
  334 + obs = gammalib.GCTAObservation()
  335 + obs.events(cube)
  336 + obs_cont = gammalib.GObservations()
  337 + obs_cont.append(obs)
  338 +
  339 + # Use ctcubemask to mask event cube pixels
  340 + cubemask = ctools.ctcubemask(obs_cont)
  341 + cubemask['ra'] = ra
  342 + cubemask['dec'] = dec
  343 + cubemask['rad'] = rad
  344 + cubemask['emin'] = emin
  345 + cubemask['emax'] = emax
  346 + cubemask['regfile'] = regfile
  347 + cubemask.run()
  348 +
  349 + # Extract copy of cube from observation container (copy is needed to
  350 + # avoid memory leaks in SWIG)
  351 + cube = cubemask.obs()[0].events().copy()
  352 +
  353 + # Return cube
  354 + return cube
  355 +
  356 + def _cube_to_spectrum(self, cube, evlist_info):
  357 + """
  358 + Derive from event cube a count spectrum
  359 +
  360 + If data come from event list use only the ROI and energy range of
  361 + the original data. Apply user defined mask if requested.
  362 +
  363 + Parameters
  364 + ----------
  365 + cube : `~gammalib.GCTAEventCube()`
  366 + Event cube
  367 + evlist_info : dict
  368 + Dictionary with information on original event list
  369 +
  370 + Returns
  371 + -------
  372 + array : `~gammalib.GNdarray()'
  373 + Counts spectrum
  374 + """
  375 + # If we started from event list mask the ROI only
  376 + # for residual computation
  377 + if evlist_info['was_list']:
  378 + msg = 'Masking ROI from original event list'
  379 + self._log_string(gammalib.EXPLICIT, msg)
  380 + cube = self._masked_cube(cube, evlist_info['roi_ra'],
  381 + evlist_info['roi_dec'],
  382 + evlist_info['roi_rad'],
  383 + emin=evlist_info['emin'],
  384 + emax=evlist_info['emax'])
  385 +
  386 + # Apply user mask
  387 + if self._mask:
  388 + msg = 'Masking ROI requested by user'
  389 + self._log_string(gammalib.EXPLICIT, msg)
  390 + cube = self._masked_cube(cube, self['ra'], self['dec'],
  391 + self['rad'],
  392 + regfile=self['regfile'])
  393 +
  394 + # Extract skymap and clip at 0 to null masked areas
  395 + counts = cube.counts().copy()
  396 + counts = counts.clip(0.)
  397 +
  398 + # Convert skymap into GNdarray count spectrum
  399 + counts = counts.counts()
  400 +
  401 + # Return
  402 + return counts
  403 +
  404 + def _residuals_table(self, obs_id, ebounds, counts, model, residuals):
  405 + """
  406 + Create a Fits Table and store counts, model, and residuals
  407 +
  408 + Parameters
  409 + ----------
  410 + obs_id : str
  411 + Observation id
  412 + ebounds : `~gammalib.GEbounds()'
  413 + Energy boundaries
  414 + counts : `~gammalib.GNdarray()'
  415 + Counts spectrum
  416 + model : `~gammalib.GNdarray()'
  417 + Model spectrum
  418 + residuals : `~gammalib.GNdarray()'
  419 + Residual spectrum
  420 +
  421 + Returns
  422 + -------
  423 + table : `~gammalib.GFitsBinTable()'
  424 + Residual spectrum as FITS binary table
  425 + """
  426 + # Create FITS table columns
  427 + nrows = ebounds.size()
  428 + energy_low = gammalib.GFitsTableDoubleCol('Emin', nrows)
  429 + energy_high = gammalib.GFitsTableDoubleCol('Emax', nrows)
  430 + counts_col = gammalib.GFitsTableDoubleCol('Counts', nrows)
  431 + model_col = gammalib.GFitsTableDoubleCol('Model', nrows)
  432 + resid_col = gammalib.GFitsTableDoubleCol('Residuals', nrows)
  433 + energy_low.unit('TeV')
  434 + energy_high.unit('TeV')
  435 +
  436 + # Fill FITS table columns
  437 + for i in range(nrows):
  438 + energy_low[i] = ebounds.emin(i).TeV()
  439 + energy_high[i] = ebounds.emax(i).TeV()
  440 + counts_col[i] = counts[i]
  441 + model_col[i] = model[i]
  442 + resid_col[i] = residuals[i]
  443 +
  444 + # Initialise FITS Table with extension set to obs id
  445 + table = gammalib.GFitsBinTable(nrows)
  446 + # If name not empty add leading blank
  447 + if obs_id != '':
  448 + obs_id = ' ' + obs_id
  449 + table.extname('RESIDUALS' + obs_id)
  450 +
  451 + # Add Header card to specify algorithm used for residual computation
  452 + table.card('ALGORITHM', self['algorithm'].string(),
  453 + 'Algorithm used for computation of residuals')
  454 +
  455 + # Append filled columns to fits table
  456 + table.append(energy_low)
  457 + table.append(energy_high)
  458 + table.append(counts_col)
  459 + table.append(model_col)
  460 + table.append(resid_col)
  461 +
  462 + # Return binary table
  463 + return table
  464 +
  465 + def _append_column(self, table, name, data):
  466 + """
  467 + Append optional column to residual table
  468 +
  469 + Parameters
  470 + ----------
  471 + table : `~gammalib.GFitsBinTable()'
  472 + FITS binary table
  473 + name : str
  474 + Column name
  475 + data : `~gammalib.GEbounds()'
  476 + Data to be filled into new column
  477 +
  478 + Returns
  479 + -------
  480 + table : `~gammalib.GFitsBinTable()'
  481 + FITS binary table
  482 + """
  483 + # Check size compatibility
  484 + if table.nrows() == data.size():
  485 + pass
  486 + # Otherwise throw error
  487 + else:
  488 + msg = 'csresspec._append_column: FITS table and data have ' \
  489 + 'incompatible size.'
  490 + raise RuntimeError(msg)
  491 +
  492 + # Create column
  493 + column = gammalib.GFitsTableDoubleCol(name, table.nrows())
  494 +
  495 + # Fill data
  496 + for i, value in enumerate(data):
  497 + column[i] = value
  498 +
  499 + # Append new column to table
  500 + table.append(column)
  501 +
  502 + # Return modified table
  503 + return table
  504 +
  505 +
  506 + # Public methods
  507 + def run(self):
  508 + """
  509 + Run the script
  510 + """
  511 + # Switch screen logging on in debug mode
  512 + if self._logDebug():
  513 + self._log.cout(True)
  514 +
  515 + # Get parameters
  516 + self._get_parameters()
  517 +
  518 + # Write observation into logger
  519 + self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
  520 +
  521 + # Stack observations if requested
  522 + if self._stack:
  523 + self._log_header1(gammalib.TERSE, 'Stacking Observations')
  524 + self.obs(self._stack_observations())
  525 +
  526 + # Log processing header
  527 + self._log_header1(gammalib.TERSE, 'Processing Observations')
  528 +
  529 + # Loop over observations and calculate residuals
  530 + for s, observation in enumerate(self.obs()):
  531 +
  532 + # Retrieve and store obs id
  533 + obs_id = observation.id()
  534 +
  535 + # If observation id is empty and there is more than one observation
  536 + # replace with incremental number
  537 + if obs_id == '' and self.obs().size() > 1:
  538 + obs_id = str(s)
  539 +
  540 + # Log processing of observation
  541 + if self.obs().size() > 1:
  542 + self._log_header2(gammalib.NORMAL,
  543 + 'Processing observation %s' % obs_id)
  544 +
  545 + # Turn into observation container and assign models
  546 + obs = gammalib.GObservations()
  547 + obs.append(observation)
  548 + obs.models(self.obs().models())
  549 +
  550 + # if 3D observation
  551 + if obs[0].classname() == 'GCTAObservation':
  552 +
  553 + # Prepare Observations
  554 +
  555 + # If already binned set the evlist_info dictionary to have
  556 + # attribute was_list False
  557 + if obs[0].eventtype() == 'CountsCube':
  558 + evlist_info = {'was_list': False}
  559 + # Otherwise bin now
  560 + else:
  561 + # we remember if we binned an event list
  562 + # so that we can mask only the ROI for residual calculation
  563 + msg = 'Setting up binned observation'
  564 + self._log_string(gammalib.NORMAL, msg)
  565 + obs, evlist_info = self._bin_evlist(obs)
  566 +
  567 + # Calculate Model and residuals
  568 + # If model cube is provided load it
  569 + if self._use_maps:
  570 + modcube = gammalib.GCTAEventCube(self['inmodel'].filename())
  571 + # Otherwise calculate it now
  572 + else:
  573 + msg = 'Computing model cube'
  574 + self._log_string(gammalib.NORMAL, msg)
  575 + modelcube = ctools.ctmodel(obs)
  576 + modelcube.run()
  577 + modcube = modelcube.cube().copy()
  578 +
  579 + # Extract cntcube for residual computation
  580 + cntcube = obs[0].events().copy()
  581 +
  582 + # Derive count spectra from cubes
  583 + msg = 'Computing counts, model, and residual spectra'
  584 + self._log_string(gammalib.NORMAL, msg)
  585 + counts = self._cube_to_spectrum(cntcube, evlist_info)
  586 + model = self._cube_to_spectrum(modcube, evlist_info)
  587 +
  588 + # Calculate residuals
  589 + residuals = obsutils.residuals(self, counts, model)
  590 +
  591 + # Extract energy bounds
  592 + ebounds = cntcube.ebounds()
  593 +
  594 + # Fill results table
  595 + msg = 'Filling residual table'
  596 + self._log_string(gammalib.NORMAL, msg)
  597 + table = self._residuals_table(obs_id, ebounds, counts, model,
  598 + residuals)
  599 +
  600 + # Calculate models of individual components if requested
  601 + if self['components'].boolean():
  602 + for component in obs.models():
  603 + self._log_value(gammalib.NORMAL,
  604 + 'Computing model for component',
  605 + component.name())
  606 + # Set model cube models to
  607 + # individual component
  608 + model_cont = gammalib.GModels()
  609 + model_cont.append(component)
  610 + modelcube.obs().models(model_cont)
  611 +
  612 + # Run model cube
  613 + modelcube.run()
  614 +
  615 + # Extract spectrum of individual component
  616 + modcube = modelcube.cube().copy()
  617 + model = self._cube_to_spectrum(modcube, evlist_info)
  618 +
  619 + # append component to table
  620 + table = self._append_column(table, component.name(),
  621 + model)
  622 +
  623 + # otherwise, if On/Off
  624 + elif obs[0].classname() == 'GCTAOnOffObservation':
  625 +
  626 + onoff = obs[0]
  627 +
  628 + # Calculate Counts, Model and residuals
  629 + msg = 'Computing counts, model, and residual spectra'
  630 + self._log_string(gammalib.NORMAL, msg)
  631 +
  632 + # On spectrum
  633 + counts = onoff.on_spec().counts_spectrum()
  634 +
  635 + # Model for On region
  636 + background = onoff.model_background(obs.models()).counts_spectrum()
  637 + alpha = onoff.on_spec().backscal_spectrum()
  638 + model = background.copy()
  639 + model *= alpha
  640 + model += onoff.model_gamma(obs.models()).counts_spectrum()
  641 +
  642 + # On Residuals
  643 + residuals = obsutils.residuals(self, counts, model)
  644 +
  645 + # Extract energy bounds
  646 + ebounds = onoff.on_spec().ebounds()
  647 +
  648 + # Fill results table
  649 + msg = 'Filling residual table'
  650 + self._log_string(gammalib.NORMAL, msg)
  651 + table = self._residuals_table(obs_id, ebounds, counts, model,
  652 + residuals)
  653 +
  654 + # Get Off spectrum and add to table
  655 + msg = 'Computing counts, model, and residual spectra for Off regions'
  656 + self._log_string(gammalib.NORMAL, msg)
  657 + counts_off = onoff.off_spec().counts_spectrum()
  658 + table = self._append_column(table, 'Counts_Off',
  659 + counts_off)
  660 +
  661 + # Add background/Off model to table
  662 + table = self._append_column(table, 'Model_Off',
  663 + background)
  664 +
  665 + # Calculate Off residuals and add to table
  666 + residuals_off = obsutils.residuals(self, counts_off, background)
  667 + table = self._append_column(table, 'Residuals_Off',
  668 + residuals_off)
  669 +
  670 + # Calculate models of individual components if requested
  671 + if self['components'].boolean():
  672 + for component in obs.models():
  673 + # If background pass
  674 + # We always add the background at the end so that
  675 + # we accommodate WSTAT for which the background is not
  676 + # mandatory in the model
  677 + if component.classname() == 'GCTAModelIrfBackground':
  678 + pass
  679 + # Otherwise calculate gamma component and append to Table
  680 + else:
  681 + self._log_value(gammalib.NORMAL,
  682 + 'Computing model for component',
  683 + component.name())
  684 + # Create observation container for individual components
  685 + model_cont = gammalib.GModels()
  686 + model_cont.append(component)
  687 + # Calculate gamma model
  688 + model = onoff.model_gamma(model_cont)
  689 + model = model.counts_spectrum()
  690 + # Append to table
  691 + table = self._append_column(table, component.name(),
  692 + model)
  693 + # Add background already calculated
  694 + self._log_value(gammalib.NORMAL,
  695 + 'Computing model for component',
  696 + 'Background')
  697 + background *= onoff.on_spec().backscal_spectrum()
  698 + table = self._append_column(table, 'Background',
  699 + background)
  700 +
  701 + # Append results table to output file
  702 + self._fits.append(table)
  703 +
  704 + # Return
  705 + return
  706 +
  707 + def save(self):
  708 + """
  709 + Save residuals
  710 + """
  711 + # Write header
  712 + self._log_header1(gammalib.TERSE, 'Save residuals')
  713 +
  714 + # Continue only if FITS file is valid
  715 + if self._fits != None:
  716 +
  717 + # Get outfile parameter
  718 + outfile = self['outfile'].filename()
  719 +
  720 + # Log file name
  721 + self._log_value(gammalib.NORMAL, 'Residuals file', outfile.url())
  722 +
  723 + # Save residuals
  724 + self._fits.saveto(outfile, self['clobber'].boolean())
  725 +
  726 + # Return
  727 + return
  728 +
  729 +
  730 +# ======================== #
  731 +# Main routine entry point #
  732 +# ======================== #
  733 +if __name__ == '__main__':
  734 +
  735 + # Create instance of application
  736 + app = csresspec(sys.argv)
  737 +
  738 + # Execute application
  739 + app.execute()
... ...
cscripts/obsutils.py
  View file @47cd665
1 1 # ==========================================================================
2 2 # Utility functions for observation handling
3 3 #
4   -# Copyright (C) 2011-2017 Juergen Knoedlseder
  4 +# Copyright (C) 2011-2018 Juergen Knoedlseder
5 5 #
6 6 # This program is free software: you can redistribute it and/or modify
7 7 # it under the terms of the GNU General Public License as published by
... ... @@ -829,3 +829,102 @@ def get_onoff_obs(cls, obs):
829 829  
830 830 # Return On/Off oberservation container
831 831 return onoff_obs
  832 +
  833 +
  834 +# ========================================= #
  835 +# Calculate residuals from counts and model #
  836 +# ========================================= #
  837 +def residuals(cls, counts, model):
  838 + """
  839 + Calculate residuals given counts and models, according to algorithm
  840 + specified by user.
  841 +
  842 + Can handle GSkyMap or GNdarray objects
  843 +
  844 +
  845 + Parameters
  846 + ----------
  847 + cls : `~ctools.cscript`
  848 + cscript class
  849 + counts : `~gammalib.GSkyMap/~gammalib.GNdarray'
  850 + Data counts
  851 + model : `~gammalib.GSkyMap/~gammalib.GNdarray'
  852 + Model counts
  853 +
  854 + Returns
  855 + -------
  856 + residuals : `~gammalib.GSkyMap/~gammalib.GNdarray'
  857 + Residuals
  858 + """
  859 + # Find type of objects we are manipulating and set size to iterate later
  860 +
  861 + # If GNdarray
  862 + if counts.classname() == 'GNdarray':
  863 + nelem = counts.size()
  864 + elif counts.classname() == 'GSkyMap':
  865 + nelem = counts.npix()
  866 + else:
  867 + msg = 'cscripts.obsutils.residuals only handles ' \
  868 + + 'gammalib.GNdarray or gammalib.GSkyMap objects.\n'
  869 + raise RuntimeError(msg)
  870 +
  871 + # Copy counts to initialise residuals
  872 + residuals = counts.copy()
  873 +
  874 + # Get residual map algorithm type
  875 + algorithm = cls['algorithm'].string()
  876 +
  877 + # Subtract
  878 + if algorithm == 'SUB':
  879 + residuals -= model
  880 +
  881 + # Subtract and divide by model
  882 + elif algorithm == 'SUBDIV':
  883 + residuals -= model
  884 + residuals /= model
  885 +
  886 + # Subtract and divide by sqrt of model
  887 + elif algorithm == 'SUBDIVSQRT':
  888 + residuals -= model
  889 + residuals /= model.sqrt()
  890 +
  891 + # Calculate significance from Li&Ma
  892 + elif algorithm == 'SIGNIFICANCE':
  893 +
  894 + # Compute sign
  895 + sign = (residuals - model).sign()
  896 +
  897 + # Loop over every bin
  898 + for i in range(nelem):
  899 +
  900 + # If the model value > 0.0 do the computation as normal ...
  901 + model_val = model[i]
  902 + if model_val > 0.0:
  903 +
  904 + # If the data value is also > 0 then compute the
  905 + # significance^2 and save it ...
  906 + data_val = residuals[i]
  907 + if data_val > 0.0:
  908 + log_val = math.log(data_val / model_val)
  909 + residuals[i] = (data_val * log_val) + model_val - data_val
  910 +
  911 + # ... otherwise compute the reduced value of the above
  912 + # expression. This is necessary to avoid computing log(0).
  913 + else:
  914 + residuals[i] = model_val
  915 +
  916 + # ... otherwise hard-code the significance to 0
  917 + else:
  918 + residuals[i] = 0.0
  919 +
  920 + # Compute significance map
  921 + residuals *= 2.0
  922 + residuals = residuals.sqrt()
  923 + residuals *= sign
  924 +
  925 + # Raise exception if algorithm is unknown
  926 + else:
  927 + raise TypeError('Algorithm "' + algorithm + '" not known')
  928 +
  929 + # Return
  930 + return residuals
... ...
doc/source/users/reference_manual/csresspec.rst 0 → 100644
  View file @47cd665
  1 +.. _csresspec:
  2 +
  3 +csresspec
  4 +=========
  5 +
  6 +Generates residual spectrum.
  7 +
  8 +
  9 +Synopsis
  10 +--------
  11 +
  12 +This scripts generates a residual spectrum for a given model. It works for
  13 +event lists, counts cubes, On/Off observations or observation definition files.
  14 +For event lists parameters that define the spectral binning need to
  15 +be provided so that the script can bin the data internally. The model
  16 +is then convolved with the instrumental response function for the chosen
  17 +binning (either specified by the user or intrinsic to the observations)
  18 +and used for residual computation. It is possible to compute the model
  19 +contributions from individual model components.
  20 +
  21 +If starting from event lists or count cubes,
  22 +before residual computation the 3D counts and model cubes are collapsed into
  23 +spectra by summing over all spatial direction. It is possible to specify a mask
  24 +to compute counts, model, and residuals only for a region of interest.
  25 +
  26 +For On/Off observations, counts, model, and residual spectra for the Off regions
  27 +are computed as well.
  28 +
  29 +For an observation definition file, when all observations are event lists or
  30 +On/Off observations it is possible to stack them before residual computation. To
  31 +do so spatial binning parameters for event lists are required. Otherwise the
  32 +residuals are calculated separately for each observation.
  33 +
  34 +Four options exist then for residual computation:
  35 +
  36 +* the subtraction of the model from the counts (SUB)
  37 +* the subtraction and division by the model (SUBDIV)
  38 +* the subtraction and division by the square root of the model (SUBDIVSQRT)
  39 +* the computation of the significance using Li & Ma (SIGNIFICANCE)
  40 +
  41 +The counts, model, and residual spectra are written into a FITS file.
  42 +
  43 +
  44 +General parameters
  45 +------------------
  46 +
  47 +``inobs [file]``
  48 + Input event list, counts cube or observation definition XML file.
  49 +
  50 +``modcube [file]``
  51 + Input model cube file (generated with ctmodel).
  52 +
  53 +``inmodel [file]``
  54 + Input model XML file.
  55 +
  56 +``expcube [file]``
  57 + Input exposure cube file.
  58 +
  59 +``psfcube [file]``
  60 + Input PSF cube file.
  61 +
  62 +``edispcube [file]``
  63 + Input energy dispersion cube file.
  64 +
  65 +``bkgcube [file]``
  66 + Input background cube file.
  67 +
  68 +``caldb [string]``
  69 + Calibration database.
  70 +
  71 +``irf [string]``
  72 + Instrument response function.
  73 +
  74 +``(edisp = no) [boolean]``
  75 + Apply energy dispersion?
  76 +
  77 +``outfile [file]``
  78 + Output file.
  79 +
  80 +``(statistic = DEFAULT) <DEFAULT|CSTAT|WSTAT> [string]``
  81 + Optimization statistic. For On/Off observations determines how
  82 + the background counts are computed. ``DEFAULT``, i.e., ``CSTAT``,
  83 + will use the input background model. ``WSTAT`` uses the background
  84 + counts from likelihood profiling as defined for the W statistic.
  85 +
  86 +``ebinalg <FILE|LIN|LOG> [string]``
  87 + Algorithm for defining energy bins.
  88 +
  89 +``emin [real]``
  90 + Lower energy value for first energy bin (in TeV).
  91 +
  92 +``emax [real]``
  93 + Upper energy value for last energy bin (in TeV).
  94 +
  95 +``enumbins [integer]``
  96 + Number of energy bins.
  97 +
  98 +``ebinfile [file]``
  99 + Name of the file containing the energy bin definition.
  100 +
  101 +``(stack = no) [boolean]``
  102 + Stack multiple observations for residual calculation.
  103 +
  104 +``coordsys <CEL|GAL> [string]``
  105 + Coordinate system (CEL - celestial, GAL - galactic) for stacking
  106 + multiple event lists.
  107 +
  108 +``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
  109 + Projection method for stacking multiple event lists.
  110 +
  111 +``xref [real]``
  112 + Right Ascension / Galactic longitude of cube centre (J2000, in degrees)
  113 + for stacking multiple event lists.
  114 +
  115 +``yref [real]``
  116 + Declination / Galactic latitude of cube centre (J2000, in degrees) for stacking
  117 + multiple event lists.
  118 +
  119 +``nxpix [integer]``
  120 + Number of cube bins in Right Ascension or Galactic longitude for stacking
  121 + multiple event lists.
  122 +
  123 +``nypix [integer]``
  124 + Number of cube bins in Declination or Galactic latitude for stacking
  125 + multiple event lists.
  126 +
  127 +``binsz [real]``
  128 + Cube bin size (in degrees/pixel) for stacking
  129 + multiple event lists.
  130 +
  131 +``ra [real]``
  132 + Right Ascension of circular selection region centre (J2000, in degrees).
  133 +
  134 +``dec [real]``
  135 + Declination of circular selection region centre (J2000, in degrees).
  136 +
  137 +``rad [real]``
  138 + Radius of circular selection region (in degrees).
  139 +
  140 +``regfile [file]``
  141 + Input exclusion region file in ds9 format.
  142 +
  143 +``algorithm <SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE> [string]``
  144 + Algorithm used to generate the residual map:
  145 +
  146 + ``SUB``: :math:`DATA - MODEL`
  147 +
  148 + ``SUBDIV``: :math:`(DATA - MODEL) / MODEL`
  149 +
  150 + ``SUBDIVSQRT``: :math:`(DATA - MODEL) / \sqrt{MODEL}`
  151 +
  152 + ``SIGNIFICANCE``: :math:`{\rm sign}(DATA-MODEL) \times \sqrt{ 2 \times ( DATA \times \ln \left(\frac{DATA}{MODEL} \right) + MODEL - DATA ) }`
  153 +
  154 +
  155 +Standard parameters
  156 +-------------------
  157 +
  158 +``(publish = no) [boolean]``
  159 + Specifies whether the residual map should be published on VO Hub.
  160 +
  161 +``(chatter = 2) [integer]``
  162 + Verbosity of the executable:
  163 + ``chatter = 0``: no information will be logged
  164 +
  165 + ``chatter = 1``: only errors will be logged
  166 +
  167 + ``chatter = 2``: errors and actions will be logged
  168 +
  169 + ``chatter = 3``: report about the task execution
  170 +
  171 + ``chatter = 4``: detailed report about the task execution
  172 +
  173 +``(clobber = yes) [boolean]``
  174 + Specifies whether an existing residual map file should be overwritten.
  175 +
  176 +``(debug = no) [boolean]``
  177 + Enables debug mode. In debug mode the executable will dump any log file output to the console.
  178 +
  179 +``(mode = ql) [string]``
  180 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
  181 +
  182 +``(logfile = csresmap.log) [string]``
  183 + Log filename.
  184 +
  185 +
  186 +Related tools or scripts
  187 +------------------------
  188 +
  189 +None
... ...
doc/source/users/reference_manual/index.rst
  View file @47cd665
... ... @@ -57,6 +57,7 @@ cscripts
57 57 csphasecrv --- Computes phase curve <csphasecrv>
58 58 cspull --- Generates pull distribution <cspull>
59 59 csresmap --- Generates residual map <csresmap>
  60 + csresspec --- Generates residual spectrum <csresspec>
60 61 csroot2caldb --- Creates a caldb entry from a ROOT file <csroot2caldb>
61 62 cssens --- Computes CTA sensitivity <cssens>
62 63 csspec --- Computes spectral points <csspec>
... ...
doc/source/users/tutorials/beyond/beyond_residual_spectrum.rst 0 → 100644
  View file @47cd665
  1 +.. _sec_residual_spectrum:
  2 +
  3 +Generating a residual spectrum
  4 +~~~~~~~~~~~~~~~~~~~~~~~~~
  5 +
  6 +In the previous tutorial you have used the :ref:`csresmap` script to
  7 +inspect the spatial residuals. You can also inspect the spectral
  8 +residuals using the :ref:`csresspec` script.
  9 +
  10 +The :ref:`csresspec` script can deal with any observation type (unbinned,
  11 +binned, stacked, On/Off). For the example case with unbinned
  12 +observations the script is then used as follows:
  13 +
  14 +.. code-block:: bash
  15 +
  16 + $csresspec components=yes
  17 + Input event list, counts cube, or observation definition XML file [obs.xml]
  18 + Calibration database [prod2]
  19 + Instrument response function [South_0.5h]
  20 + Lower energy limit (TeV) [0.1]
  21 + Upper energy limit (TeV) [100.0]
  22 + Number of model cube energy bins [30]
  23 + Stack observations [yes]
  24 + Coordinate System (CEL|GAL) [CEL]
  25 + Projection method (AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN) [CAR]
  26 + First coordinate of image/source region center in degrees (RA or galactic l) (0-360) [83.63]
  27 + Second coordinate of image/source region center in degrees (DEC or galactic b) (-90-90) [22.01]
  28 + Size of the X axis in pixels [200] 51
  29 + Size of the Y axis in pixels [200] 51
  30 + Pixel size (deg/pixel) [0.02]
  31 + Input model definition XML file [crab_results.xml]
  32 + Mask data to calculate residuals in ROI [yes] no
  33 + Residuals computation algorithm (SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE) [SIGNIFICANCE]
  34 + Output residual spectrum file [residual_spectrum.fits]
  35 +
  36 +In our example, the residuals will cover energies from 100 GeV to 100
  37 +TeV, in 30 energy bins. Since we have decided to stack multiple
  38 +unbinned observations we have also entered parameters to perform a
  39 +spatial binning. Note that we have defined the spatial grid such that
  40 +we will obtain the residuals in a 1 degree square region around the
  41 +Crab. Optionally, a mask can be applied to either unbinned or binned
  42 +3D observations. We have used the hidden ``components`` parameter to
  43 +include the individual model components in the results file.
  44 +
  45 +The :ref:`csresspec` script created the FITS file ``residual_spectrum.fits`` that contains
  46 +the results of the residual computation. A graphical display of the fit results is
  47 +shown below.
  48 +
  49 +.. figure:: residual_spectrum.png
  50 + :width: 600px
  51 + :align: center
  52 +
  53 + *Data, model, and residuals from the spectral fit of the Crab nebula*
  54 +
  55 +..
  56 +
  57 +
  58 +.. note::
  59 +
  60 + The figure was created by typing:
  61 +
  62 + .. code-block:: bash
  63 +
  64 + $ $CTOOLS/share/examples/python/show_spectrum.py residual_spectrum.fits
  65 +
  66 + Use the option ``-h`` to select the table to display if there are
  67 + multiple observations.
  68 +
  69 +
  70 +
  71 +
  72 +
... ...
doc/source/users/tutorials/beyond/index.rst
  View file @47cd665
... ... @@ -17,6 +17,7 @@ analyses.
17 17 stacked-analysis
18 18 beyond_tsmap
19 19 beyond_residual
  20 + beyond_residual_spectrum
20 21 advanced-python-fitting
21 22 beyond_irf
22 23 beyond_model
... ...
doc/source/users/tutorials/beyond/residual_spectrum.png 0 → 100644
View file @47cd665

70.6 KB
examples/Makefile.am
  View file @47cd665
1 1 #############################################################################
2 2 # Makefile for ctools examples #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2017 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2018 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -28,6 +28,7 @@ scripts = README.md \
28 28 show_pha.py \
29 29 show_pull_evolution.py \
30 30 show_pull_histogram.py \
  31 + show_residuals.py \
31 32 show_response.py \
32 33 show_sensitivity.py \
33 34 show_spectrum.py \
... ...
examples/show_residuals.py 0 → 100755
  View file @47cd665
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Display residuals generated by csresspec
  4 +#
  5 +# Copyright (C) 2017-2018 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 sys
  22 +import math
  23 +import gammalib
  24 +import cscripts
  25 +try:
  26 + import matplotlib.pyplot as plt
  27 + from matplotlib import gridspec
  28 + plt.figure()
  29 + plt.close()
  30 +except (ImportError, RuntimeError):
  31 + print('This script needs the "matplotlib" module')
  32 + sys.exit()
  33 +
  34 +
  35 +# ================================= #
  36 +# Fill counts/model/residual arrays #
  37 +# ================================= #
  38 +def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
  39 + c_resid, algorithm):
  40 + """
  41 + Helper function that fills in the counts, model, residuals lists and
  42 + calculate relative errors for given table row
  43 +
  44 + Parameters
  45 + ----------
  46 + row : int
  47 + Table row index
  48 + counts : list of float
  49 + Counts
  50 + model : list of float
  51 + Model
  52 + resid : list of float
  53 + Residual
  54 + e_counts : list of float
  55 + Counts error
  56 + e_resid : list of float
  57 + Residual error
  58 + c_counts : `~gammalib.GFitsTableCol()'
  59 + Counts FITS table column
  60 + c_model : `~gammalib.GFitsTableCol()'
  61 + Model FITS table column
  62 + c_resid : `~gammalib.GFitsTableCol()'
  63 + Residual FITS table column
  64 + algorithm : str
  65 + Algorithm
  66 +
  67 + Returns
  68 + -------
  69 + counts, model, resid, e_counts, e_resid : tuple of float
  70 + Values
  71 + """
  72 + # Extract values
  73 + counts.append(c_counts.real(row))
  74 + model.append(c_model.real(row))
  75 + resid.append(c_resid.real(row))
  76 +
  77 + # Calculate count error
  78 + err = math.sqrt(c_counts.real(row))
  79 + if err == 1:
  80 + # this prevents visualization problem in matplotlib
  81 + err = 0.99
  82 + else:
  83 + pass
  84 + e_counts.append(err)
  85 +
  86 + # Calculate residual error
  87 + if algorithm == 'SUB':
  88 + e_resid.append(err)
  89 + elif algorithm == 'SUBDIV':
  90 + e_resid.append(err / c_model.real(row))
  91 + elif algorithm == 'SUBDIVSQRT':
  92 + e_resid.append(err / math.sqrt(c_model.real(row)))
  93 + elif algorithm == 'SIGNIFICANCE':
  94 + e_resid.append(1.)
  95 +
  96 + # Return tuple
  97 + return counts, model, resid, e_counts, e_resid
  98 +
  99 +
  100 +# ============= #
  101 +# Plot spectrum #
  102 +# ============= #
  103 +def plot_residuals(filename, plotfile, hdu):
  104 + """
  105 + Plot spectrum
  106 +
  107 + Parameters
  108 + ----------
  109 + filename : str
  110 + Name of spectrum FITS file
  111 + plotfile : str
  112 + Plot file name
  113 + hdu : int
  114 + Number of observation to plot
  115 + """
  116 + # Read spectrum file
  117 + fits = gammalib.GFits(filename)
  118 + table = fits.table(1 + hdu)
  119 + c_emin = table['Emin']
  120 + c_emax = table['Emax']
  121 + c_counts = table['Counts']
  122 + c_model = table['Model']
  123 + c_resid = table['Residuals']
  124 + try:
  125 + c_counts_off = table['Counts_Off']
  126 + c_model_off = table['Model_Off']
  127 + c_resid_off = table['Residuals_Off']
  128 + is_onoff = True
  129 + except:
  130 + is_onoff = False
  131 +
  132 + # Initialise arrays to be filled
  133 + ebounds = []
  134 + ed_engs = []
  135 + eu_engs = []
  136 + em_engs = []
  137 + counts = []
  138 + e_counts = []
  139 + model = []
  140 + resid = []
  141 + e_resid = []
  142 + counts_off = []
  143 + e_counts_off = []
  144 + model_off = []
  145 + resid_off = []
  146 + e_resid_off = []
  147 +
  148 + # Residual algorithm
  149 + algorithm = table.header()['ALGORITHM'].string()
  150 +
  151 + # Loop over rows of the file
  152 + nrows = table.nrows()
  153 +
  154 + # Add first energy boundary
  155 + ebounds.append(c_emin.real(0))
  156 + for row in range(nrows):
  157 +
  158 + # Boundaries
  159 + ebounds.append(c_emax.real(row))
  160 +
  161 + # Geometrical mean energy and errors
  162 + emean = math.sqrt(c_emin.real(row) * c_emax.real(row))
  163 + em_engs.append(emean)
  164 + ed_engs.append(emean - c_emin.real(row))
  165 + eu_engs.append(c_emax.real(row) - emean)
  166 +
  167 + # counts, model, residuals
  168 + counts, model, resid, e_counts, e_resid = fill_cmr(row, counts, model,
  169 + resid, e_counts,
  170 + e_resid, c_counts,
  171 + c_model, c_resid,
  172 + algorithm)
  173 + # On/off
  174 + if is_onoff:
  175 + counts_off, model_off, resid_off, e_counts_off, e_resid_off = fill_cmr(
  176 + row, counts_off,
  177 + model_off, resid_off,
  178 + e_counts_off,
  179 + e_resid_off,
  180 + c_counts_off,
  181 + c_model_off, c_resid_off,
  182 + algorithm)
  183 +
  184 + # Add model value to be compatible with plt.step
  185 + model = [model[0]] + model
  186 + if is_onoff:
  187 + model_off = [model_off[0]] + model_off
  188 +
  189 + # Plot
  190 + axarr = []
  191 + if is_onoff:
  192 + f = plt.figure(figsize=(8, 4))
  193 + gs = gridspec.GridSpec(2, 2, height_ratios=[2, 1])
  194 + gs.update(hspace=0, left=0.09, right=0.97)
  195 + axarr.append(f.add_subplot(gs[0]))
  196 + axarr.append(f.add_subplot(gs[2]))
  197 + axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
  198 + axarr.append(f.add_subplot(gs[3], sharex=axarr[1]))
  199 + else:
  200 + f = plt.figure()
  201 + gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
  202 + gs.update(hspace=0)
  203 + axarr.append(f.add_subplot(gs[0]))
  204 + axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
  205 + plt.setp(axarr[1].get_xticklabels(), visible=False)
  206 + axarr[0].set_xscale('log')
  207 + axarr[1].set_xscale('log')
  208 + axarr[0].set_yscale('log')
  209 +
  210 + # Counts and model
  211 + axarr[0].errorbar(em_engs, counts, yerr=e_counts, xerr=[ed_engs, eu_engs],
  212 + fmt='ko', capsize=0, linewidth=2, zorder=2, label='Data')
  213 + axarr[0].step(ebounds, model, color='0.5', linewidth=2, zorder=1,
  214 + label='Model')
  215 + axarr[0].set_ylabel('Counts')
  216 +
  217 + # Residuals
  218 + axarr[1].errorbar(em_engs, resid, yerr=e_resid, xerr=[ed_engs, eu_engs],
  219 + fmt='ko', capsize=0, linewidth=2, zorder=2)
  220 + axarr[1].axhline(0, color='0.5', linestyle='--')
  221 + axarr[1].set_xlabel('Energy (TeV)')
  222 + if algorithm == 'SUB':
  223 + axarr[1].set_ylabel('Residuals (counts)')
  224 + elif algorithm == 'SUBDIV':
  225 + axarr[1].set_ylabel('Residuals (fraction)')
  226 + elif algorithm == 'SUBDIVSQRT' or algorithm == 'SIGNIFICANCE':
  227 + axarr[1].set_ylabel(r'Residuals ($\sigma$)')
  228 +
  229 + # On/Off
  230 + if is_onoff:
  231 + axarr[0].set_title('ON')
  232 + axarr[2].set_title('OFF')
  233 + plt.setp(axarr[3].get_xticklabels(), visible=False)
  234 + axarr[2].set_xscale('log')
  235 + axarr[3].set_xscale('log')
  236 + axarr[2].set_yscale('log')
  237 +
  238 + # Counts and model
  239 + axarr[2].errorbar(em_engs, counts_off, yerr=e_counts_off,
  240 + xerr=[ed_engs, eu_engs],
  241 + fmt='ko', capsize=0, linewidth=2, zorder=2,
  242 + label='Data')
  243 + axarr[2].step(ebounds, model_off, color='0.5', linewidth=2, zorder=1,
  244 + label='Model')
  245 +
  246 + # Residuals
  247 + axarr[3].errorbar(em_engs, resid_off, yerr=e_resid_off,
  248 + xerr=[ed_engs, eu_engs],
  249 + fmt='ko', capsize=0, linewidth=2, zorder=2)
  250 + axarr[3].axhline(0, color='0.5', linestyle='--')
  251 + axarr[3].set_xlabel('Energy (TeV)')
  252 +
  253 + # Add spectra of individual components
  254 + skiplist = ['Counts', 'Model', 'Residuals',
  255 + 'Counts_Off', 'Model_Off', 'Residuals_Off',
  256 + 'Emin', 'Emax']
  257 + for s in range(table.ncols()):
  258 + if table[s].name() in skiplist:
  259 + pass
  260 + else:
  261 + component = []
  262 + for row in range(nrows):
  263 + component.append(table[s].real(row))
  264 + component = [component[0]] + component
  265 + axarr[0].step(ebounds, component, zorder=0,
  266 + label=table[s].name())
  267 +
  268 + # Add legend
  269 + axarr[0].legend(loc='best')
  270 +
  271 + # Show figure
  272 + if len(plotfile) > 0:
  273 + plt.savefig(plotfile)
  274 + else:
  275 + plt.show()
  276 +
  277 + # Return
  278 + return
  279 +
  280 +
  281 +# =============== #
  282 +# Show residuals #
  283 +# =============== #
  284 +def show_residuals():
  285 + """
  286 + Show residuals
  287 + """
  288 + # Set usage string
  289 + usage = 'show_spectrum.py [-p plotfile] [-h hdu] [file]'
  290 +
  291 + # Set default options
  292 + options = [{'option': '-p', 'value': ''},
  293 + {'option': '-h', 'value': 0}]
  294 +
  295 + # Get arguments and options from command line arguments
  296 + args, options = cscripts.ioutils.get_args_options(options, usage)
  297 +
  298 + # Extract script parameters from options
  299 + plotfile = options[0]['value']
  300 + hdu = options[1]['value']
  301 +
  302 + # Show residuals
  303 + plot_residuals(args[0], plotfile, hdu)
  304 +
  305 + # Return
  306 + return
  307 +
  308 +
  309 +# ======================== #
  310 +# Main routine entry point #
  311 +# ======================== #
  312 +if __name__ == '__main__':
  313 +
  314 + # Show residuals
  315 + show_residuals()
... ...
src/ctcubemask/ctcubemask.hpp
  View file @47cd665
... ... @@ -33,6 +33,7 @@
33 33 #include "GammaLib.hpp"
34 34 #include "GCTALib.hpp"
35 35 #include "ctobservation.hpp"
  36 +#include "GCTAEventCube.hpp"
36 37  
37 38 /* __Definitions _________________________________________________________ */
38 39 #define CTCUBEMASK_NAME "ctcubemask"
... ... @@ -87,5 +88,4 @@ protected:
87 88 bool m_select_energy; //!< Perform energy selection
88 89 };
89 90  
90   -
91 91 #endif /* CTCUBEMASK_HPP */
... ...
test/test_csresspec.py 0 → 100755
  View file @47cd665
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the csresspec script.
  4 +#
  5 +# Copyright (C) 2017-2018 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 csresspec script #
  28 +# =============================== #
  29 +class Test(test):
  30 + """
  31 + Test class for csresspec script
  32 +
  33 + This test class makes unit tests for the csresspec 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 +
  45 + # Set members
  46 + self._inonoff = self._datadir + '/onoff_obs.xml'
  47 + self._onoff_model = self._datadir + '/onoff_model.xml'
  48 +
  49 + # Return
  50 + return
  51 +
  52 + # Set test functions
  53 + def set(self):
  54 + """
  55 + Set all test functions
  56 + """
  57 + # Set test name
  58 + self.name('csresspec')
  59 +
  60 + # Append tests
  61 + self.append(self._test_cmd, 'Test csresspec on command line')
  62 + self.append(self._test_python, 'Test csresspec from Python')
  63 +
  64 + # Return
  65 + return
  66 +
  67 + # Test csresspec on command line
  68 + def _test_cmd(self):
  69 + """
  70 + Test csresspec on the command line.
  71 + """
  72 + # Set script name
  73 + csresspec = self._script('csresspec')
  74 +
  75 + # Setup csresspec command
  76 + cmd = csresspec + ' inobs="' + self._events + '"' + \
  77 + ' outfile="csresspec_cmd1.fits"' + \
  78 + ' inmodel="' + self._model + '"' + \
  79 + ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
  80 + ' emin=0.1 emax=100.0 enumbins=20' + \
  81 + ' mask=no algorithm="SUB"' + \
  82 + ' logfile="csresspec_cmd1.log" chatter=4'
  83 +
  84 + # Check if execution of wrong command fails
  85 + self.test_assert(self._execute('command_that_does_not_exist') != 0,
  86 + 'Self test of test script')
  87 +
  88 + # Check if execution was successful
  89 + self.test_assert(self._execute(cmd) == 0,
  90 + 'Check successful execution from command line')
  91 +
  92 + # Check result file
  93 + self._check_result_file('csresspec_cmd1.fits', 20, 5)
  94 +
  95 + # Setup csresspec command
  96 + cmd = csresspec + ' inobs="event_file_that_does_not_exist.fits"' + \
  97 + ' outfile="csresspec_cmd2.fits"' + \
  98 + ' inmodel="' + self._model + '"' + \
  99 + ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
  100 + ' emin=0.1 emax=100.0 enumbins=20' + \
  101 + ' mask=no algorithm="SUBDIVSQRT"' + \
  102 + ' logfile="csresspec_cmd2.log" debug=yes chatter=2'
  103 +
  104 + # Check if execution failed
  105 + self.test_assert(self._execute(cmd) != 0,
  106 + 'Check invalid input file when executed from command line')
  107 +
  108 + # Check csresspec --help
  109 + self._check_help(csresspec)
  110 +
  111 + # Return
  112 + return
  113 +
  114 + # Test csresspec from Python
  115 + def _test_python(self):
  116 + """
  117 + Test csresspec from Python
  118 + """
  119 + # Set-up csresspec for event list
  120 + resspec = cscripts.csresspec()
  121 + resspec['inobs'] = self._events
  122 + resspec['outfile'] = 'csresspec_py1.fits'
  123 + resspec['inmodel'] = self._model
  124 + resspec['caldb'] = self._caldb
  125 + resspec['irf'] = self._irf
  126 + resspec['emin'] = 0.1
  127 + resspec['emax'] = 100.0
  128 + resspec['enumbins'] = 20
  129 + resspec['mask'] = False
  130 + resspec['algorithm'] = 'SUB'
  131 + resspec['logfile'] = 'csresspec_py1.log'
  132 + resspec['chatter'] = 2
  133 +
  134 + # Run csresspec script
  135 + resspec.logFileOpen() # Make sure we get a log file
  136 + resspec.run()
  137 + resspec.save()
  138 +
  139 + # Check result file
  140 + self._check_result_file('csresspec_py1.fits', 20, 5)
  141 +
  142 + # Set-up csresspec for event list
  143 + # with computation of individual model components
  144 + resspec = cscripts.csresspec()
  145 + resspec['inobs'] = self._events
  146 + resspec['outfile'] = 'csresspec_py2.fits'
  147 + resspec['inmodel'] = self._model
  148 + resspec['caldb'] = self._caldb
  149 + resspec['irf'] = self._irf
  150 + resspec['emin'] = 0.1
  151 + resspec['emax'] = 100.0
  152 + resspec['enumbins'] = 20
  153 + resspec['mask'] = False
  154 + resspec['algorithm'] = 'SUBDIV'
  155 + resspec['components'] = True
  156 + resspec['logfile'] = 'csresspec_py2.log'
  157 + resspec['chatter'] = 2
  158 +
  159 + # Run csresspec script
  160 + resspec.execute()
  161 +
  162 + # Check result file
  163 + self._check_result_file('csresspec_py2.fits', 20, 7)
  164 +
  165 + # Set-up csresspec for counts cube
  166 + resspec = cscripts.csresspec()
  167 + resspec['inobs'] = self._cntcube
  168 + resspec['modcube'] = 'NONE'
  169 + resspec['expcube'] = 'NONE'
  170 + resspec['psfcube'] = 'NONE'
  171 + resspec['edispcube'] = 'NONE'
  172 + resspec['bkgcube'] = 'NONE'
  173 + resspec['caldb'] = self._caldb
  174 + resspec['irf'] = self._irf
  175 + resspec['inmodel'] = self._model
  176 + resspec['outfile'] = 'csresspec_py3.fits'
  177 + resspec['mask'] = False
  178 + resspec['algorithm'] = 'SUBDIVSQRT'
  179 + resspec['logfile'] = 'csresspec_py3.log'
  180 + resspec['chatter'] = 3
  181 +
  182 + # Run csresspec script
  183 + resspec.execute()
  184 +
  185 + # Check result file
  186 + self._check_result_file('csresspec_py3.fits', 20, 5)
  187 +
  188 + # Set-up csresspec for counts cube
  189 + # with source mask
  190 + resspec = cscripts.csresspec()
  191 + resspec['inobs'] = self._cntcube
  192 + resspec['modcube'] = 'NONE'
  193 + resspec['expcube'] = 'NONE'
  194 + resspec['psfcube'] = 'NONE'
  195 + resspec['edispcube'] = 'NONE'
  196 + resspec['bkgcube'] = 'NONE'
  197 + resspec['caldb'] = self._caldb
  198 + resspec['irf'] = self._irf
  199 + resspec['inmodel'] = self._model
  200 + resspec['outfile'] = 'csresspec_py4.fits'
  201 + resspec['mask'] = True
  202 + resspec['ra'] = 83.63
  203 + resspec['dec'] = 22.01
  204 + resspec['rad'] = 0.2
  205 + resspec['regfile'] = 'NONE'
  206 + resspec['algorithm'] = 'SUBDIVSQRT'
  207 + resspec['logfile'] = 'csresspec_py4.log'
  208 + resspec['chatter'] = 3
  209 +
  210 + # Run csresspec script
  211 + resspec.execute()
  212 +
  213 + # Check result file
  214 + self._check_result_file('csresspec_py4.fits', 20, 5)
  215 +
  216 + # Set-up csresspec for On/Off observations
  217 + resspec = cscripts.csresspec()
  218 + resspec['inobs'] = self._inonoff
  219 + resspec['inmodel'] = self._onoff_model
  220 + resspec['outfile'] = 'csresspec_py5.fits'
  221 + resspec['algorithm'] = 'SIGNIFICANCE'
  222 + resspec['stack'] = False
  223 + resspec['logfile'] = 'csresspec_py5.log'
  224 + resspec['chatter'] = 3
  225 +
  226 + # Run csresspec script
  227 + resspec.execute()
  228 +
  229 + # Check result file
  230 + self._check_result_file('csresspec_py5.fits', 20, 8, nobs=2)
  231 +
  232 + # Set-up csresspec for On/Off observations
  233 + # in stacked mode and with WSTAT background
  234 + resspec = cscripts.csresspec()
  235 + resspec['inobs'] = self._inonoff
  236 + resspec['inmodel'] = self._onoff_model
  237 + resspec['outfile'] = 'csresspec_py6.fits'
  238 + resspec['algorithm'] = 'SIGNIFICANCE'
  239 + resspec['statistic'] = 'WSTAT'
  240 + resspec['stack'] = True
  241 + resspec['logfile'] = 'csresspec_py6.log'
  242 + resspec['chatter'] = 3
  243 +
  244 + # Run csresspec script
  245 + resspec.execute()
  246 +
  247 + # Check result file
  248 + self._check_result_file('csresspec_py6.fits', 20, 8)
  249 +
  250 + # Return
  251 + return
  252 +
  253 + # Check result file
  254 + def _check_result_file(self, filename, enumbins, ncols, nobs=1):
  255 + """
  256 + Check result file
  257 + """
  258 + # Load residual file
  259 + fits = gammalib.GFits(filename)
  260 +
  261 + # Check number of HDUs = observations (+ first empty image)
  262 + self.test_value(fits.size(), 1 + nobs, 'Number of observations')
  263 +
  264 + # Check first hdu with residuals
  265 + table = fits[1]
  266 +
  267 + # Check number of rows = number of energy bins
  268 + self.test_value(table.nrows(), enumbins,
  269 + 'Check for %d rows in residuals' % enumbins)
  270 +
  271 + # Check number of columns, 2 for energy bounds + 3 for 3D or 6 for
  272 + # On/Off + 1 for each model component if components = yes
  273 + self.test_value(table.ncols(), ncols,
  274 + 'Check for %d columns in spectrum' % ncols)
  275 +
  276 + # Return
  277 + return
... ...
test/test_python_cscripts.py
  View file @47cd665
... ... @@ -41,6 +41,7 @@ import test_csphagen
41 41 import test_csphasecrv
42 42 import test_cspull
43 43 import test_csresmap
  44 +import test_csresspec
44 45 import test_csroot2caldb
45 46 import test_cssens
46 47 import test_csspec
... ... @@ -133,6 +134,7 @@ def test(installed=False, debug=False):
133 134 test_csphasecrv.Test(),
134 135 test_cspull.Test(),
135 136 test_csresmap.Test(),
  137 + test_csresspec.Test(),
136 138 test_csroot2caldb.Test(),
137 139 test_cssens.Test(),
138 140 test_csspec.Test(),
... ...