Luigi Tibaldo / ctools

Commit 2bab4887a400df1cd994d08888cc3bac34a40a1b

Authored by Jürgen Knödlseder
2 parents e25a31e4 ab72243e

Merge branch 'integration' into devel

Inline Side-by-side
Showing 4 changed files with 257 additions and 39 deletions
.gitignore
  View file @2bab488
... ... @@ -82,6 +82,7 @@ test/*.reg
82 82 test/*.trs
83 83 test/*.txt
84 84 test/*.eps
  85 +test/*.jpg
85 86 test/test_cscripts.sh
86 87 test/test_examples.py
87 88 test/reports/*.xml
... ...
cscripts/cssens.py
  View file @2bab488
... ... @@ -53,6 +53,7 @@ class cssens(ctools.cscript):
53 53 # Initialise class members
54 54 self._obs = gammalib.GObservations()
55 55 self._ebounds = gammalib.GEbounds()
  56 + self._obs_ebounds = []
56 57 self._srcname = ""
57 58 self._outfile = gammalib.GFilename()
58 59 self._ra = None
... ... @@ -270,8 +271,41 @@ class cssens(ctools.cscript):
270 271 # Loop over all observations in container
271 272 for obs in self._obs:
272 273  
  274 + # Get observation energy boundaries
  275 + obs_ebounds = obs.events().ebounds()
  276 +
  277 + # Get minimum and maximum energy
  278 + obs_emin = obs_ebounds.emin()
  279 + obs_emax = obs_ebounds.emax()
  280 +
  281 + # Case A: bin fully contained in observation ebounds
  282 + if obs_emin <= emin and obs_emax >= emax:
  283 + ebounds = gammalib.GEbounds(emin, emax)
  284 +
  285 + # Case B: bin fully outside of obs ebounds
  286 + elif emax < obs_emin or emin > obs_emax:
  287 +
  288 + # Set zero range (inspired by ctselect)
  289 + e0 = gammalib.GEnergy(0.0, "TeV")
  290 + ebounds = gammalib.GEbounds(e0, e0)
  291 +
  292 + # Case C: bin partly overlapping with observation ebounds
  293 + else:
  294 +
  295 + # Set energy range as obs ebounds were fully contained inside energy bin
  296 + set_emin = emin
  297 + set_emax = emax
  298 +
  299 + # Adjust energy bin to respect observation energy boundary
  300 + if emin < obs_emin:
  301 + set_emin = obs_emin
  302 + if emax > obs_emax:
  303 + set_emax = obs_emax
  304 +
  305 + #Set energy boundaries
  306 + ebounds = gammalib.GEbounds(set_emin, set_emax)
  307 +
273 308 # Set energy boundaries
274   - ebounds = gammalib.GEbounds(emin, emax)
275 309 obs.events().ebounds(ebounds)
276 310  
277 311 # Return
... ... @@ -624,6 +658,10 @@ class cssens(ctools.cscript):
624 658 for par in model:
625 659 self._log(str(par)+"\n")
626 660  
  661 + # Restore energy boundaries of observation container
  662 + for i, obs in enumerate(self._obs):
  663 + obs.events().ebounds(self._obs_ebounds[i])
  664 +
627 665 # Store result
628 666 result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
629 667 'crab_flux': crab_flux, 'photon_flux': photon_flux, \
... ... @@ -645,7 +683,11 @@ class cssens(ctools.cscript):
645 683  
646 684 # Get parameters
647 685 self._get_parameters()
648   -
  686 +
  687 + # Loop over observations and store ebounds
  688 + for obs in self._obs:
  689 + self._obs_ebounds.append(obs.events().ebounds())
  690 +
649 691 # Initialise script
650 692 colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux',
651 693 'energy_flux', 'sensitivity']
... ...
examples/show_obs.py 0 → 100755
  View file @2bab488
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Display summary of observation definition XML file
  4 +#
  5 +# Copyright (C) 2016 Michael Mayer
  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 gammalib
  23 +import ctools
  24 +import cscripts
  25 +try:
  26 + import matplotlib.pyplot as plt
  27 + plt.figure()
  28 + plt.close()
  29 +except:
  30 + print('This script needs the "matplotlib" module')
  31 + sys.exit()
  32 +
  33 +
  34 +# ==================== #
  35 +# Run csobsinfo script #
  36 +# ==================== #
  37 +def run_csobsinfo(filename, ra=None, dec=None, debug=True):
  38 + """
  39 + Run csobsinfo script
  40 +
  41 + Parameters
  42 + ----------
  43 + filename : str
  44 + File name of observation definition XML file
  45 + ra : float, optional
  46 + Target Right Ascension (deg)
  47 + dec : float, optional
  48 + Target declination of pointing (deg)
  49 + debug : bool, optional
  50 + Switch on debugging in csobsinfo run
  51 +
  52 + Returns
  53 + -------
  54 + info : `~cscripts.csobsinfo`
  55 + csobsinfo instance
  56 + """
  57 + # Get observation definition XML filename
  58 + obsdef = gammalib.GFilename(filename)
  59 +
  60 + # Setup csobsinfo script
  61 + info = cscripts.csobsinfo()
  62 + info['inobs'] = obsdef.url()
  63 + info['debug'] = debug
  64 +
  65 + # Set offset
  66 + if ra != None and dec != None:
  67 + info['ra'] = ra
  68 + info['dec'] = dec
  69 + info['offset'] = True
  70 + else:
  71 + info['offset'] = False
  72 +
  73 + # Run csobsinfo
  74 + info.run()
  75 +
  76 + # Return
  77 + return info
  78 +
  79 +
  80 +# ================ #
  81 +# Plot information #
  82 +# ================ #
  83 +def plot_information(info, plotfile=''):
  84 + """
  85 + Plot information
  86 +
  87 + Parameters
  88 + ----------
  89 + info : `~cscripts.csobsinfo`
  90 + csobsinfo instance
  91 + plotfile : str, optional
  92 + Plot filename
  93 + """
  94 + # Retrieve observation info
  95 + zeniths = info.zeniths()
  96 + azimuths = info.azimuths()
  97 + offsets = info.offsets()
  98 + ebounds = info.ebounds()
  99 + gti = info.gti()
  100 +
  101 + # Create figure with subplots
  102 + nrows = 2
  103 + ncols = 2
  104 + if info['offset'].boolean() and ebounds.size() > 0:
  105 + ncols += 1
  106 + iplot = 1
  107 + plt.figure()
  108 + plt.subplot(nrows, ncols, iplot)
  109 +
  110 + # Plot zenith angle distribution
  111 + zmin = min(zeniths)
  112 + zmax = max(zeniths)
  113 + plt.hist(zeniths, bins=30, range=(zmin, zmax), fc='blue')
  114 + plt.xlabel('Zenith Angle (deg)')
  115 + plt.ylabel('Frequency')
  116 + plt.title('Zenith angle distribution')
  117 +
  118 + # Plot azimuth angle distribution
  119 + iplot += 1
  120 + plt.subplot(nrows, ncols, iplot)
  121 + amin = min(azimuths)
  122 + amax = max(azimuths)
  123 + plt.hist(azimuths, bins=30, range=(amin, amax), fc='blue')
  124 + plt.xlabel('Azimuth Angle (deg)')
  125 + plt.ylabel('Frequency')
  126 + plt.title('Azimuth distribution')
  127 +
  128 + # Plot offset if possible
  129 + if info['offset'].boolean():
  130 + iplot += 1
  131 + plt.subplot(nrows, ncols, iplot)
  132 + omin = min(offsets)
  133 + omax = max(offsets)
  134 + plt.hist(offsets, bins=30, range=(omin, omax), fc='blue')
  135 + plt.xlabel('Offset from (RA, DEC)=('+str(ra)+','+str(dec)+') (deg)')
  136 + plt.ylabel('Frequency')
  137 + plt.title('Offset distribution')
  138 +
  139 + # Plot energy thresholds if possible
  140 + if ebounds.size() > 0:
  141 + emin = []
  142 + emax = []
  143 + for i in range(ebounds.size()):
  144 + emin.append(ebounds.emin(i).log10TeV())
  145 + emax.append(ebounds.emax(i).log10TeV())
  146 + iplot += 1
  147 + plt.subplot(nrows, ncols, iplot)
  148 + plt.hist(emin, bins=80, range=(-1.0, 2.0), fc='red', label='emin')
  149 + plt.hist(emax, bins=80, range=(-1.0, 2.0), fc='blue', label='emax')
  150 + plt.xlabel('Energy threshold (log10 (E/TeV))')
  151 + plt.ylabel('Frequency')
  152 + plt.legend(loc='upper left')
  153 + plt.title('Energy threshold')
  154 +
  155 + # Plot observation point in time wrt zenith angle
  156 + iplot += 1
  157 + plt.subplot(nrows, ncols, iplot)
  158 + times = []
  159 + for i in range(gti.size()):
  160 + tmean = gti.tstart(i) + 0.5*(gti.tstart(i)-gti.tstart(i))
  161 + times.append(tmean.mjd())
  162 + plt.plot(times, zeniths, 'o', lw=2.0, color='black')
  163 + plt.xlabel('Time (MJD)')
  164 + plt.ylabel('Zenith Angle (deg)')
  165 + plt.title('Observation time')
  166 +
  167 + # Show plots or save it into file
  168 + if len(plotfile) > 0:
  169 + plt.savefig(plotfile)
  170 + else:
  171 + plt.show()
  172 +
  173 + # Return
  174 + return
  175 +
  176 +
  177 +# ======================== #
  178 +# Main routine entry point #
  179 +# ======================== #
  180 +if __name__ == '__main__':
  181 +
  182 + # Print usage message if wrong number of arguments is provided
  183 + if len(sys.argv) < 2:
  184 + sys.exit('Usage: inspect_obs.py inobs.xml [ra] [dec] [file]')
  185 +
  186 + # Extract optional parameters
  187 + plotfile = ''
  188 + ra = None
  189 + dec = None
  190 + if len(sys.argv) == 3:
  191 + plotfile = sys.argv[2]
  192 + elif len(sys.argv) == 4:
  193 + ra = float(sys.argv[2])
  194 + dec = float(sys.argv[3])
  195 +
  196 + # Run csobsinfo
  197 + info = run_csobsinfo(sys.argv[1], ra=ra, dec=dec)
  198 +
  199 + # Plot information
  200 + plot_information(info, plotfile=plotfile)
... ...
test/test_examples.py.in
  View file @2bab488
... ... @@ -108,6 +108,7 @@ class Test(gammalib.GPythonTestSuite):
108 108 self.append(self.test_show_sensitivity, 'Test test_show_sensitivity.py')
109 109 self.append(self.test_show_spectrum, 'Test test_show_spectrum.py')
110 110 self.append(self.test_show_ts_distribution, 'Test show_ts_distribution.py')
  111 + self.append(self.test_show_obs, 'Test show_obs.py')
111 112  
112 113 # Return
113 114 return
... ... @@ -128,8 +129,6 @@ class Test(gammalib.GPythonTestSuite):
128 129 self._execute_python(script, args='gc', log='make_pointings_gc')
129 130 self._execute_python(script, args='lmc', log='make_pointings_lmc')
130 131  
131   - #TODO: Do any testing
132   -
133 132 # Return
134 133 return
135 134  
... ... @@ -141,8 +140,6 @@ class Test(gammalib.GPythonTestSuite):
141 140 # Execute script
142 141 self._execute_python('generate_prod3_irfs')
143 142  
144   - #TODO: Do any testing
145   -
146 143 # Return
147 144 return
148 145  
... ... @@ -154,8 +151,6 @@ class Test(gammalib.GPythonTestSuite):
154 151 # Execute script
155 152 self._execute_python('make_spectrum', args='eps')
156 153  
157   - #TODO: Do any testing
158   -
159 154 # Return
160 155 return
161 156  
... ... @@ -167,8 +162,6 @@ class Test(gammalib.GPythonTestSuite):
167 162 # Execute script
168 163 self._execute_python('make_ts_distributions', args='-n 2 -e 0 -d 1800')
169 164  
170   - #TODO: Do any testing
171   -
172 165 # Return
173 166 return
174 167  
... ... @@ -180,8 +173,6 @@ class Test(gammalib.GPythonTestSuite):
180 173 # Execute script
181 174 self._execute_python('pipeline_binned_disk')
182 175  
183   - #TODO: Do any testing
184   -
185 176 # Return
186 177 return
187 178  
... ... @@ -193,8 +184,6 @@ class Test(gammalib.GPythonTestSuite):
193 184 # Execute script
194 185 self._execute_python('pipeline_binned_mem')
195 186  
196   - #TODO: Do any testing
197   -
198 187 # Return
199 188 return
200 189  
... ... @@ -206,8 +195,6 @@ class Test(gammalib.GPythonTestSuite):
206 195 # Execute script
207 196 self._execute_python('pipeline_stacked_disk')
208 197  
209   - #TODO: Do any testing
210   -
211 198 # Return
212 199 return
213 200  
... ... @@ -219,8 +206,6 @@ class Test(gammalib.GPythonTestSuite):
219 206 # Execute script
220 207 self._execute_python('pipeline_stacked_mem')
221 208  
222   - #TODO: Do any testing
223   -
224 209 # Return
225 210 return
226 211  
... ... @@ -232,8 +217,6 @@ class Test(gammalib.GPythonTestSuite):
232 217 # Execute script
233 218 self._execute_python('pipeline_unbinned_disk')
234 219  
235   - #TODO: Do any testing
236   -
237 220 # Return
238 221 return
239 222  
... ... @@ -245,8 +228,6 @@ class Test(gammalib.GPythonTestSuite):
245 228 # Execute script
246 229 self._execute_python('pipeline_unbinned_mem')
247 230  
248   - #TODO: Do any testing
249   -
250 231 # Return
251 232 return
252 233  
... ... @@ -260,8 +241,6 @@ class Test(gammalib.GPythonTestSuite):
260 241 args='data/butterfly.txt '
261 242 'example_show_butterfly.eps')
262 243  
263   - #TODO: Do any testing
264   -
265 244 # Return
266 245 return
267 246  
... ... @@ -274,8 +253,6 @@ class Test(gammalib.GPythonTestSuite):
274 253 self._execute_python('show_model',
275 254 args='data/crab.xml Crab example_show_model.eps')
276 255  
277   - #TODO: Do any testing
278   -
279 256 # Return
280 257 return
281 258  
... ... @@ -288,8 +265,6 @@ class Test(gammalib.GPythonTestSuite):
288 265 self._execute_python('show_pha',
289 266 args='data/pha.fits example_show_pha.eps')
290 267  
291   - #TODO: Do any testing
292   -
293 268 # Return
294 269 return
295 270  
... ... @@ -303,8 +278,6 @@ class Test(gammalib.GPythonTestSuite):
303 278 args='data/pull.dat Pull_Crab_Prefactor '
304 279 'example_show_pull_evolution.eps')
305 280  
306   - #TODO: Do any testing
307   -
308 281 # Return
309 282 return
310 283  
... ... @@ -318,8 +291,6 @@ class Test(gammalib.GPythonTestSuite):
318 291 args='data/pull.dat Pull_Crab_Prefactor 50 '
319 292 'example_show_pull_histogram.eps')
320 293  
321   - #TODO: Do any testing
322   -
323 294 # Return
324 295 return
325 296  
... ... @@ -331,8 +302,6 @@ class Test(gammalib.GPythonTestSuite):
331 302 # Execute script
332 303 self._execute_python('show_response', args='example_show_response.eps')
333 304  
334   - #TODO: Do any testing
335   -
336 305 # Return
337 306 return
338 307  
... ... @@ -346,8 +315,6 @@ class Test(gammalib.GPythonTestSuite):
346 315 args='data/sensitivity.dat '
347 316 'example_show_sensitivity.eps')
348 317  
349   - #TODO: Do any testing
350   -
351 318 # Return
352 319 return
353 320  
... ... @@ -361,8 +328,6 @@ class Test(gammalib.GPythonTestSuite):
361 328 args='data/spectrum.fits '
362 329 'example_show_spectrum.eps')
363 330  
364   - #TODO: Do any testing
365   -
366 331 # Return
367 332 return
368 333  
... ... @@ -376,7 +341,17 @@ class Test(gammalib.GPythonTestSuite):
376 341 args='data/ts.dat -n 100 '
377 342 '-f example_show_ts_distribution.eps')
378 343  
379   - #TODO: Do any testing
  344 + # Return
  345 + return
  346 +
  347 + # Test show_obs
  348 + def test_show_obs(self):
  349 + """
  350 + Test show_obs
  351 + """
  352 + # Execute script
  353 + self._execute_python('show_obs',
  354 + args='obs_selected.xml example_show_obs.jpg')
380 355  
381 356 # Return
382 357 return
... ...