jean perrot / ctools

Commit 6049bf07b4d985a8e80d55dcfcdfab710e1c4ba5

Authored by Jürgen Knödlseder
1 parent dd903c61

Add cssrcdetect script (#1617)

Inline Side-by-side
Showing 13 changed files with 1127 additions and 191 deletions
ChangeLog
  View file @6049bf0
1 -2016-11-25 1 +2016-12-01
2 2
3 * Version 1.2.0 released 3 * Version 1.2.0 released
4 ======================== 4 ========================
5 5
  6 + Add cssrcdetect script (#1617)
6 Add background subtraction to ctskymap tool (#1882) 7 Add background subtraction to ctskymap tool (#1882)
7 Add csobsselect script (#1881) 8 Add csobsselect script (#1881)
8 Add csmodelselect script (#1879) 9 Add csmodelselect script (#1879)
1 New Features and Important Changes in ctools 1.2.0 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 Introduction 6 Introduction
@@ -206,30 +206,33 @@ The class now derives from ctlikelihood (#1828). @@ -206,30 +206,33 @@ The class now derives from ctlikelihood (#1828).
206 206
207 cscripts 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 "obsutils" module (#1830). 211 "obsutils" module (#1830).
212 212
213 The csroot2caldb script now takes the zenith and azimuth angles as parameters 213 The csroot2caldb script now takes the zenith and azimuth angles as parameters
214 and writes them into the calibration database (#1835). 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 of hours that a given position on the sky is observable under a given zenith 217 of hours that a given position on the sky is observable under a given zenith
218 angle by an IACT array during a given time period. The tool implements so far 218 angle by an IACT array during a given time period. The tool implements so far
219 a Sun constraint, but no Moon constraint (#1835). 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 Fix energy boundary handling in cssens that has been corrupted during code 224 Fix energy boundary handling in cssens that has been corrupted during code
225 refactoring (#1873). 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 XML file (#1879). 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 definition XML file (#1881). 231 definition XML file (#1881).
232 232
  233 +Add cssrcdetect script to detect sources in a sky map (#1617).
  234 +
  235 +
233 236
234 Examples 237 Examples
235 -------- 238 --------
cscripts/Makefile.am
  View file @6049bf0
@@ -36,6 +36,7 @@ pythonscripts = $(builddir)/__init__.py \ @@ -36,6 +36,7 @@ pythonscripts = $(builddir)/__init__.py \
36 $(srcdir)/csroot2caldb.py \ 36 $(srcdir)/csroot2caldb.py \
37 $(srcdir)/cssens.py \ 37 $(srcdir)/cssens.py \
38 $(srcdir)/csspec.py \ 38 $(srcdir)/csspec.py \
  39 + $(srcdir)/cssrcdetect.py \
39 $(srcdir)/cstsdist.py \ 40 $(srcdir)/cstsdist.py \
40 $(srcdir)/cstsmapmerge.py \ 41 $(srcdir)/cstsmapmerge.py \
41 $(srcdir)/cstsmapsplit.py \ 42 $(srcdir)/cstsmapsplit.py \
@@ -65,6 +66,7 @@ parfiles = $(srcdir)/cscaldb.par \ @@ -65,6 +66,7 @@ parfiles = $(srcdir)/cscaldb.par \
65 $(srcdir)/csroot2caldb.par \ 66 $(srcdir)/csroot2caldb.par \
66 $(srcdir)/cssens.par \ 67 $(srcdir)/cssens.par \
67 $(srcdir)/csspec.par \ 68 $(srcdir)/csspec.par \
  69 + $(srcdir)/cssrcdetect.par \
68 $(srcdir)/cstsdist.par \ 70 $(srcdir)/cstsdist.par \
69 $(srcdir)/cstsmapmerge.par \ 71 $(srcdir)/cstsmapmerge.par \
70 $(srcdir)/cstsmapsplit.par \ 72 $(srcdir)/cstsmapsplit.par \
@@ -91,6 +93,7 @@ linkscripts = cscaldb \ @@ -91,6 +93,7 @@ linkscripts = cscaldb \
91 csroot2caldb \ 93 csroot2caldb \
92 cssens \ 94 cssens \
93 csspec \ 95 csspec \
  96 + cssrcdetect \
94 cstsdist \ 97 cstsdist \
95 cstsmapmerge \ 98 cstsmapmerge \
96 cstsmapsplit \ 99 cstsmapsplit \
@@ -119,6 +122,7 @@ testscripts = $(srcdir)/tests/__init__.py \ @@ -119,6 +122,7 @@ testscripts = $(srcdir)/tests/__init__.py \
119 $(top_srcdir)/test/test_csroot2caldb.py \ 122 $(top_srcdir)/test/test_csroot2caldb.py \
120 $(top_srcdir)/test/test_cssens.py \ 123 $(top_srcdir)/test/test_cssens.py \
121 $(top_srcdir)/test/test_csspec.py \ 124 $(top_srcdir)/test/test_csspec.py \
  125 + $(top_srcdir)/test/test_cssrcdetect.py \
122 $(top_srcdir)/test/test_cstsdist.py \ 126 $(top_srcdir)/test/test_cstsdist.py \
123 $(top_srcdir)/test/test_cstsmapmerge.py \ 127 $(top_srcdir)/test/test_cstsmapmerge.py \
124 $(top_srcdir)/test/test_cstsmapsplit.py \ 128 $(top_srcdir)/test/test_cstsmapsplit.py \
@@ -160,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ @@ -160,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
160 $(top_srcdir)/test/data/pull.dat \ 164 $(top_srcdir)/test/data/pull.dat \
161 $(top_srcdir)/test/data/sensitivity.dat \ 165 $(top_srcdir)/test/data/sensitivity.dat \
162 $(top_srcdir)/test/data/sensitivity_integral.dat \ 166 $(top_srcdir)/test/data/sensitivity_integral.dat \
  167 + $(top_srcdir)/test/data/skymap.fits \
163 $(top_srcdir)/test/data/spectrum.fits \ 168 $(top_srcdir)/test/data/spectrum.fits \
164 $(top_srcdir)/test/data/ts.dat \ 169 $(top_srcdir)/test/data/ts.dat \
165 $(top_srcdir)/test/data/models_spatial.xml \ 170 $(top_srcdir)/test/data/models_spatial.xml \
cscripts/__init__.py.in
  View file @6049bf0
@@ -39,6 +39,7 @@ __all__ = [ @@ -39,6 +39,7 @@ __all__ = [
39 "csroot2caldb", 39 "csroot2caldb",
40 "cssens", 40 "cssens",
41 "csspec", 41 "csspec",
  42 + "cssrcdetect",
42 "cstsdist", 43 "cstsdist",
43 "cstsmapmerge", 44 "cstsmapmerge",
44 "cstsmapsplit", 45 "cstsmapsplit",
@@ -71,6 +72,7 @@ from cscripts.csresmap import csresmap @@ -71,6 +72,7 @@ from cscripts.csresmap import csresmap
71 from cscripts.csroot2caldb import csroot2caldb 72 from cscripts.csroot2caldb import csroot2caldb
72 from cscripts.cssens import cssens 73 from cscripts.cssens import cssens
73 from cscripts.csspec import csspec 74 from cscripts.csspec import csspec
  75 +from cscripts.cssrcdetect import cssrcdetect
74 from cscripts.cstsdist import cstsdist 76 from cscripts.cstsdist import cstsdist
75 from cscripts.cstsmapmerge import cstsmapmerge 77 from cscripts.cstsmapmerge import cstsmapmerge
76 from cscripts.cstsmapsplit import cstsmapsplit 78 from cscripts.cstsmapsplit import cstsmapsplit
cscripts/csmodelinfo.py
  View file @6049bf0
@@ -21,6 +21,7 @@ @@ -21,6 +21,7 @@
21 import sys 21 import sys
22 import gammalib 22 import gammalib
23 import ctools 23 import ctools
  24 +from cscripts import modutils
24 25
25 26
26 # ================= # 27 # ================= #
@@ -86,155 +87,6 @@ class csmodelinfo(ctools.cscript): @@ -86,155 +87,6 @@ class csmodelinfo(ctools.cscript):
86 # Return 87 # Return
87 return 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 # Public methods 91 # Public methods
240 def run(self): 92 def run(self):
@@ -407,36 +259,34 @@ class csmodelinfo(ctools.cscript): @@ -407,36 +259,34 @@ class csmodelinfo(ctools.cscript):
407 259
408 # Log filename 260 # Log filename
409 self._log_value(gammalib.NORMAL, 'DS9 filename', ds9file.url()) 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 # Return 290 # Return
441 return 291 return
442 292
cscripts/cssrcdetect.par 0 → 100644
  View file @6049bf0
  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 @6049bf0
  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 @6049bf0
@@ -99,3 +99,298 @@ def test_source(models, srcname, ra=None, dec=None, fitspat=False, @@ -99,3 +99,298 @@ def test_source(models, srcname, ra=None, dec=None, fitspat=False,
99 # Return model container 99 # Return model container
100 return outmodels 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 @6049bf0
  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/index.rst
  View file @6049bf0
@@ -40,7 +40,6 @@ cscripts @@ -40,7 +40,6 @@ cscripts
40 :maxdepth: 1 40 :maxdepth: 1
41 41
42 cscaldb --- Lists available instrument response functions <cscaldb> 42 cscaldb --- Lists available instrument response functions <cscaldb>
43 - csviscube --- Computes visibility cube <csviscube>  
44 cslightcrv --- Computes lightcurve <cslightcrv> 43 cslightcrv --- Computes lightcurve <cslightcrv>
45 csmodelinfo --- Shows model container content <csmodelinfo> 44 csmodelinfo --- Shows model container content <csmodelinfo>
46 csmodelmerge --- Merges several model containers into one file <csmodelmerge> 45 csmodelmerge --- Merges several model containers into one file <csmodelmerge>
@@ -52,11 +51,12 @@ cscripts @@ -52,11 +51,12 @@ cscripts
52 csresmap --- Generates residual map <csresmap> 51 csresmap --- Generates residual map <csresmap>
53 cssens --- Computes CTA sensitivity <cssens> 52 cssens --- Computes CTA sensitivity <cssens>
54 csspec --- Computes spectral points <csspec> 53 csspec --- Computes spectral points <csspec>
  54 + cssrcdetect --- Detects sources in sky map <cssrcdetect>
55 cstsdist --- Generates Test Statistic distribution <cstsdist> 55 cstsdist --- Generates Test Statistic distribution <cstsdist>
56 cstsmapsplit --- Creates commands to split the Test Statistic map computations <cstsmapsplit> 56 cstsmapsplit --- Creates commands to split the Test Statistic map computations <cstsmapsplit>
57 cstsmapmerge --- Merges slices from Test Statistic map computations <cstsmapmerge> 57 cstsmapmerge --- Merges slices from Test Statistic map computations <cstsmapmerge>
58 -  
59 - 58 + csviscube --- Computes visibility cube <csviscube>
  59 +
60 60
61 Scripts to manage an IACT database 61 Scripts to manage an IACT database
62 ---------------------------------- 62 ----------------------------------
test/data/skymap.fits 0 → 100644
View file @6049bf0
No preview for this file type
test/test_cssrcdetect.py 0 → 100755
  View file @6049bf0
  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_python_cscripts.py
  View file @6049bf0
@@ -40,6 +40,7 @@ import test_csresmap @@ -40,6 +40,7 @@ import test_csresmap
40 import test_csroot2caldb 40 import test_csroot2caldb
41 import test_cssens 41 import test_cssens
42 import test_csspec 42 import test_csspec
  43 +import test_cssrcdetect
43 import test_cstsdist 44 import test_cstsdist
44 import test_cstsmapmerge 45 import test_cstsmapmerge
45 import test_cstsmapsplit 46 import test_cstsmapsplit
@@ -126,6 +127,7 @@ def test(installed=False, debug=False): @@ -126,6 +127,7 @@ def test(installed=False, debug=False):
126 test_csroot2caldb.Test(), 127 test_csroot2caldb.Test(),
127 test_cssens.Test(), 128 test_cssens.Test(),
128 test_csspec.Test(), 129 test_csspec.Test(),
  130 + test_cssrcdetect.Test(),
129 test_cstsdist.Test(), 131 test_cstsdist.Test(),
130 test_cstsmapmerge.Test(), 132 test_cstsmapmerge.Test(),
131 test_cstsmapsplit.Test(), 133 test_cstsmapsplit.Test(),