Maria Benito / ctools

Commit 6dfd7517a60b8d2cec45243ce92e8427a816ba9a

Authored by Jürgen Knödlseder
2 parents 6d6e7f25 6049bf07

Merge branch 'devel' into 1886-dc-tutorial

Inline Side-by-side
Showing 18 changed files with 1290 additions and 259 deletions
ChangeLog
  View file @6dfd751
1   -2016-11-25
  1 +2016-12-01
2 2  
3 3 * Version 1.2.0 released
4 4 ========================
5 5  
  6 + Add cssrcdetect script (#1617)
6 7 Add background subtraction to ctskymap tool (#1882)
7 8 Add csobsselect script (#1881)
8 9 Add csmodelselect script (#1879)
... ...
1 1 New Features and Important Changes in ctools 1.2.0
2 2  
3   -25 November 2016
  3 +1 December 2016
4 4  
5 5  
6 6 Introduction
... ... @@ -206,30 +206,33 @@ The class now derives from ctlikelihood (#1828).
206 206  
207 207 cscripts
208 208 --------
209   -Added the support modules "modutils" and "ioutils" to the cscripts module and
210   -removed the obsolete functions fit(), cterror(), cntmap() and modmap() from the
  209 +Add the support modules "modutils" and "ioutils" to the cscripts module and
  210 +remove the obsolete functions fit(), cterror(), cntmap() and modmap() from the
211 211 "obsutils" module (#1830).
212 212  
213 213 The csroot2caldb script now takes the zenith and azimuth angles as parameters
214 214 and writes them into the calibration database (#1835).
215 215  
216   -Added a csviscube script that computes the visibility cube providing the number
  216 +Add csviscube script that computes the visibility cube providing the number
217 217 of hours that a given position on the sky is observable under a given zenith
218 218 angle by an IACT array during a given time period. The tool implements so far
219 219 a Sun constraint, but no Moon constraint (#1835).
220 220  
221   -A publish() method has been added to cslightcrv and csspec that publishes the
222   -results on a VO hub (#1871).
  221 +Add publish() method to cslightcrv and csspec that publishes the results on a
  222 +VO hub (#1871).
223 223  
224 224 Fix energy boundary handling in cssens that has been corrupted during code
225 225 refactoring (#1873).
226 226  
227   -Added csmodelselect script that allows selecting models from a model definition
  227 +Add csmodelselect script that allows selecting models from a model definition
228 228 XML file (#1879).
229 229  
230   -Added csobsselect script that allows selecting observations from an observation
  230 +Add csobsselect script that allows selecting observations from an observation
231 231 definition XML file (#1881).
232 232  
  233 +Add cssrcdetect script to detect sources in a sky map (#1617).
  234 +
  235 +
233 236  
234 237 Examples
235 238 --------
... ...
cscripts/Makefile.am
  View file @6dfd751
... ... @@ -36,6 +36,7 @@ pythonscripts = $(builddir)/__init__.py \
36 36 $(srcdir)/csroot2caldb.py \
37 37 $(srcdir)/cssens.py \
38 38 $(srcdir)/csspec.py \
  39 + $(srcdir)/cssrcdetect.py \
39 40 $(srcdir)/cstsdist.py \
40 41 $(srcdir)/cstsmapmerge.py \
41 42 $(srcdir)/cstsmapsplit.py \
... ... @@ -65,6 +66,7 @@ parfiles = $(srcdir)/cscaldb.par \
65 66 $(srcdir)/csroot2caldb.par \
66 67 $(srcdir)/cssens.par \
67 68 $(srcdir)/csspec.par \
  69 + $(srcdir)/cssrcdetect.par \
68 70 $(srcdir)/cstsdist.par \
69 71 $(srcdir)/cstsmapmerge.par \
70 72 $(srcdir)/cstsmapsplit.par \
... ... @@ -91,6 +93,7 @@ linkscripts = cscaldb \
91 93 csroot2caldb \
92 94 cssens \
93 95 csspec \
  96 + cssrcdetect \
94 97 cstsdist \
95 98 cstsmapmerge \
96 99 cstsmapsplit \
... ... @@ -119,6 +122,7 @@ testscripts = $(srcdir)/tests/__init__.py \
119 122 $(top_srcdir)/test/test_csroot2caldb.py \
120 123 $(top_srcdir)/test/test_cssens.py \
121 124 $(top_srcdir)/test/test_csspec.py \
  125 + $(top_srcdir)/test/test_cssrcdetect.py \
122 126 $(top_srcdir)/test/test_cstsdist.py \
123 127 $(top_srcdir)/test/test_cstsmapmerge.py \
124 128 $(top_srcdir)/test/test_cstsmapsplit.py \
... ... @@ -160,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
160 164 $(top_srcdir)/test/data/pull.dat \
161 165 $(top_srcdir)/test/data/sensitivity.dat \
162 166 $(top_srcdir)/test/data/sensitivity_integral.dat \
  167 + $(top_srcdir)/test/data/skymap.fits \
163 168 $(top_srcdir)/test/data/spectrum.fits \
164 169 $(top_srcdir)/test/data/ts.dat \
165 170 $(top_srcdir)/test/data/models_spatial.xml \
... ...
cscripts/__init__.py.in
  View file @6dfd751
... ... @@ -39,6 +39,7 @@ __all__ = [
39 39 "csroot2caldb",
40 40 "cssens",
41 41 "csspec",
  42 + "cssrcdetect",
42 43 "cstsdist",
43 44 "cstsmapmerge",
44 45 "cstsmapsplit",
... ... @@ -71,6 +72,7 @@ from cscripts.csresmap import csresmap
71 72 from cscripts.csroot2caldb import csroot2caldb
72 73 from cscripts.cssens import cssens
73 74 from cscripts.csspec import csspec
  75 +from cscripts.cssrcdetect import cssrcdetect
74 76 from cscripts.cstsdist import cstsdist
75 77 from cscripts.cstsmapmerge import cstsmapmerge
76 78 from cscripts.cstsmapsplit import cstsmapsplit
... ...
cscripts/csmodelinfo.py
  View file @6dfd751
... ... @@ -21,6 +21,7 @@
21 21 import sys
22 22 import gammalib
23 23 import ctools
  24 +from cscripts import modutils
24 25  
25 26  
26 27 # ================= #
... ... @@ -86,155 +87,6 @@ class csmodelinfo(ctools.cscript):
86 87 # Return
87 88 return
88 89  
89   - def _ds9attributes(self, model):
90   - """
91   - Returns attributes for DS9 string
92   - Parameters
93   - ----------
94   - model : `~gammalib.GModel`
95   - Model
96   -
97   - Returns
98   - -------
99   - attributes : str
100   - DS9 attribute string
101   - """
102   - # Determine region color. A model where all parameters are fixed will
103   - # get the color defined by "fixed_color", a model with at least one
104   - # free parameter will get the color defined by "free_color".
105   - color = self['fixed_color'].string()
106   - for par in model:
107   - if par.is_free():
108   - color = self['free_color'].string()
109   - break
110   -
111   - # Set DS9 attributes
112   - attributes = (' color=%s width=%d font=%s %d %s %s' %
113   - (color, self['width'].integer(),
114   - self['fontfamily'].string(), self['fontsize'].integer(),
115   - self['fontweight'].string(), self['fontslant'].string()))
116   -
117   - # Return attributes
118   - return attributes
119   -
120   - def _model2ds9string(self, model):
121   - """
122   - Converts model into a DS9 region string
123   -
124   - Parameters
125   - ----------
126   - model : `~gammalib.GModel`
127   - Model
128   -
129   - Returns
130   - -------
131   - ds9string : str
132   - DS9 region string
133   - """
134   - # Initialise DS9 region string
135   - ds9string = ''
136   -
137   - # Retrieve model sky direction. The model is skipped in case that
138   - # it does not provide a sky direction.
139   - is_valid = True
140   - try:
141   - modelpos = model.spatial().dir()
142   - except AttributeError:
143   - msg = ('Skip model "%s" since it has no sky direction.\n' %
144   - model.name())
145   - self._log(msg)
146   - is_valid = False
147   -
148   - # Continue only if sky direction was found
149   - if is_valid:
150   -
151   - # Retrieve model type and name
152   - modeltype = model.type()
153   - modelname = model.name()
154   -
155   - # Handle point source
156   - if modeltype == 'PointSource':
157   -
158   - # Append point with Right Ascension and Declination to the DS9
159   - # string. The type of the point is specified by the "pnt_type"
160   - # parameter, the size of the point by the "pnt_mark_size"
161   - # parameter
162   - ds9string += ('point(%.6f,%.6f) # point=%s %d' %
163   - (modelpos.ra_deg(), modelpos.dec_deg(),
164   - self['pnt_type'].string(),
165   - self['pnt_mark_size'].integer()))
166   -
167   - # Handle extended sources
168   - elif modeltype == "ExtendedSource":
169   -
170   - # Retrieve spatial model
171   - spatial = model.spatial()
172   - classname = spatial.classname()
173   -
174   - # Handle radial sources
175   - if 'Radial' in classname:
176   -
177   - # Retrieve short name of model class (e.g. "Gauss"
178   - # or "Disk)
179   - shorttype = classname.split('Radial')[-1]
180   -
181   - # Handle Disk and Shell model
182   - if (classname == 'GModelSpatialRadialDisk' or
183   - classname == 'GModelSpatialRadialShell'):
184   - size = spatial.radius()
185   -
186   - # Handle Gauss Model
187   - elif classname == 'GModelSpatialRadialGauss':
188   - size = spatial.sigma()
189   -
190   - # Skip if source is unknown
191   - else:
192   - msg = ('Skip model "%s" since the radial model "%s" '
193   - 'is unknown.\n' % (model.name(), classname))
194   - self._log(msg)
195   - is_valid = False
196   -
197   - # Append circle to DS9 string
198   - if is_valid:
199   - ds9string += ('circle(%.6f,%.6f,%.6f) #' %
200   - (modelpos.ra_deg(), modelpos.dec_deg(),
201   - size*3600.0))
202   -
203   - # Handle elliptical sources
204   - elif 'Elliptical' in classname:
205   -
206   - # Retrieve short name and source size
207   - shorttype = classname.split('Elliptical')[-1]
208   - size1 = spatial.semimajor()
209   - size2 = spatial.semiminor()
210   - angle = spatial.posangle()
211   -
212   - # Append ellipse to DS9 string
213   - ds9string += ('ellipse(%.6f,%.6f,%.6f,%.6f,%.6f) #' %
214   - (modelpos.ra_deg(), modelpos.dec_deg(),
215   - size1*3600.0, size2*3600.0, angle+90.0))
216   -
217   - # Skip if source is neither radial nor elliptical
218   - else:
219   - msg = ('Skip model "%s" since the model "%s" is neither '
220   - 'a point source, a radial source, nor an elliptical '
221   - 'source.\n' % (model.name(), classname))
222   - self._log(msg)
223   - is_valid = False
224   -
225   - # Add short model type to modelname
226   - if self['show_ext_type'].boolean():
227   - modelname +=' ('+shorttype+')'
228   -
229   - # Add DS9 attributes
230   - if is_valid:
231   - ds9string += self._ds9attributes(model)
232   - if self['show_labels'].boolean():
233   - ds9string += ' text={'+modelname+'}'
234   -
235   - # Return string
236   - return ds9string
237   -
238 90  
239 91 # Public methods
240 92 def run(self):
... ... @@ -407,36 +259,34 @@ class csmodelinfo(ctools.cscript):
407 259  
408 260 # Log filename
409 261 self._log_value(gammalib.NORMAL, 'DS9 filename', ds9file.url())
410   -
411   - # Open file
412   - f = open(ds9file.url(), 'w')
413   -
414   - # Write coordinate system
415   - f.write('fk5\n')
416   -
417   - # Loop over models
418   - for model in self._models:
419   -
420   - # Continue only if point source or extended source model
421   - if (model.type() == 'PointSource' or
422   - model.type() == 'ExtendedSource'):
423   - line = self._model2ds9string(model)
424   - if len(line):
425   - f.write(line+'\n')
426   -
427   - # Logging for diffuse components
428   - elif model.type() == 'DiffuseSource':
429   - self._log('Skipping diffuse model "'+model.name()+'".\n')
430   -
431   - # Logging for background components
432   - else:
433   - if self._logExplicit():
434   - self._log('Skipping background model "'+model.name()+
435   - '".\n')
436   -
437   - # Close file
438   - f.close()
439   -
  262 +
  263 + # Get ds9 parameters
  264 + pnt_type = self['pnt_type'].string()
  265 + pnt_mark_size = self['pnt_mark_size'].integer()
  266 + free_color = self['free_color'].string()
  267 + fixed_color = self['fixed_color'].string()
  268 + width = self['width'].integer()
  269 + fontfamily = self['fontfamily'].string()
  270 + fontsize = self['fontsize'].integer()
  271 + fontweight = self['fontweight'].string()
  272 + fontslant = self['fontslant'].string()
  273 +
  274 + # Save DS9 region file
  275 + errors = modutils.models2ds9file(self._models, ds9file.url(),
  276 + pnt_type=pnt_type,
  277 + pnt_mark_size=pnt_mark_size,
  278 + free_color=free_color,
  279 + fixed_color=fixed_color,
  280 + width=width,
  281 + fontfamily=fontfamily,
  282 + fontsize=fontsize,
  283 + fontweight=fontweight,
  284 + fontslant=fontslant)
  285 +
  286 + # Log errors
  287 + if len(errors):
  288 + self._log(errors)
  289 +
440 290 # Return
441 291 return
442 292  
... ...
cscripts/cssrcdetect.par 0 → 100644
  View file @6dfd751
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: cssrcdetect.par #
  8 +# #
  9 +# Author: Juergen Knoedlseder #
  10 +# jknodlseder@irap.omp.eu #
  11 +# IRAP, Toulouse, France #
  12 +# #
  13 +# Purpose: Parameter file for the cssrcdetect script #
  14 +# #
  15 +#############################################################################
  16 +#
  17 +# Input and output parameters
  18 +#============================
  19 +inmap, f, a, skymap.fits,,, "Input sky map file"
  20 +outmodel, f, a, models.xml,,, "Output model definition XML file"
  21 +outds9file, f, a, ds9.reg,,, "Output DS9 region file"
  22 +
  23 +#
  24 +# Script parameters
  25 +#==================
  26 +srcmodel, s, a, POINT,POINT,, "Source model type"
  27 +bkgmodel, s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type"
  28 +threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)"
  29 +maxsrcs, i, h, 20,1,1000, "Maximum number of sources"
  30 +exclrad, r, h, 0.2,,, "Exclusion radius (deg)"
  31 +fit_pos, b, h, yes,,, "Fit source positions?"
  32 +fit_shape, b, h, yes,,, "Fit source shapes?"
  33 +
  34 +#
  35 +# Standard parameters
  36 +#====================
  37 +chatter, i, h, 2,0,4, "Chattiness of output"
  38 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  39 +debug, b, h, no,,, "Debugging mode activated"
  40 +mode, s, h, ql,,, "Mode of automatic parameters"
  41 +logfile, f, h, cssrcdetect.log,,, "Log filename"
... ...
cscripts/cssrcdetect.py 0 → 100755
  View file @6dfd751
  1 +#!/usr/bin/env python
  2 +# ==========================================================================
  3 +# Detects sources in a sky map
  4 +#
  5 +# Copyright (C) 2016 Juergen Knoedlseder
  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 ctools
  25 +from cscripts import modutils
  26 +
  27 +
  28 +# ================= #
  29 +# cssrcdetect class #
  30 +# ================= #
  31 +class cssrcdetect(ctools.cscript):
  32 + """
  33 + Detects sources in a sky map
  34 + """
  35 +
  36 + # Constructor
  37 + def __init__(self, *argv):
  38 + """
  39 + Constructor
  40 + """
  41 + # Set name
  42 + self._name = 'cssrcdetect'
  43 + self._version = '1.2.0'
  44 +
  45 + # Set protected members
  46 + self._models = gammalib.GModels()
  47 +
  48 + # Initialise application by calling the appropriate class constructor
  49 + self._init_cscript(argv)
  50 +
  51 + # Return
  52 + return
  53 +
  54 +
  55 + # Private methods
  56 + def _get_parameters(self):
  57 + """
  58 + Get parameters from parfile
  59 + """
  60 + # Query input parameters
  61 + self['inmap'].filename()
  62 +
  63 + # Query further parameters
  64 + self['srcmodel'].string()
  65 + self['bkgmodel'].string()
  66 + self['threshold'].real()
  67 + self['maxsrcs'].integer()
  68 + self['exclrad'].real()
  69 + self['fit_pos'].boolean()
  70 + self['fit_shape'].boolean()
  71 +
  72 + # Query ahead output model filename
  73 + if self._read_ahead():
  74 + self['outmodel'].filename()
  75 + self['outds9file'].filename()
  76 +
  77 + # Write input parameters into logger
  78 + self._log_parameters(gammalib.TERSE)
  79 +
  80 + # Return
  81 + return
  82 +
  83 + def _detect_sources(self, counts):
  84 + """
  85 + Detect sources in counts map
  86 +
  87 + Parameters
  88 + ----------
  89 + counts : `~gammalib.GSkyMap()`
  90 + Counts map
  91 + """
  92 + # Source detection iterations
  93 + for i in range(self['maxsrcs'].integer()):
  94 +
  95 + # Write header
  96 + self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1))
  97 +
  98 + # Get map moments
  99 + mean, std = self._map_moments(counts)
  100 +
  101 + # Log map moments
  102 + self._log_value(gammalib.NORMAL, 'Map mean', mean)
  103 + self._log_value(gammalib.NORMAL, 'Map standard deviation', std)
  104 +
  105 + # Compute threshold
  106 + threshold = mean + self['threshold'].real() * std
  107 +
  108 + # Log map threshold
  109 + self._log_value(gammalib.NORMAL, 'Map threshold', threshold)
  110 +
  111 + # Get maximum value and corresponding sky direction
  112 + value, dir = self._find_maximum(counts, threshold=threshold)
  113 +
  114 + # If maximum found then log maximum and add model
  115 + if dir is not None:
  116 +
  117 + # Set source name
  118 + name = 'Src%3.3d' % (i+1)
  119 +
  120 + # Log maximum
  121 + self._log_value(gammalib.NORMAL, 'Map maximum', str(value))
  122 + self._log_value(gammalib.NORMAL, name+' position', str(dir))
  123 +
  124 + # Add model
  125 + self._add_model(dir, name)
  126 +
  127 + # Remove maximum from map
  128 + counts = self._remove_maximum(counts, dir, mean,
  129 + radius=self['exclrad'].real())
  130 +
  131 + # ... otherwise log that no maximum was found and break iterations
  132 + else:
  133 +
  134 + # Log than no maximum was found
  135 + self._log_value(gammalib.NORMAL, 'Map maximum',
  136 + 'None above threshold')
  137 +
  138 + # Break
  139 + break
  140 +
  141 + # Return
  142 + return
  143 +
  144 + def _find_maximum(self, map, threshold=0.0):
  145 + """
  146 + Find maximum in a sky map
  147 +
  148 + Parameters
  149 + ----------
  150 + map : `~gammalib.GSkyMap()`
  151 + Sky map
  152 + threshold : float, optional
  153 + Threshold for maximum value
  154 +
  155 + Returns
  156 + -------
  157 + value, dir : tuple of float and `~gammalib.GSkyDir()`
  158 + Maximum sky map value and corresponding sky direction
  159 + """
  160 + # Initialise maximum pixel value and sky direction
  161 + value = threshold
  162 + dir = None
  163 +
  164 + # Loop over all pixels and find maximum
  165 + for i in range(map.npix()):
  166 + if map[i] > value:
  167 + value = map[i]
  168 + dir = map.inx2dir(i)
  169 +
  170 + # Return sky direction of maximum
  171 + return value, dir
  172 +
  173 + def _remove_maximum(self, map, dir, value=0.0, radius=0.1):
  174 + """
  175 + Remove maximum from sky map by replacing pixels with a given value
  176 +
  177 + Parameters
  178 + ----------
  179 + map : `~gammalib.GSkyMap()`
  180 + Sky map
  181 + dir : `~gammalib.GSkyDir()`
  182 + Sky direction of maximum
  183 + value : float, optional
  184 + Replacement value
  185 + radius : float, optional
  186 + Radius within which pixel values are replaced
  187 +
  188 + Returns
  189 + -------
  190 + map : `~gammalib.GSkyMap()`
  191 + Sky map with maximum removed
  192 + """
  193 + # Copy skymap
  194 + map_copy = map.copy()
  195 +
  196 + # Loop over all pixels and find maximum
  197 + for i in range(map_copy.npix()):
  198 + mapdir = map_copy.inx2dir(i)
  199 + if mapdir.dist_deg(dir) < radius:
  200 + map_copy[i] = value
  201 +
  202 + # Return copy of map
  203 + return map_copy
  204 +
  205 + def _map_moments(self, map):
  206 + """
  207 + Determine moments of sky map pixels
  208 +
  209 + Parameters
  210 + ----------
  211 + map : `~gammalib.GSkyMap()`
  212 + Sky map
  213 +
  214 + Returns
  215 + -------
  216 + mean, std : tuple of float
  217 + Mean and standard deviation of pixel values
  218 + """
  219 + # Initialise mean and standard deviation
  220 + mean = 0.0
  221 + std = 0.0
  222 +
  223 + # Loop over all pixels
  224 + for i in range(map.npix()):
  225 + mean += map[i]
  226 + std += map[i]*map[i]
  227 +
  228 + # Compute mean and standard deviation
  229 + mean /= float(map.npix())
  230 + std /= float(map.npix())
  231 + std -= mean
  232 + std = math.sqrt(std)
  233 +
  234 + # Return mean and standard deviation
  235 + return mean, std
  236 +
  237 + def _add_model(self, dir, name):
  238 + """
  239 + Add model to model container
  240 +
  241 + Parameters
  242 + ----------
  243 + dir : `~gammalib.GSkyDir()`
  244 + Sky direction of model
  245 + name : str
  246 + Model name
  247 + """
  248 + # Set point source model
  249 + model = self._set_ptsrc(dir)
  250 +
  251 + # Set model name
  252 + model.name(name)
  253 +
  254 + # Append model to container
  255 + self._models.append(model)
  256 +
  257 + # Return
  258 + return
  259 +
  260 + def _set_bkg(self, type):
  261 + """
  262 + Set background model
  263 +
  264 + Parameters
  265 + ----------
  266 + type : str
  267 + Model type ('IRF', 'AEFF' or 'CUBE')
  268 +
  269 + Returns
  270 + -------
  271 + model : `~gammalib.GModelData()`
  272 + Background model
  273 + """
  274 + # Set spectral component
  275 + spectral = gammalib.GModelSpectralPlaw(1.0, 0.0,
  276 + gammalib.GEnergy(1.0, 'TeV'))
  277 +
  278 + # Set background model
  279 + if type == 'IRF':
  280 + model = gammalib.GCTAModelIrfBackground(spectral)
  281 + elif type == 'AEFF':
  282 + model = gammalib.GCTAModelAeffBackground(spectral)
  283 + elif type == 'CUBE':
  284 + model = gammalib.GCTAModelCubeBackground(spectral)
  285 + else:
  286 + model = None
  287 +
  288 + # Set background name
  289 + if model is not None:
  290 + model.name('Background')
  291 +
  292 + # Return model
  293 + return model
  294 +
  295 + def _set_ptsrc(self, dir):
  296 + """
  297 + Set point source model
  298 +
  299 + Parameters
  300 + ----------
  301 + dir : `~gammalib.GSkyDir()`
  302 + Sky direction of model
  303 +
  304 + Returns
  305 + -------
  306 + model : `~gammalib.GModelSky()`
  307 + Point source model
  308 + """
  309 + # Set spatial component
  310 + spatial = gammalib.GModelSpatialPointSource(dir)
  311 +
  312 + # Get fit parameters
  313 + fit_pos = self['fit_pos'].boolean()
  314 + fit_shape = self['fit_shape'].boolean()
  315 +
  316 + # Loop over all parameters
  317 + for par in spatial:
  318 +
  319 + # Handle position parameters
  320 + if par.name() == 'RA' or par.name() == 'DEC':
  321 + if fit_pos:
  322 + par.free()
  323 + else:
  324 + par.fix()
  325 +
  326 + # Set spectral component (powerlaw with 1% of Crab)
  327 + spectral = gammalib.GModelSpectralPlaw(5.7e-18, -2.48,
  328 + gammalib.GEnergy(0.3, 'TeV'))
  329 +
  330 + # Set sky model
  331 + model = gammalib.GModelSky(spatial, spectral)
  332 +
  333 + # Return model
  334 + return model
  335 +
  336 +
  337 + # Public methods
  338 + def run(self):
  339 + """
  340 + Run the script
  341 + """
  342 + # Clear model container
  343 + self._models.clear()
  344 +
  345 + # Switch screen logging on in debug mode
  346 + if self._logDebug():
  347 + self._log.cout(True)
  348 +
  349 + # Get parameters
  350 + self._get_parameters()
  351 +
  352 + # Write header
  353 + self._log_header1(gammalib.NORMAL, 'Source detection')
  354 +
  355 + # Get skymap filename
  356 + inmap = self['inmap'].filename()
  357 +
  358 + # Open sky map file
  359 + fits = gammalib.GFits(inmap)
  360 +
  361 + # Extract primary extension as counts map
  362 + counts = gammalib.GSkyMap(fits.image(0))
  363 +
  364 + # Detect sources
  365 + self._detect_sources(counts)
  366 +
  367 + # Close sky map file
  368 + fits.close()
  369 +
  370 + # Write detected sources as models into logger
  371 + self._log_models(gammalib.NORMAL, self._models, 'Detected source model')
  372 +
  373 + # Write header
  374 + self._log_header1(gammalib.NORMAL, 'Add background model')
  375 +
  376 + # Add background model if requested
  377 + bkgmodel = self._set_bkg(self['bkgmodel'].string())
  378 + if bkgmodel is not None:
  379 +
  380 + # Append model
  381 + self._models.append(bkgmodel)
  382 +
  383 + # Log model
  384 + self._log_string(gammalib.NORMAL, str(bkgmodel))
  385 +
  386 + # ... otherwise notify that no background model is added
  387 + else:
  388 + self._log_value(gammalib.NORMAL, 'Background model', 'None')
  389 +
  390 + # Return
  391 + return
  392 +
  393 + def save(self):
  394 + """
  395 + Save sources into model definition XML file and ds9 file
  396 + """
  397 + # Write header
  398 + self._log_header1(gammalib.TERSE, 'Save sources')
  399 +
  400 + # Get output filenames
  401 + outmodel = self['outmodel'].filename()
  402 + outds9file = self['outds9file'].filename()
  403 +
  404 + # If model definition XML file exists and clobber flag is false then
  405 + # raise an exception
  406 + if outmodel.exists() and not self._clobber:
  407 + msg = ('Cannot save "'+outmodel.url()+'": File already exists. '
  408 + 'Use parameter clobber=yes to allow overwriting of files.')
  409 + raise RuntimeError(msg)
  410 +
  411 + # If ds9 file exists and clobber flag is false then raise an exception
  412 + elif outds9file.exists() and not self._clobber:
  413 + msg = ('Cannot save "'+outds9file.url()+'": File already exists. '
  414 + 'Use parameter clobber=yes to allow overwriting of files.')
  415 + raise RuntimeError(msg)
  416 +
  417 + # Otherwise log filename and save file
  418 + else:
  419 + # Log model definition XML filename
  420 + self._log_value(gammalib.NORMAL, 'Model definition XML file',
  421 + outmodel.url())
  422 +
  423 + # Save models
  424 + self._models.save(outmodel)
  425 +
  426 + # Log ds9 filename
  427 + self._log_value(gammalib.NORMAL, 'DS9 region file',
  428 + outds9file.url())
  429 +
  430 + # Save models as ds9 region file
  431 + modutils.models2ds9file(self._models, outds9file.url())
  432 +
  433 + # Return
  434 + return
  435 +
  436 + def execute(self):
  437 + """
  438 + Execute the script
  439 + """
  440 + # Open logfile
  441 + self.logFileOpen()
  442 +
  443 + # Read ahead output parameters
  444 + self._read_ahead(True)
  445 +
  446 + # Run the script
  447 + self.run()
  448 +
  449 + # Save observation definition file
  450 + self.save()
  451 +
  452 + # Return
  453 + return
  454 +
  455 +
  456 +# ======================== #
  457 +# Main routine entry point #
  458 +# ======================== #
  459 +if __name__ == '__main__':
  460 +
  461 + # Create instance of application
  462 + app = cssrcdetect(sys.argv)
  463 +
  464 + # Execute application
  465 + app.execute()
... ...
cscripts/modutils.py
  View file @6dfd751
... ... @@ -99,3 +99,298 @@ def test_source(models, srcname, ra=None, dec=None, fitspat=False,
99 99 # Return model container
100 100 return outmodels
101 101  
  102 +
  103 +# ================================= #
  104 +# Returns attributes for ds9 string #
  105 +# ================================= #
  106 +def ds9attributes(model, free_color='green', fixed_color='magenta',
  107 + width=2, fontfamily='helvetica', fontsize=12,
  108 + fontweight='normal', fontslant='roman'):
  109 + """
  110 + Returns attributes for DS9 string
  111 +
  112 + Parameters
  113 + ----------
  114 + model : `~gammalib.GModel`
  115 + Model
  116 + free_color : str, optional
  117 + Color for sources with free parameters (any DS9 color or hex code)
  118 + fixed_color : str, optional
  119 + Color for source without free parameters (any DS9 color or hex code)
  120 + width : int, optional
  121 + Line width for regions
  122 + fontfamily : str, optional
  123 + Font family for source labels (helvetica,times,courier)
  124 + fontsize : int, optional
  125 + Font size for source labels
  126 + fontweight : str, optional
  127 + Font weight for source labels (normal,bold)
  128 + fontslant : str, optional
  129 + Font slant for source labels (roman,italic)
  130 +
  131 + Returns
  132 + -------
  133 + attributes : str
  134 + DS9 attribute string
  135 + """
  136 + # Determine region color. A model where all parameters are fixed will
  137 + # get the color defined by "fixed_color", a model with at least one
  138 + # free parameter will get the color defined by "free_color".
  139 + color = fixed_color
  140 + for par in model:
  141 + if par.is_free():
  142 + color = free_color
  143 + break
  144 +
  145 + # Set DS9 attributes
  146 + attributes = (' color=%s width=%d font="%s %d %s %s"' %
  147 + (color, width, fontfamily, fontsize, fontweight, fontslant))
  148 +
  149 + # Return attributes
  150 + return attributes
  151 +
  152 +
  153 +# =========================== #
  154 +# Convert model to ds9 string #
  155 +# =========================== #
  156 +def model2ds9string(model, pnt_type='cross', pnt_mark_size=12,
  157 + show_labels=True, show_ext_type=True,
  158 + free_color='green', fixed_color='magenta',
  159 + width=2, fontfamily='helvetica', fontsize=12,
  160 + fontweight='normal', fontslant='roman'):
  161 + """
  162 + Converts model into a DS9 region string
  163 +
  164 + Parameters
  165 + ----------
  166 + model : `~gammalib.GModel`
  167 + Model
  168 + pnt_type : str, optional
  169 + Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle)
  170 + pnt_mark_size : integer, optional
  171 + Marker size for point sources
  172 + show_labels : boolean, optional
  173 + Add source labels?
  174 + show_ext_type : boolean, optional
  175 + Show type of extended model in source name?
  176 + free_color : str, optional
  177 + Color for sources with free parameters (any DS9 color or hex code)
  178 + fixed_color : str, optional
  179 + Color for source without free parameters (any DS9 color or hex code)
  180 + width : int, optional
  181 + Line width for regions
  182 + fontfamily : str, optional
  183 + Font family for source labels (helvetica,times,courier)
  184 + fontsize : int, optional
  185 + Font size for source labels
  186 + fontweight : str, optional
  187 + Font weight for source labels (normal,bold)
  188 + fontslant : str, optional
  189 + Font slant for source labels (roman,italic)
  190 +
  191 + Returns
  192 + -------
  193 + ds9string : str
  194 + DS9 region string
  195 + msg : str
  196 + Error message
  197 + """
  198 + # Initialise DS9 region string and message string
  199 + ds9string = ''
  200 + msg = ''
  201 +
  202 + # Retrieve model sky direction. The model is skipped in case that
  203 + # it does not provide a sky direction.
  204 + is_valid = True
  205 + try:
  206 + modelpos = model.spatial().dir()
  207 + except AttributeError:
  208 + msg = ('Skip model "%s" since it has no sky direction.\n' %
  209 + model.name())
  210 + is_valid = False
  211 +
  212 + # Continue only if sky direction was found
  213 + if is_valid:
  214 +
  215 + # Retrieve model type and name
  216 + modeltype = model.type()
  217 + modelname = model.name()
  218 +
  219 + # Handle point source
  220 + if modeltype == 'PointSource':
  221 +
  222 + # Append point with Right Ascension and Declination to the DS9
  223 + # string. The type of the point is specified by the "pnt_type"
  224 + # parameter, the size of the point by the "pnt_mark_size"
  225 + # parameter
  226 + ds9string += ('point(%.6f,%.6f) # point=%s %d' %
  227 + (modelpos.ra_deg(), modelpos.dec_deg(),
  228 + pnt_type, pnt_mark_size))
  229 +
  230 + # Handle extended sources
  231 + elif modeltype == "ExtendedSource":
  232 +
  233 + # Retrieve spatial model
  234 + spatial = model.spatial()
  235 + classname = spatial.classname()
  236 +
  237 + # Handle radial sources
  238 + if 'Radial' in classname:
  239 +
  240 + # Retrieve short name of model class (e.g. "Gauss"
  241 + # or "Disk)
  242 + shorttype = classname.split('Radial')[-1]
  243 +
  244 + # Handle Disk and Shell model
  245 + if (classname == 'GModelSpatialRadialDisk' or
  246 + classname == 'GModelSpatialRadialShell'):
  247 + size = spatial.radius()
  248 +
  249 + # Handle Gauss Model
  250 + elif classname == 'GModelSpatialRadialGauss':
  251 + size = spatial.sigma()
  252 +
  253 + # Skip if source is unknown
  254 + else:
  255 + msg = ('Skip model "%s" since the radial model "%s" '
  256 + 'is unknown.\n' % (model.name(), classname))
  257 + is_valid = False
  258 +
  259 + # Append circle to DS9 string
  260 + if is_valid:
  261 + ds9string += ('circle(%.6f,%.6f,%.6f) #' %
  262 + (modelpos.ra_deg(), modelpos.dec_deg(),
  263 + size*3600.0))
  264 +
  265 + # Handle elliptical sources
  266 + elif 'Elliptical' in classname:
  267 +
  268 + # Retrieve short name and source size
  269 + shorttype = classname.split('Elliptical')[-1]
  270 + size1 = spatial.semimajor()
  271 + size2 = spatial.semiminor()
  272 + angle = spatial.posangle()
  273 +
  274 + # Append ellipse to DS9 string
  275 + ds9string += ('ellipse(%.6f,%.6f,%.6f,%.6f,%.6f) #' %
  276 + (modelpos.ra_deg(), modelpos.dec_deg(),
  277 + size1*3600.0, size2*3600.0, angle+90.0))
  278 +
  279 + # Skip if source is neither radial nor elliptical
  280 + else:
  281 + msg = ('Skip model "%s" since the model "%s" is neither '
  282 + 'a point source, a radial source, nor an elliptical '
  283 + 'source.\n' % (model.name(), classname))
  284 + is_valid = False
  285 +
  286 + # Add short model type to modelname
  287 + if show_ext_type:
  288 + modelname +=' ('+shorttype+')'
  289 +
  290 + # Add DS9 attributes
  291 + if is_valid:
  292 + ds9string += ds9attributes(model, free_color=free_color,
  293 + fixed_color=fixed_color,
  294 + width=width,
  295 + fontfamily=fontfamily,
  296 + fontsize=fontsize,
  297 + fontweight=fontweight,
  298 + fontslant=fontslant)
  299 + if show_labels:
  300 + ds9string += ' text={'+modelname+'}'
  301 + else:
  302 + ds9string = ''
  303 +
  304 + # Return ds9 and message strings
  305 + return ds9string, msg
  306 +
  307 +
  308 +# ========================= #
  309 +# Save models into ds9 file #
  310 +# ========================= #
  311 +def models2ds9file(models, filename,
  312 + pnt_type='cross', pnt_mark_size=12,
  313 + show_labels=True, show_ext_type=True,
  314 + free_color='green', fixed_color='magenta',
  315 + width=2, fontfamily='helvetica', fontsize=12,
  316 + fontweight='normal', fontslant='roman'):
  317 + """
  318 + Save models in a DS9 region file
  319 +
  320 + Parameters
  321 + ----------
  322 + model : `~gammalib.GModel`
  323 + Model
  324 + filename : str
  325 + DS9 region file name
  326 + pnt_type : str, optional
  327 + Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle)
  328 + pnt_mark_size : integer, optional
  329 + Marker size for point sources
  330 + show_labels : boolean, optional
  331 + Add source labels?
  332 + show_ext_type : boolean, optional
  333 + Show type of extended model in source name?
  334 + free_color : str, optional
  335 + Color for sources with free parameters (any DS9 color or hex code)
  336 + fixed_color : str, optional
  337 + Color for source without free parameters (any DS9 color or hex code)
  338 + width : int, optional
  339 + Line width for regions
  340 + fontfamily : str, optional
  341 + Font family for source labels (helvetica,times,courier)
  342 + fontsize : int, optional
  343 + Font size for source labels
  344 + fontweight : str, optional
  345 + Font weight for source labels (normal,bold)
  346 + fontslant : str, optional
  347 + Font slant for source labels (roman,italic)
  348 +
  349 + Returns
  350 + -------
  351 + errors : str
  352 + Error message
  353 + """
  354 + # Initialise error string
  355 + errors = ''
  356 +
  357 + # Open file
  358 + f = open(filename, 'w')
  359 +
  360 + # Write coordinate system
  361 + f.write('fk5\n')
  362 +
  363 + # Loop over models
  364 + for model in models:
  365 +
  366 + # Continue only if point source or extended source model
  367 + if (model.type() == 'PointSource' or
  368 + model.type() == 'ExtendedSource'):
  369 + line, msg = model2ds9string(model,
  370 + pnt_type=pnt_type,
  371 + pnt_mark_size=pnt_mark_size,
  372 + free_color=free_color,
  373 + fixed_color=fixed_color,
  374 + width=width,
  375 + fontfamily=fontfamily,
  376 + fontsize=fontsize,
  377 + fontweight=fontweight,
  378 + fontslant=fontslant)
  379 + if len(line):
  380 + f.write(line+'\n')
  381 + elif len(msg):
  382 + errors += msg+'\n'
  383 +
  384 + # Logging for diffuse components
  385 + elif model.type() == 'DiffuseSource':
  386 + errors += 'Skipping diffuse model "'+model.name()+'".\n'
  387 +
  388 + # Logging for background components
  389 + else:
  390 + errors += 'Skipping background model "'+model.name()+'".\n'
  391 +
  392 + # Close file
  393 + f.close()
  394 +
  395 + # Return error message string
  396 + return errors
... ...
doc/source/users/reference_manual/cssrcdetect.rst 0 → 100644
  View file @6dfd751
  1 +.. _cssrcdetect:
  2 +
  3 +cssrcdetect
  4 +===========
  5 +
  6 +Detects sources in sky map and writes results into a model definition XML file.
  7 +
  8 +
  9 +Synopsis
  10 +--------
  11 +
  12 +This script detects sources in a sky map using a peak detection method and
  13 +writes the detected sources into a model definition XML file. In addition,
  14 +the script produces also a DS9 region file of the detected sources that can
  15 +be overlayed over a sky map displayed with DS9.
  16 +
  17 +For the moment, the script writes out all sources as point sources. Optionally,
  18 +a background model can be added to the model definition XML file so that the
  19 +output file can be used directly as input of a :doc:`ctlike` run.
  20 +
  21 +
  22 +General parameters
  23 +------------------
  24 +
  25 +``inmap [file]``
  26 + Input sky map file.
  27 +
  28 +``outmodel [file]``
  29 + Output model definition XML file.
  30 +
  31 +``outds9file [file]``
  32 + Output DS9 region file.
  33 +
  34 +``srcmodel <POINT> [string]``
  35 + Source model type. For the moment, all detected sources will be added as
  36 + point sources with power law spectral shapes to the model definition
  37 + XML file.
  38 +
  39 +``bkgmodel <NONE|IRF|AEFF|CUBE> [string]``
  40 + Background model type. Using ``NONE`` no background model will be added
  41 + to the model definition XML file. Using ``IRF`` a background model based
  42 + on the template information in the Instrument Response Function will be
  43 + added. Using ``AEFF`` a background model based on the shape of the effective
  44 + area will be added. Using ``CUBE`` a background model for stacked analysis
  45 + will be added.
  46 +
  47 +``threshold [real]``
  48 + Detection threshold (Gaussian sigma). Only sources above this detection
  49 + threshold will be added to the model definiton XML file.
  50 +
  51 +``(maxsrcs = 20) [integer]``
  52 + Maximum number of sources that should be detected.
  53 +
  54 +``(exclrad = 0.2) [real]``
  55 + Radius around a detected source that is excluded from further source
  56 + detection (degrees).
  57 +
  58 +``(fit_pos = yes) [boolean]``
  59 + Enable fitting of source positions in model definiton XML file?
  60 +
  61 +``(fit_shape = yes) [boolean]``
  62 + Enable fitting of source shapes in model definiton XML file?
  63 +
  64 +
  65 +Standard parameters
  66 +-------------------
  67 +
  68 +``(chatter = 2) [integer]``
  69 + Verbosity of the executable:
  70 + chatter = 0: no information will be logged
  71 +
  72 + chatter = 1: only errors will be logged
  73 +
  74 + chatter = 2: errors and actions will be logged
  75 +
  76 + chatter = 3: report about the task execution
  77 +
  78 + chatter = 4: detailed report about the task execution
  79 +
  80 +``(clobber = yes) [boolean]``
  81 + Specifies whether an existing output model file should be overwritten.
  82 +
  83 +``(debug = no) [boolean]``
  84 + Enables debug mode. In debug mode the executable will dump any log file
  85 + output to the console.
  86 +
  87 +``(mode = ql) [string]``
  88 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
  89 +
  90 +``(logfile = cssrcdetect.log) [filename]``
  91 + Log filename.
  92 +
  93 +
  94 +Related tools or scripts
  95 +------------------------
  96 +
  97 +None
... ...
doc/source/users/reference_manual/ctskymap.rst
  View file @6dfd751
... ... @@ -21,56 +21,68 @@ If ``IRF`` is selected, the background template that are shipped with the
21 21 Instrument Response Functions will be used for background subtraction.
22 22  
23 23 ctskymap generates a FITS file comprising a sky map as primary extension.
  24 +If a background subtraction method was selected, the FITS file will contain
  25 +the additional extensions ``BACKGRAOUND`` and ``SIGNIFICANCE`` that contain
  26 +the background map and a significance map, respectively. For the latter, the
  27 +significance of the signal is computed for each sky map pixel.
24 28  
25 29  
26 30 General parameters
27 31 ------------------
28 32  
29 33 ``inobs [file]``
30   - Input event list or observation definition XML file
31   -
  34 + Input event list or observation definition XML file.
  35 +
  36 +``caldb [string]``
  37 + Calibration database (only required for IRF background subtraction if no
  38 + response information is provided by ``inobs``).
  39 +
  40 +``irf [string]``
  41 + Instrument response function (only required for IRF background subtraction
  42 + if no response information is provided by ``inobs``).
  43 +
32 44 ``outmap [file]``
33   - Output sky map file
  45 + Output sky map file.
34 46  
35 47 ``emin [real]``
36   - Minimum energy in map (in TeV)
  48 + Minimum energy in map (in TeV).
37 49  
38 50 ``emax [real]``
39   - Maximum energy in map (in TeV)
  51 + Maximum energy in map (in TeV).
40 52  
41 53 ``(usepnt = no) [boolean]``
42 54 Use pointing instead of xref/yref parameters?
43 55  
44 56 ``nxpix [integer]``
45   - Size of the Right Ascension / Galactic longitude axis (in pixels)
  57 + Size of the Right Ascension / Galactic longitude axis (in pixels).
46 58  
47 59 ``nypix [integer]``
48   - Size of the Declination / Galactic latitude axis (in pixels)
  60 + Size of the Declination / Galactic latitude axis (in pixels).
49 61  
50 62 ``binsz [real]``
51 63 Pixel size (in degrees/pixel).
52 64  
53 65 ``coordsys <CEL|GAL> [string]``
54   - Coordinate system (CEL - celestial, GAL - galactic)
  66 + Coordinate system (CEL - celestial, GAL - galactic).
55 67  
56 68 ``proj <AIT|AZP|CAR|MER|STG|TAN> [string]``
57   - Projection method
  69 + Projection method.
58 70  
59 71 ``xref [real]``
60   - Right Ascension / Galactic longitude of image centre (J2000, in degrees)
  72 + Right Ascension / Galactic longitude of image centre (J2000, in degrees).
61 73  
62 74 ``yref [real]``
63   - Declination / Galactic latitude of image centre (J2000, in degrees)
  75 + Declination / Galactic latitude of image centre (J2000, in degrees).
64 76  
65 77 ``bkgsubtract <NONE|IRF> [string]``
66   - Background subtraction method
  78 + Background subtraction method.
67 79  
68 80  
69 81 Standard parameters
70 82 -------------------
71 83  
72 84 ``(publish = no) [boolean]``
73   - Specifies whether the sky map should be published on VO Hub
  85 + Specifies whether the sky map should be published on VO Hub.
74 86  
75 87 ``(chatter = 2) [integer]``
76 88 Verbosity of the executable:
... ... @@ -85,16 +97,17 @@ Standard parameters
85 97 chatter = 4: detailed report about the task execution
86 98  
87 99 ``(clobber = yes) [boolean]``
88   - Specifies whether an existing output sky map file should be overwritten
  100 + Specifies whether an existing output sky map file should be overwritten.
89 101  
90 102 ``(debug = no) [boolean]``
91   - Enables debug mode. In debug mode the executable will dump any log file output to the console.
  103 + Enables debug mode. In debug mode the executable will dump any log file
  104 + output to the console.
92 105  
93 106 ``(mode = ql) [string]``
94   - Mode of automatic parameters (default is "ql", i.e. "query and learn")
  107 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
95 108  
96 109 ``(logfile = ctskymap.log) [string]``
97   - Name of log file
  110 + Name of log file.
98 111  
99 112  
100 113 Related tools or scripts
... ...
doc/source/users/reference_manual/index.rst
  View file @6dfd751
... ... @@ -40,7 +40,6 @@ cscripts
40 40 :maxdepth: 1
41 41  
42 42 cscaldb --- Lists available instrument response functions <cscaldb>
43   - csviscube --- Computes visibility cube <csviscube>
44 43 cslightcrv --- Computes lightcurve <cslightcrv>
45 44 csmodelinfo --- Shows model container content <csmodelinfo>
46 45 csmodelmerge --- Merges several model containers into one file <csmodelmerge>
... ... @@ -52,11 +51,12 @@ cscripts
52 51 csresmap --- Generates residual map <csresmap>
53 52 cssens --- Computes CTA sensitivity <cssens>
54 53 csspec --- Computes spectral points <csspec>
  54 + cssrcdetect --- Detects sources in sky map <cssrcdetect>
55 55 cstsdist --- Generates Test Statistic distribution <cstsdist>
56 56 cstsmapsplit --- Creates commands to split the Test Statistic map computations <cstsmapsplit>
57 57 cstsmapmerge --- Merges slices from Test Statistic map computations <cstsmapmerge>
58   -
59   -
  58 + csviscube --- Computes visibility cube <csviscube>
  59 +
60 60  
61 61 Scripts to manage an IACT database
62 62 ----------------------------------
... ...
src/ctskymap/ctskymap.cpp
  View file @6dfd751
... ... @@ -233,11 +233,23 @@ void ctskymap::run(void)
233 233 // Map events into sky map
234 234 map_events(obs);
235 235  
236   - // Subtract background from sky map
237   - bkg_subtract(obs);
  236 + // Compute background sky map
  237 + map_background(obs);
238 238  
239 239 } // endfor: looped over observations
240 240  
  241 + // If background subtract was required then compute significance map
  242 + // and subtract background map from counts map
  243 + if (m_bkgsubtract != "NONE") {
  244 +
  245 + // Compute significance map
  246 + m_sigmap = (m_skymap - m_bkgmap) / sqrt(m_skymap);
  247 +
  248 + // Subtract background map from counts map
  249 + m_skymap -= m_bkgmap;
  250 +
  251 + } // endif: background subtraction was requested
  252 +
241 253 // Optionally publish sky map
242 254 if (m_publish) {
243 255 publish();
... ... @@ -271,6 +283,24 @@ void ctskymap::save(void)
271 283 // Write sky map into FITS file
272 284 m_skymap.write(fits);
273 285  
  286 + // If background subtraction is requested then write background map
  287 + // and significance map to FITS file
  288 + if (m_bkgsubtract != "NONE") {
  289 +
  290 + // Write background map into FITS file
  291 + m_bkgmap.write(fits);
  292 +
  293 + // Set background map extension name
  294 + fits[fits.size()-1]->extname("BACKGROUND");
  295 +
  296 + // Write significance map into FITS file
  297 + m_sigmap.write(fits);
  298 +
  299 + // Set significance map extension name
  300 + fits[fits.size()-1]->extname("SIGNIFICANCE");
  301 +
  302 + } // endif: background subtraction was requested
  303 +
274 304 // Save FITS file to disk
275 305 fits.saveto(m_outmap, clobber());
276 306  
... ... @@ -328,6 +358,8 @@ void ctskymap::init_members(void)
328 358 {
329 359 // Initialise members
330 360 m_skymap.clear();
  361 + m_bkgmap.clear();
  362 + m_sigmap.clear();
331 363 m_emin = 0.0;
332 364 m_emax = 0.0;
333 365 m_bkgsubtract = "NONE";
... ... @@ -348,6 +380,8 @@ void ctskymap::copy_members(const ctskymap&amp; app)
348 380 {
349 381 // Copy attributes
350 382 m_skymap = app.m_skymap;
  383 + m_bkgmap = app.m_bkgmap;
  384 + m_sigmap = app.m_sigmap;
351 385 m_emin = app.m_emin;
352 386 m_emax = app.m_emax;
353 387 m_bkgsubtract = app.m_bkgsubtract;
... ... @@ -387,10 +421,18 @@ void ctskymap::get_parameters(void)
387 421 // Create sky map based on task parameters
388 422 m_skymap = create_map(m_obs);
389 423  
390   - // Get remaining parameters
  424 + // Get further parameters
391 425 m_emin = (*this)["emin"].real();
392 426 m_emax = (*this)["emax"].real();
393 427 m_bkgsubtract = (*this)["bkgsubtract"].string();
  428 +
  429 + // If IRF background subtraction is requested then make sure that the
  430 + // CTA observations in the observation container have response information
  431 + if (m_bkgsubtract == "IRF") {
  432 + set_response(m_obs);
  433 + }
  434 +
  435 + // Get remaining parameters
394 436 m_publish = (*this)["publish"].boolean();
395 437 m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
396 438  
... ... @@ -399,6 +441,13 @@ void ctskymap::get_parameters(void)
399 441 m_outmap = (*this)["outmap"].filename();
400 442 }
401 443  
  444 + // Create background map and significance map if background subtraction
  445 + // is requested
  446 + if (m_bkgsubtract != "NONE") {
  447 + m_bkgmap = create_map(m_obs);
  448 + m_sigmap = create_map(m_obs);
  449 + }
  450 +
402 451 // Write parameters into logger
403 452 log_parameters(TERSE);
404 453  
... ... @@ -487,22 +536,23 @@ void ctskymap::map_events(GCTAObservation* obs)
487 536  
488 537  
489 538 /***********************************************************************//**
490   - * @brief Subtract background events from skymap
  539 + * @brief Estimates the background in sky map for observation
491 540 *
492 541 * @param[in] obs CTA observation.
493 542 *
494   - * Subtracts background events from the sky map using various methods. The
495   - * background subtraction method is specified by the "bkgsubtract" parameter
496   - * which can take the following values:
  543 + * Estimates the background in the sky map for a given observation and adds
  544 + * this estimate to the background sky map. The background estimation method
  545 + * is specified by the "bkgsubtract" parameter which can take the following
  546 + * values:
497 547 *
498   - * NONE - No background subtraction
499   - * IRF - Background subtraction based on IRF template
  548 + * NONE - No background estimation
  549 + * IRF - Background estimation based on IRF template
500 550 ***************************************************************************/
501   -void ctskymap::bkg_subtract(GCTAObservation* obs)
  551 +void ctskymap::map_background(GCTAObservation* obs)
502 552 {
503   - // Dispatch to apporpriate background subtraction method
  553 + // Dispatch to appropriate background estimation method
504 554 if (m_bkgsubtract == "IRF") {
505   - bkg_subtract_irf(obs);
  555 + map_background_irf(obs);
506 556 }
507 557  
508 558 // Return
... ... @@ -511,7 +561,7 @@ void ctskymap::bkg_subtract(GCTAObservation* obs)
511 561  
512 562  
513 563 /***********************************************************************//**
514   - * @brief Subtract background events based on IRF template from sky map
  564 + * @brief Estimates the background in sky map based on IRF template
515 565 *
516 566 * @param[in] obs CTA observation.
517 567 *
... ... @@ -519,9 +569,10 @@ void ctskymap::bkg_subtract(GCTAObservation* obs)
519 569 * No response information available for observation.
520 570 * No background template available in instrument response function.
521 571 *
522   - * Subtracts background events from the sky map using the IRF template.
  572 + * Estimates the background in the sky map using the IRF template for a given
  573 + * observation and adds this estimate to the background sky map.
523 574 ***************************************************************************/
524   -void ctskymap::bkg_subtract_irf(GCTAObservation* obs)
  575 +void ctskymap::map_background_irf(GCTAObservation* obs)
525 576 {
526 577 // Get IRF response
527 578 const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
... ... @@ -563,11 +614,11 @@ void ctskymap::bkg_subtract_irf(GCTAObservation* obs)
563 614 // Setup energy nodes for computation
564 615 GEnergies energies(nebin, emin, emax); // More complicated nbins later
565 616  
566   - // Loop over all sky map pixels
567   - for (int i = 0; i < m_skymap.npix(); ++i) {
  617 + // Loop over all background map pixels
  618 + for (int i = 0; i < m_bkgmap.npix(); ++i) {
568 619  
569 620 // Get sky direction of pixel
570   - GSkyDir skydir = m_skymap.inx2dir(i);
  621 + GSkyDir skydir = m_bkgmap.inx2dir(i);
571 622  
572 623 // Convert sky direction in instrument direction
573 624 GCTAInstDir instdir = obs->pointing().instdir(skydir);
... ... @@ -587,12 +638,12 @@ void ctskymap::bkg_subtract_irf(GCTAObservation* obs)
587 638 }
588 639  
589 640 // Multiply background rate with livetime and solid angle
590   - value *= obs->livetime() * m_skymap.solidangle(i);
  641 + value *= obs->livetime() * m_bkgmap.solidangle(i);
591 642  
592   - // Subtract background value from sky map
593   - m_skymap(i,0) -= value;
  643 + // Add number of background events to background map
  644 + m_bkgmap(i,0) += value;
594 645  
595   - } // endfor: looped over sky map pixels
  646 + } // endfor: looped over background map pixels
596 647  
597 648 // Return
598 649 return;
... ...
src/ctskymap/ctskymap.hpp
  View file @6dfd751
... ... @@ -70,8 +70,8 @@ protected:
70 70 void free_members(void);
71 71 void get_parameters(void);
72 72 void map_events(GCTAObservation* obs);
73   - void bkg_subtract(GCTAObservation* obs);
74   - void bkg_subtract_irf(GCTAObservation* obs);
  73 + void map_background(GCTAObservation* obs);
  74 + void map_background_irf(GCTAObservation* obs);
75 75  
76 76 // User parameters
77 77 GFilename m_outmap; //!< Output file name
... ... @@ -83,6 +83,8 @@ protected:
83 83  
84 84 // Protected members
85 85 GSkyMap m_skymap; //!< Sky map
  86 + GSkyMap m_bkgmap; //!< Background map
  87 + GSkyMap m_sigmap; //!< Significance map
86 88 };
87 89  
88 90  
... ...
src/ctskymap/ctskymap.par
  View file @6dfd751
... ... @@ -17,6 +17,8 @@
17 17 # General parameters
18 18 #===================
19 19 inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
  20 +caldb, s, a, "prod2",,, "Calibration database"
  21 +irf, s, a, "South_0.5h",,, "Instrument response function"
20 22 outmap, f, a, "skymap.fits",,, "Output skymap file"
21 23  
22 24 #
... ...
test/data/skymap.fits 0 → 100644
View file @6dfd751
No preview for this file type
test/test_cssrcdetect.py 0 → 100755
  View file @6dfd751
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the cssrcdetect script
  4 +#
  5 +# Copyright (C) 2016 Juergen Knoedlseder
  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 os
  22 +import gammalib
  23 +import cscripts
  24 +from testing import test
  25 +
  26 +
  27 +# ================================= #
  28 +# Test class for cssrcdetect script #
  29 +# ================================= #
  30 +class Test(test):
  31 + """
  32 + Test class for cssrcdetect script
  33 +
  34 + This test class makes unit tests for the cssrcdetect script by using it
  35 + from the command line and from Python.
  36 + """
  37 +
  38 + # Constructor
  39 + def __init__(self):
  40 + """
  41 + Constructor.
  42 + """
  43 + # Call base class constructor
  44 + test.__init__(self)
  45 +
  46 + # Set members
  47 + self._map = self._datadir + '/skymap.fits'
  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('cssrcdetect')
  59 +
  60 + # Append tests
  61 + self.append(self._test_cmd, 'Test cssrcdetect on command line')
  62 + self.append(self._test_python, 'Test cssrcdetect from Python')
  63 +
  64 + # Return
  65 + return
  66 +
  67 + # Test cssrcdetect on command line
  68 + def _test_cmd(self):
  69 + """
  70 + Test cssrcdetect on the command line
  71 + """
  72 + # Set script name
  73 + cssrcdetect = self._script('cssrcdetect')
  74 +
  75 + # Setup cssrcdetect command
  76 + cmd = cssrcdetect+' inmap="'+self._map+'"'+ \
  77 + ' outmodel="cssrcdetect_cmd1.xml"'+ \
  78 + ' outds9file="cssrcdetect_cmd1.reg"'+ \
  79 + ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
  80 + ' logfile="cssrcdetect_cmd1.log" chatter=1'
  81 +
  82 + # Check if execution of wrong command fails
  83 + self.test_assert(self._execute('command_that_does_not_exist') != 0,
  84 + 'Self test of test script')
  85 +
  86 + # Check if execution was successful
  87 + self.test_assert(self._execute(cmd) == 0,
  88 + 'Check successful execution from command line')
  89 +
  90 + # Check model file
  91 + self._check_model_file('cssrcdetect_cmd1.xml', 2)
  92 +
  93 + # Setup cssrcdetect command
  94 + cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \
  95 + ' outmodel="cssrcdetect_cmd1.xml"'+ \
  96 + ' outds9file="cssrcdetect_cmd1.reg"'+ \
  97 + ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
  98 + ' logfile="cssrcdetect_cmd2.log" chatter=2'
  99 +
  100 + # Check if execution failed
  101 + self.test_assert(self._execute(cmd) != 0,
  102 + 'Check invalid input file when executed from command line')
  103 +
  104 + # Return
  105 + return
  106 +
  107 + # Test cssrcdetect from Python
  108 + def _test_python(self):
  109 + """
  110 + Test cssrcdetect from Python
  111 + """
  112 + # Allocate empty cssrcdetect script
  113 + srcdetect = cscripts.cssrcdetect()
  114 +
  115 + # Check that saving does not nothing
  116 + srcdetect['outmodel'] = 'cssrcdetect_py0.xml'
  117 + srcdetect['outds9file'] = 'cssrcdetect_py0.reg'
  118 + srcdetect['logfile'] = 'cssrcdetect_py0.log'
  119 + srcdetect.logFileOpen()
  120 + srcdetect.save()
  121 +
  122 + # Check than an empty model definition file has been created
  123 + self._check_model_file('cssrcdetect_py0.xml', 0)
  124 +
  125 + # Check that clearing does not lead to an exception or segfault
  126 + srcdetect.clear()
  127 +
  128 + # Now set cssrcdetect parameters
  129 + srcdetect['inmap'] = self._map
  130 + srcdetect['outmodel'] = 'cssrcdetect_py1.xml'
  131 + srcdetect['outds9file'] = 'cssrcdetect_py1.reg'
  132 + srcdetect['srcmodel'] = 'POINT'
  133 + srcdetect['bkgmodel'] = 'NONE'
  134 + srcdetect['threshold'] = 10.0
  135 + srcdetect['logfile'] = 'cssrcdetect_py1.log'
  136 + srcdetect['chatter'] = 2
  137 +
  138 + # Run cssrcdetect script and save models
  139 + srcdetect.logFileOpen() # Make sure we get a log file
  140 + srcdetect.run()
  141 + srcdetect.save()
  142 +
  143 + # Check model file
  144 + self._check_model_file('cssrcdetect_py1.xml', 2)
  145 +
  146 + # Execute cssrcdetect script again, now with a higher chatter level
  147 + # than before
  148 + srcdetect['outmodel'] = 'cssrcdetect_py2.xml'
  149 + srcdetect['outds9file'] = 'cssrcdetect_py2.reg'
  150 + srcdetect['bkgmodel'] = 'IRF'
  151 + srcdetect['logfile'] = 'cssrcdetect_py2.log'
  152 + srcdetect['chatter'] = 3
  153 + srcdetect.logFileOpen() # Needed to get a new log file
  154 + srcdetect.execute()
  155 +
  156 + # Check model file
  157 + self._check_model_file('cssrcdetect_py2.xml', 3)
  158 +
  159 + # Return
  160 + return
  161 +
  162 + # Check model file
  163 + def _check_model_file(self, filename, number):
  164 + """
  165 + Check model file
  166 + """
  167 + # Open models
  168 + models = gammalib.GModels(filename)
  169 +
  170 + # Check number of models
  171 + self.test_value(models.size(), number, 'Check for '+str(number)+
  172 + ' models in XML file')
  173 +
  174 + # Return
  175 + return
... ...
test/test_ctskymap.py
  View file @6dfd751
... ... @@ -70,7 +70,7 @@ class Test(test):
70 70 ' outmap="ctskymap_cmd1.fits"'+ \
71 71 ' emin=0.1 emax=100.0 nxpix=200 nypix=200'+ \
72 72 ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \
73   - ' xref=83.63 yref=22.01'+ \
  73 + ' xref=83.63 yref=22.01 bkgsubtract="NONE"'+ \
74 74 ' logfile="ctskymap_cmd1.log" chatter=1'
75 75  
76 76 # Check if execution of wrong command fails
... ... @@ -89,7 +89,7 @@ class Test(test):
89 89 ' outmap="ctskymap_cmd2.fits"'+ \
90 90 ' emin=0.1 emax=100.0 nxpix=200 nypix=200'+ \
91 91 ' binsz=0.02 coordsys="CEL" proj="CAR"'+ \
92   - ' xref=83.63 yref=22.01'+ \
  92 + ' xref=83.63 yref=22.01 bkgsubtract="NONE"'+ \
93 93 ' logfile="ctskymap_cmd2.log" chatter=2'
94 94  
95 95 # Check if execution failed
... ... @@ -134,19 +134,20 @@ class Test(test):
134 134 skymap.clear()
135 135  
136 136 # Now set ctskymap parameters
137   - skymap['inobs'] = self._events
138   - skymap['emin'] = 0.1
139   - skymap['emax'] = 100
140   - skymap['nxpix'] = 200
141   - skymap['nypix'] = 200
142   - skymap['binsz'] = 0.02
143   - skymap['coordsys'] = 'CEL'
144   - skymap['proj'] = 'CAR'
145   - skymap['xref'] = 83.63
146   - skymap['yref'] = 22.01
147   - skymap['outmap'] = 'ctskymap_py1.fits'
148   - skymap['logfile'] = 'ctskymap_py1.log'
149   - skymap['chatter'] = 2
  137 + skymap['inobs'] = self._events
  138 + skymap['emin'] = 0.1
  139 + skymap['emax'] = 100
  140 + skymap['nxpix'] = 200
  141 + skymap['nypix'] = 200
  142 + skymap['binsz'] = 0.02
  143 + skymap['coordsys'] = 'CEL'
  144 + skymap['proj'] = 'CAR'
  145 + skymap['xref'] = 83.63
  146 + skymap['yref'] = 22.01
  147 + skymap['bkgsubtract'] = 'NONE'
  148 + skymap['outmap'] = 'ctskymap_py1.fits'
  149 + skymap['logfile'] = 'ctskymap_py1.log'
  150 + skymap['chatter'] = 2
150 151  
151 152 # Run ctskymap tool
152 153 skymap.logFileOpen() # Make sure we get a log file
... ... @@ -190,18 +191,19 @@ class Test(test):
190 191  
191 192 # Allocate ctskymap tool from observation container
192 193 skymap = ctools.ctskymap(obs)
193   - skymap['emin'] = 0.1
194   - skymap['emax'] = 100
195   - skymap['nxpix'] = 200
196   - skymap['nypix'] = 200
197   - skymap['binsz'] = 0.02
198   - skymap['coordsys'] = 'GAL'
199   - skymap['proj'] = 'CAR'
200   - skymap['xref'] = 184.5575
201   - skymap['yref'] = -5.7844
202   - skymap['outmap'] = 'ctskymap_py3.fits'
203   - skymap['logfile'] = 'ctskymap_py3.log'
204   - skymap['chatter'] = 4
  194 + skymap['emin'] = 0.1
  195 + skymap['emax'] = 100
  196 + skymap['nxpix'] = 200
  197 + skymap['nypix'] = 200
  198 + skymap['binsz'] = 0.02
  199 + skymap['coordsys'] = 'GAL'
  200 + skymap['proj'] = 'CAR'
  201 + skymap['xref'] = 184.5575
  202 + skymap['yref'] = -5.7844
  203 + skymap['bkgsubtract'] = 'NONE'
  204 + skymap['outmap'] = 'ctskymap_py3.fits'
  205 + skymap['logfile'] = 'ctskymap_py3.log'
  206 + skymap['chatter'] = 4
205 207  
206 208 # Execute tool
207 209 skymap.logFileOpen() # Needed to get a new log file
... ... @@ -213,6 +215,31 @@ class Test(test):
213 215 # Publish with name
214 216 skymap.publish('My sky map')
215 217  
  218 + # Allocate ctskymap tool from observation container
  219 + skymap = ctools.ctskymap(obs)
  220 + skymap['emin'] = 0.1
  221 + skymap['emax'] = 100
  222 + skymap['nxpix'] = 200
  223 + skymap['nypix'] = 200
  224 + skymap['binsz'] = 0.02
  225 + skymap['coordsys'] = 'GAL'
  226 + skymap['proj'] = 'CAR'
  227 + skymap['xref'] = 184.5575
  228 + skymap['yref'] = -5.7844
  229 + skymap['bkgsubtract'] = 'IRF'
  230 + skymap['caldb'] = self._caldb
  231 + skymap['irf'] = self._irf
  232 + skymap['outmap'] = 'ctskymap_py4.fits'
  233 + skymap['logfile'] = 'ctskymap_py4.log'
  234 + skymap['chatter'] = 4
  235 +
  236 + # Execute tool
  237 + skymap.logFileOpen() # Needed to get a new log file
  238 + skymap.execute()
  239 +
  240 + # Check result file
  241 + self._check_result_file('ctskymap_py3.fits')
  242 +
216 243 # Return
217 244 return
218 245  
... ...
test/test_python_cscripts.py
  View file @6dfd751
... ... @@ -40,6 +40,7 @@ import test_csresmap
40 40 import test_csroot2caldb
41 41 import test_cssens
42 42 import test_csspec
  43 +import test_cssrcdetect
43 44 import test_cstsdist
44 45 import test_cstsmapmerge
45 46 import test_cstsmapsplit
... ... @@ -126,6 +127,7 @@ def test(installed=False, debug=False):
126 127 test_csroot2caldb.Test(),
127 128 test_cssens.Test(),
128 129 test_csspec.Test(),
  130 + test_cssrcdetect.Test(),
129 131 test_cstsdist.Test(),
130 132 test_cstsmapmerge.Test(),
131 133 test_cstsmapsplit.Test(),
... ...