Eirik Hatlen / ctools

Commit 3f05501519d97f8151e8942485ead3ca39959e4f

Authored by Jürgen Knödlseder
1 parent ef7bd2a2

Add obsutils unit tests

Inline Side-by-side
Showing 2 changed files with 310 additions and 22 deletions
cscripts/obsutils.py
  View file @3f05501
... ... @@ -849,9 +849,7 @@ def get_onoff_obs(cls, obs, nthreads=0):
849 849  
850 850 # On/Off observations are created with CSTAT as default statistic
851 851 # We will change this to the user choice
852   - if cls['statistic'].string() == 'DEFAULT':
853   - pass
854   - else:
  852 + if cls['statistic'].string() != 'DEFAULT':
855 853 for observation in onoff_obs:
856 854 observation.statistic(cls['statistic'].string())
857 855  
... ...
test/test_obsutils.py
  View file @3f05501
... ... @@ -19,13 +19,14 @@
19 19 #
20 20 # ==========================================================================
21 21 import gammalib
  22 +import cscripts
22 23 from cscripts import obsutils
23 24 from testing import test
24 25  
25 26  
26   -# =============================== #
27   -# Test class for csviscube script #
28   -# =============================== #
  27 +# ============================== #
  28 +# Test class for obsutils module #
  29 +# ============================== #
29 30 class Test(test):
30 31 """
31 32 Test class for csviscube script
... ... @@ -79,6 +80,8 @@ class Test(test):
79 80 'Test obsutils.get_stacked_response() function')
80 81 self.append(self._test_get_stacked_obs,
81 82 'Test obsutils.get_stacked_obs() function')
  83 + self.append(self._test_get_onoff_obs,
  84 + 'Test obsutils.get_onoff_obs() function')
82 85  
83 86 # Return
84 87 return
... ... @@ -122,7 +125,7 @@ class Test(test):
122 125 self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
123 126 'Check minimum energy')
124 127 self.test_value(res[0].events().ebounds().emax().TeV(), 10.0,
125   - 'Check minimum energy')
  128 + 'Check maximum energy')
126 129 self.test_value(res[0].events().number(), 4,
127 130 'Check number of events in list')
128 131  
... ... @@ -153,7 +156,7 @@ class Test(test):
153 156 self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
154 157 'Check minimum energy')
155 158 self.test_value(res[0].events().ebounds().emax().TeV(), 10.0,
156   - 'Check minimum energy')
  159 + 'Check maximum energy')
157 160 self.test_value(res[0].events().ebounds().size(), 5,
158 161 'Check number of energy bins')
159 162 self.test_value(res[0].events().number(), 4,
... ... @@ -186,7 +189,7 @@ class Test(test):
186 189 self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
187 190 'Check minimum energy')
188 191 self.test_value(res[0].events().ebounds().emax().TeV(), 10.0,
189   - 'Check minimum energy')
  192 + 'Check maximum energy')
190 193 self.test_value(res[0].events().ebounds().size(), 5,
191 194 'Check number of energy bins')
192 195 self.test_value(res[0].events().number(), 15,
... ... @@ -220,7 +223,7 @@ class Test(test):
220 223 self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
221 224 'Check minimum energy')
222 225 self.test_value(res[0].events().ebounds().emax().TeV(), 10.0,
223   - 'Check minimum energy')
  226 + 'Check maximum energy')
224 227 self.test_value(res[0].events().ebounds().size(), 5,
225 228 'Check number of energy bins')
226 229 self.test_value(res[0].events().number(), 20,
... ... @@ -263,7 +266,7 @@ class Test(test):
263 266 self.test_value(res[0].on_spec().ebounds().emin().TeV(), 1.0,
264 267 'Check minimum energy of On spectrum')
265 268 self.test_value(res[0].on_spec().ebounds().emax().TeV(), 10.0,
266   - 'Check minimum energy of On spectrum')
  269 + 'Check maximum energy of On spectrum')
267 270 self.test_value(res[0].on_spec().ebounds().size(), 5,
268 271 'Check number of energy bins of On spectrum')
269 272 self.test_value(res[0].on_spec().counts(), 24,
... ... @@ -271,7 +274,7 @@ class Test(test):
271 274 self.test_value(res[0].off_spec().ebounds().emin().TeV(), 1.0,
272 275 'Check minimum energy of Off spectrum')
273 276 self.test_value(res[0].off_spec().ebounds().emax().TeV(), 10.0,
274   - 'Check minimum energy of Off spectrum')
  277 + 'Check maximum energy of Off spectrum')
275 278 self.test_value(res[0].off_spec().ebounds().size(), 5,
276 279 'Check number of energy bins of Off spectrum')
277 280 self.test_value(res[0].off_spec().counts(), 2,
... ... @@ -279,19 +282,19 @@ class Test(test):
279 282 self.test_value(res[0].arf().ebounds().emin().TeV(), 0.5,
280 283 'Check minimum energy of ARF')
281 284 self.test_value(res[0].arf().ebounds().emax().TeV(), 12.0,
282   - 'Check minimum energy of ARF')
  285 + 'Check maximum energy of ARF')
283 286 self.test_value(res[0].arf().ebounds().size(), 41,
284 287 'Check number of energy bins of ARF')
285 288 self.test_value(res[0].rmf().etrue().emin().TeV(), 0.5,
286 289 'Check minimum true energy of RMF')
287 290 self.test_value(res[0].rmf().etrue().emax().TeV(), 12.0,
288   - 'Check minimum true energy of RMF')
  291 + 'Check maximum true energy of RMF')
289 292 self.test_value(res[0].rmf().etrue().size(), 41,
290 293 'Check number of true energy bins of RMF')
291 294 self.test_value(res[0].rmf().emeasured().emin().TeV(), 1.0,
292 295 'Check minimum reconstructed energy of RMF')
293 296 self.test_value(res[0].rmf().emeasured().emax().TeV(), 10.0,
294   - 'Check minimum reconstructed energy of RMF')
  297 + 'Check maximum reconstructed energy of RMF')
295 298 self.test_value(res[0].rmf().emeasured().size(), 5,
296 299 'Check number of reconstructed energy bins of RMF')
297 300  
... ... @@ -315,7 +318,7 @@ class Test(test):
315 318 self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
316 319 'Check minimum energy')
317 320 self.test_value(res[0].events().ebounds().emax().TeV(), 10.0,
318   - 'Check minimum energy')
  321 + 'Check maximum energy')
319 322  
320 323 # Check energy dispersion flag
321 324 self.test_assert(not (res[0].response().use_edisp()),
... ... @@ -332,7 +335,39 @@ class Test(test):
332 335 """
333 336 Test set_obs() function
334 337 """
335   - # TODO: implement
  338 + # Setup pointing direction
  339 + pnt = gammalib.GSkyDir()
  340 + pnt.radec_deg(83.63, 22.51)
  341 +
  342 + # Setup one CTA observation
  343 + res = obsutils.set_obs(pnt, emin=1.0, emax=10.0)
  344 +
  345 + # Check result
  346 + self.test_value(res.eventtype(), 'EventList', 'Check event type')
  347 + self.test_value(res.events().ebounds().emin().TeV(), 1.0,
  348 + 'Check minimum energy')
  349 + self.test_value(res.events().ebounds().emax().TeV(), 10.0,
  350 + 'Check maximum energy')
  351 + self.test_value(res.pointing().dir().ra_deg(), 83.63,
  352 + 'Check pointing Right Ascension')
  353 + self.test_value(res.pointing().dir().dec_deg(), 22.51,
  354 + 'Check pointing declination')
  355 +
  356 + # Setup one CTA observation for local caldb directory
  357 + res = obsutils.set_obs(pnt, caldb=self._datadir, irf='irf_file.fits',
  358 + emin=1.0, emax=10.0)
  359 +
  360 + # Check result
  361 + self.test_value(res.eventtype(), 'EventList', 'Check event type')
  362 + self.test_value(res.events().ebounds().emin().TeV(), 1.0,
  363 + 'Check minimum energy')
  364 + self.test_value(res.events().ebounds().emax().TeV(), 10.0,
  365 + 'Check maximum energy')
  366 + self.test_value(res.pointing().dir().ra_deg(), 83.63,
  367 + 'Check pointing Right Ascension')
  368 + self.test_value(res.pointing().dir().dec_deg(), 22.51,
  369 + 'Check pointing declination')
  370 +
336 371 # Return
337 372 return
338 373  
... ... @@ -341,7 +376,23 @@ class Test(test):
341 376 """
342 377 Test set_obs_list() function
343 378 """
344   - # TODO: implement
  379 + # Setup observation definition list
  380 + obsdeflist = [{'ra': 83.63, 'dec': 22.51}]
  381 +
  382 + # Setup observation container
  383 + res = obsutils.set_obs_list(obsdeflist)
  384 +
  385 + # Check result
  386 + self.test_value(res.size(), 1, 'Check number of observations')
  387 + self.test_value(res.models().size(), 0, 'Check number of models')
  388 + self.test_value(res.nobserved(), 0, 'Check number of observed events')
  389 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  390 + self.test_value(res[0].eventtype(), 'EventList', 'Check event type')
  391 + self.test_value(res[0].pointing().dir().ra_deg(), 83.63,
  392 + 'Check pointing Right Ascension')
  393 + self.test_value(res[0].pointing().dir().dec_deg(), 22.51,
  394 + 'Check pointing declination')
  395 +
345 396 # Return
346 397 return
347 398  
... ... @@ -360,7 +411,7 @@ class Test(test):
360 411 self.test_value(len(res), 4, 'Check number of observations for "four"'
361 412 ' pattern')
362 413  
363   - # Get invalud observations pattern
  414 + # Get invalid observations pattern
364 415 self.test_try('Check invalid pointing pattern')
365 416 try:
366 417 res = obsutils.set_obs_patterns('five')
... ... @@ -376,7 +427,28 @@ class Test(test):
376 427 """
377 428 Test set_observations() function
378 429 """
379   - # TODO: implement
  430 + # Set observation
  431 + res = obsutils.set_observations(83.63, 22.51, 5.0,
  432 + 0.0, 1800.0,
  433 + 1.0, 100.0,
  434 + 'South_50h', 'prod2',
  435 + pattern='single')
  436 +
  437 + # Check result
  438 + self.test_value(res.size(), 1, 'Check number of observations')
  439 + self.test_value(res.models().size(), 0, 'Check number of models')
  440 + self.test_value(res.nobserved(), 0, 'Check number of observed events')
  441 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  442 + self.test_value(res[0].eventtype(), 'EventList', 'Check event type')
  443 + self.test_value(res[0].events().ebounds().emin().TeV(), 1.0,
  444 + 'Check minimum energy')
  445 + self.test_value(res[0].events().ebounds().emax().TeV(), 100.0,
  446 + 'Check maximum energy')
  447 + self.test_value(res[0].pointing().dir().ra_deg(), 83.63,
  448 + 'Check pointing Right Ascension')
  449 + self.test_value(res[0].pointing().dir().dec_deg(), 22.51,
  450 + 'Check pointing declination')
  451 +
380 452 # Return
381 453 return
382 454  
... ... @@ -385,7 +457,97 @@ class Test(test):
385 457 """
386 458 Test get_stacked_response() function
387 459 """
388   - # TODO: implement
  460 + # Set-up observation container
  461 + obs = self._setup_sim()
  462 +
  463 + # Get stacked response with xref/yref set
  464 + res = obsutils.get_stacked_response(obs, 83.6331, 22.0145,
  465 + coordsys='CEL',
  466 + binsz=0.1, nxpix=10, nypix=10,
  467 + emin=1.0, emax=100.0, enumbins=5,
  468 + edisp=False)
  469 +
  470 + # Check result
  471 + self.test_value(res['expcube'].cube().npix(), 100,
  472 + 'Check number of exposure cube pixels')
  473 + self.test_value(res['expcube'].cube().nmaps(), 61,
  474 + 'Check number of exposure cube maps')
  475 + self.test_value(res['psfcube'].cube().npix(), 4,
  476 + 'Check number of PSF cube pixels')
  477 + self.test_value(res['psfcube'].cube().nmaps(), 12200,
  478 + 'Check number of PSF cube maps')
  479 + self.test_value(res['bkgcube'].cube().npix(), 100,
  480 + 'Check number of background cube pixels')
  481 + self.test_value(res['bkgcube'].cube().nmaps(), 61,
  482 + 'Check number of background cube maps')
  483 +
  484 + # Get stacked response with None xref/yref
  485 + res = obsutils.get_stacked_response(obs, None, None,
  486 + coordsys='CEL',
  487 + binsz=0.1, nxpix=10, nypix=10,
  488 + emin=1.0, emax=100.0, enumbins=5,
  489 + edisp=False)
  490 +
  491 + # Check result
  492 + self.test_value(res['expcube'].cube().npix(), 100,
  493 + 'Check number of exposure cube pixels')
  494 + self.test_value(res['expcube'].cube().nmaps(), 61,
  495 + 'Check number of exposure cube maps')
  496 + self.test_value(res['psfcube'].cube().npix(), 4,
  497 + 'Check number of PSF cube pixels')
  498 + self.test_value(res['psfcube'].cube().nmaps(), 12200,
  499 + 'Check number of PSF cube maps')
  500 + self.test_value(res['bkgcube'].cube().npix(), 100,
  501 + 'Check number of background cube pixels')
  502 + self.test_value(res['bkgcube'].cube().nmaps(), 61,
  503 + 'Check number of background cube maps')
  504 +
  505 + # Get stacked response with xref/yref set and large number of energy bins
  506 + res = obsutils.get_stacked_response(obs, 83.6331, 22.0145,
  507 + coordsys='CEL',
  508 + binsz=0.1, nxpix=10, nypix=10,
  509 + emin=1.0, emax=100.0, enumbins=100,
  510 + edisp=False)
  511 +
  512 + # Check result
  513 + self.test_value(res['expcube'].cube().npix(), 100,
  514 + 'Check number of exposure cube pixels')
  515 + self.test_value(res['expcube'].cube().nmaps(), 101,
  516 + 'Check number of exposure cube maps')
  517 + self.test_value(res['psfcube'].cube().npix(), 4,
  518 + 'Check number of PSF cube pixels')
  519 + self.test_value(res['psfcube'].cube().nmaps(), 20200,
  520 + 'Check number of PSF cube maps')
  521 + self.test_value(res['bkgcube'].cube().npix(), 100,
  522 + 'Check number of background cube pixels')
  523 + self.test_value(res['bkgcube'].cube().nmaps(), 101,
  524 + 'Check number of background cube maps')
  525 +
  526 + # Get stacked response with energy dispersion
  527 + res = obsutils.get_stacked_response(obs, 83.6331, 22.0145,
  528 + coordsys='CEL',
  529 + binsz=0.1, nxpix=10, nypix=10,
  530 + emin=1.0, emax=100.0, enumbins=5,
  531 + edisp=True)
  532 +
  533 + # Check result
  534 + self.test_value(res['expcube'].cube().npix(), 100,
  535 + 'Check number of exposure cube pixels')
  536 + self.test_value(res['expcube'].cube().nmaps(), 61,
  537 + 'Check number of exposure cube maps')
  538 + self.test_value(res['psfcube'].cube().npix(), 4,
  539 + 'Check number of PSF cube pixels')
  540 + self.test_value(res['psfcube'].cube().nmaps(), 12200,
  541 + 'Check number of PSF cube maps')
  542 + self.test_value(res['bkgcube'].cube().npix(), 100,
  543 + 'Check number of background cube pixels')
  544 + self.test_value(res['bkgcube'].cube().nmaps(), 61,
  545 + 'Check number of background cube maps')
  546 + self.test_value(res['edispcube'].cube().npix(), 4,
  547 + 'Check number of energy dispersion cube pixels')
  548 + self.test_value(res['edispcube'].cube().nmaps(), 6100,
  549 + 'Check number of energy dispersion cube maps')
  550 +
389 551 # Return
390 552 return
391 553  
... ... @@ -394,7 +556,135 @@ class Test(test):
394 556 """
395 557 Test get_stacked_obs() function
396 558 """
397   - # TODO: implement
  559 + # Set-up unbinned cslightcrv (is not run here!!!)
  560 + lcrv = cscripts.cslightcrv()
  561 + lcrv['inobs'] = self._events
  562 + lcrv['inmodel'] = self._model
  563 + lcrv['srcname'] = 'Crab'
  564 + lcrv['caldb'] = self._caldb
  565 + lcrv['irf'] = self._irf
  566 + lcrv['tbinalg'] = 'LIN'
  567 + lcrv['tmin'] = '2020-01-01T00:00:00'
  568 + lcrv['tmax'] = '2020-01-01T00:05:00'
  569 + lcrv['tbins'] = 2
  570 + lcrv['method'] = '3D'
  571 + lcrv['emin'] = 1.0
  572 + lcrv['emax'] = 100.0
  573 + lcrv['enumbins'] = 3
  574 + lcrv['coordsys'] = 'CEL'
  575 + lcrv['proj'] = 'TAN'
  576 + lcrv['xref'] = 83.63
  577 + lcrv['yref'] = 22.01
  578 + lcrv['binsz'] = 0.1
  579 + lcrv['nxpix'] = 10
  580 + lcrv['nypix'] = 10
  581 + lcrv['outfile'] = 'obsutils_cslightcrv_py1.fits'
  582 + lcrv['logfile'] = 'obsutils_cslightcrv_py1.log'
  583 + lcrv._get_parameters()
  584 +
  585 + # Get stacked observation container
  586 + res = obsutils.get_stacked_obs(lcrv, lcrv.obs())
  587 +
  588 + # Apply energy dispersion
  589 + res[0].response().apply_edisp(True)
  590 +
  591 + # Check result
  592 + self.test_value(res.size(), 1, 'Check number of observations')
  593 + self.test_value(res.models().size(), 2, 'Check number of models')
  594 + self.test_value(res.nobserved(), 96, 'Check number of observed events')
  595 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  596 + self.test_assert(not res[0].response().use_edisp(),
  597 + 'Check that no energy dispersion is used')
  598 +
  599 + # Set energy dispersion
  600 + lcrv['edisp'] = True
  601 +
  602 + # Get stacked observation container
  603 + res = obsutils.get_stacked_obs(lcrv, lcrv.obs())
  604 +
  605 + # Apply energy dispersion
  606 + res[0].response().apply_edisp(True)
  607 +
  608 + # Check result
  609 + self.test_value(res.size(), 1, 'Check number of observations')
  610 + self.test_value(res.models().size(), 2, 'Check number of models')
  611 + self.test_value(res.nobserved(), 96, 'Check number of observed events')
  612 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  613 + self.test_assert(res[0].response().use_edisp(),
  614 + 'Check that energy dispersion is used')
  615 +
  616 + # Increase chatter level
  617 + lcrv['edisp'] = False
  618 + lcrv['chatter'] = 4
  619 +
  620 + # Get stacked observation container
  621 + res = obsutils.get_stacked_obs(lcrv, lcrv.obs())
  622 +
  623 + # Apply energy dispersion
  624 + res[0].response().apply_edisp(True)
  625 +
  626 + # Check result
  627 + self.test_value(res.size(), 1, 'Check number of observations')
  628 + self.test_value(res.models().size(), 2, 'Check number of models')
  629 + self.test_value(res.nobserved(), 96, 'Check number of observed events')
  630 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  631 + self.test_assert(not res[0].response().use_edisp(),
  632 + 'Check that no energy dispersion is used')
  633 +
398 634 # Return
399 635 return
400 636  
  637 + # Test get_onoff_obs() function
  638 + def _test_get_onoff_obs(self):
  639 + """
  640 + Test get_onoff_obs() function
  641 + """
  642 + # Set-up unbinned cslightcrv (is not run here!!!)
  643 + lcrv = cscripts.cslightcrv()
  644 + lcrv['inobs'] = self._events
  645 + lcrv['inmodel'] = self._model
  646 + lcrv['srcname'] = 'Crab'
  647 + lcrv['caldb'] = self._caldb
  648 + lcrv['irf'] = self._irf
  649 + lcrv['tbinalg'] = 'LIN'
  650 + lcrv['tmin'] = '2020-01-01T00:00:00'
  651 + lcrv['tmax'] = '2020-01-01T00:05:00'
  652 + lcrv['tbins'] = 2
  653 + lcrv['method'] = 'ONOFF'
  654 + lcrv['emin'] = 1.0
  655 + lcrv['emax'] = 100.0
  656 + lcrv['enumbins'] = 2
  657 + lcrv['coordsys'] = 'CEL'
  658 + lcrv['xref'] = 83.63
  659 + lcrv['yref'] = 22.01
  660 + lcrv['rad'] = 0.2
  661 + lcrv['outfile'] = 'obsutils_cslightcrv_py2.fits'
  662 + lcrv['logfile'] = 'obsutils_cslightcrv_py2.log'
  663 + lcrv._get_parameters()
  664 +
  665 + # Get on/off observation container
  666 + res = obsutils.get_onoff_obs(lcrv, lcrv.obs())
  667 +
  668 + # Check result
  669 + self.test_value(res.size(), 1, 'Check number of observations')
  670 + self.test_value(res.models().size(), 2, 'Check number of models')
  671 + self.test_value(res.nobserved(), 91, 'Check number of observed events')
  672 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  673 +
  674 + # Use galactic coordinates
  675 + lcrv['coordsys'] = 'GAL'
  676 + lcrv['xref'] = 184.5597
  677 + lcrv['yref'] = -5.7892
  678 + lcrv['chatter'] = 4
  679 +
  680 + # Get on/off observation container
  681 + res = obsutils.get_onoff_obs(lcrv, lcrv.obs())
  682 +
  683 + # Check result
  684 + self.test_value(res.size(), 1, 'Check number of observations')
  685 + self.test_value(res.models().size(), 2, 'Check number of models')
  686 + self.test_value(res.nobserved(), 91, 'Check number of observed events')
  687 + self.test_value(res.npred(), 0.0, 'Check number of predicted events')
  688 +
  689 + # Return
  690 + return
... ...