Andreas Specovius / ctools

Commit 9600793a6504cb003846a40731fec91b2d33b0df

Authored by Luigi Tibaldo
1 parent e231d0c0

added test unit

Inline Side-by-side
Showing 3 changed files with 507 additions and 106 deletions
cscripts/csscs.py
  View file @9600793
... ... @@ -18,12 +18,12 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
21   -import sys
22   -import gammalib
23 21 import ctools
  22 +import gammalib
  23 +import math
  24 +import sys
24 25 from cscripts import mputils
25 26 from cscripts import obsutils
26   -import math
27 27  
28 28  
29 29 # ================= #
... ... @@ -53,7 +53,8 @@ class csscs(ctools.csobservation):
53 53 self._srcnames = []
54 54 self._method = None
55 55 self._fits = None
56   - self._excl_reg_map = None # Exclusion region map for on/off analysis
  56 + self._excl_reg_map = None # Exclusion region map for on/off analysis
  57 + self._excl_reg_name = None
57 58  
58 59 # Return
59 60 return
... ... @@ -70,12 +71,13 @@ class csscs(ctools.csobservation):
70 71 """
71 72 # Set pickled dictionary
72 73 # Set pickled dictionary
73   - state = {'base' : ctools.csobservation.__getstate__(self),
74   - 'nthreads' : self._nthreads,
75   - 'method' : self._method,
76   - 'srcnames' : self._srcnames,
77   - 'fits' : self._fits,
78   - 'excl_reg_map' : self._excl_reg_map}
  74 + state = {'base': ctools.csobservation.__getstate__(self),
  75 + 'nthreads': self._nthreads,
  76 + 'method': self._method,
  77 + 'srcnames': self._srcnames,
  78 + 'fits': self._fits,
  79 + 'excl_reg_map': self._excl_reg_map,
  80 + 'excl_reg_name' : self._excl_reg_name}
79 81  
80 82 # Return pickled dictionary
81 83 return state
... ... @@ -96,11 +98,11 @@ class csscs(ctools.csobservation):
96 98 self._method = state['method']
97 99 self._fits = state['fits']
98 100 self._excl_reg_map = state['excl_reg_map']
  101 + self._excl_reg_name = state['excl_reg_name']
99 102  
100 103 # Return
101 104 return
102 105  
103   -
104 106 # Private methods
105 107 def _get_onoff_parameters(self):
106 108 """
... ... @@ -117,8 +119,8 @@ class csscs(ctools.csobservation):
117 119 # Exclusion map set and is not empty
118 120 pass
119 121 elif self['inexclusion'].is_valid():
120   - inexclusion = self['inexclusion'].filename()
121   - self._excl_reg_map = gammalib.GSkyRegionMap(inexclusion)
  122 + self._excl_reg_name = self['inexclusion'].filename()
  123 + self._excl_reg_map = gammalib.GSkyRegionMap(self._excl_reg_name)
122 124 else:
123 125 msg = 'csscs in On/Off mode requires input exclusion region.'
124 126 raise RuntimeError(msg)
... ... @@ -150,8 +152,8 @@ class csscs(ctools.csobservation):
150 152 # Collect number of unbinned, binned and On/Off observations in
151 153 # observation container
152 154 n_unbinned = 0
153   - n_binned = 0
154   - n_onoff = 0
  155 + n_binned = 0
  156 + n_onoff = 0
155 157 for obs in self.obs():
156 158 if obs.classname() == 'GCTAObservation':
157 159 if obs.eventtype() == 'CountsCube':
... ... @@ -160,16 +162,16 @@ class csscs(ctools.csobservation):
160 162 n_unbinned += 1
161 163 elif obs.classname() == 'GCTAOnOffObservation':
162 164 n_onoff += 1
163   - n_cta = n_unbinned + n_binned + n_onoff
  165 + n_cta = n_unbinned + n_binned + n_onoff
164 166 n_other = self.obs().size() - n_cta
165 167  
166 168 # Check if there are On/Off or non-IACT observations
167 169 if n_onoff > 0 or n_other > 0:
168   - msg = 'On/Off or non-CTA observations found. '\
  170 + msg = 'On/Off or non-CTA observations found. ' \
169 171 'csscs does not support this type of observations.'
170 172 raise RuntimeError(msg)
171 173 # Otherwise check if there is a mix of binned and unbinned
172   - elif n_unbinned > 0 and n_binned>0:
  174 + elif n_unbinned > 0 and n_binned > 0:
173 175 msg = 'Mix of unbinned and binned CTA observations ' \
174 176 'found in observation container. csscs does not ' \
175 177 'support this mix.'
... ... @@ -228,7 +230,7 @@ class csscs(ctools.csobservation):
228 230 # Loop over models and verify if model is of type GModelSky
229 231 for model in self.obs().models():
230 232 if model.classname() == 'GModelSky':
231   - nsources +=1
  233 + nsources += 1
232 234 # If there are background gamma-ray sources in the model
233 235 # throw runtime error
234 236 if nsources > len(self._srcnames):
... ... @@ -327,7 +329,7 @@ class csscs(ctools.csobservation):
327 329 # Free the normalisation parameter which is assumed to be
328 330 # the first spectral parameter
329 331 if normpar.is_fixed():
330   - self._log_string(gammalib.EXPLICIT, ' Freeing "'+normpar.name()+'"')
  332 + self._log_string(gammalib.EXPLICIT, ' Freeing "' + normpar.name() + '"')
331 333 normpar.free()
332 334  
333 335 # Optionally compute Test Statistic value
... ... @@ -338,7 +340,7 @@ class csscs(ctools.csobservation):
338 340 elif self['fix_bkg'].boolean() and not model.classname() == 'GModelSky':
339 341  
340 342 if normpar.is_free():
341   - self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
  343 + self._log_string(gammalib.EXPLICIT, ' Fixing "' + normpar.name() + '"')
342 344 normpar.fix()
343 345  
344 346  
... ... @@ -346,13 +348,13 @@ class csscs(ctools.csobservation):
346 348 elif self['fix_srcs'].boolean() and model.classname() == 'GModelSky':
347 349  
348 350 if normpar.is_free():
349   - self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
  351 + self._log_string(gammalib.EXPLICIT, ' Fixing "' + normpar.name() + '"')
350 352 normpar.fix()
351 353  
352 354 # Return
353 355 return
354 356  
355   - def _mask_cube(self,obs,ra,dec,rad):
  357 + def _mask_cube(self, obs, ra, dec, rad):
356 358 """
357 359 Mask cube observation
358 360  
... ... @@ -392,7 +394,7 @@ class csscs(ctools.csobservation):
392 394 # Return
393 395 return new_obs
394 396  
395   - def _mask_evlist(self,obs,ra,dec,rad):
  397 + def _mask_evlist(self, obs, ra, dec, rad):
396 398 """
397 399 Mask event list
398 400  
... ... @@ -434,7 +436,7 @@ class csscs(ctools.csobservation):
434 436 # Return
435 437 return new_obs
436 438  
437   - def _mask_observations(self,ra,dec,rad):
  439 + def _mask_observations(self, ra, dec, rad):
438 440 """
439 441 Create observations restricted to circular ROI
440 442  
... ... @@ -457,14 +459,14 @@ class csscs(ctools.csobservation):
457 459 elif self._method == 'BINNED':
458 460 new_obs = self._mask_cube(self.obs(), ra, dec, rad)
459 461 elif self._method == 'ONOFF':
460   - new_obs = obsutils.get_onoff_obs(self,self.obs(),nthreads=1,
461   - ra = ra, dec = dec,
462   - srcname = self._srcnames[0])
  462 + new_obs = obsutils.get_onoff_obs(self, self.obs(), nthreads=1,
  463 + ra=ra, dec=dec,
  464 + srcname=self._srcnames[0])
463 465  
464 466 # Return
465 467 return new_obs
466 468  
467   - def _pixel_analysis(self,inx):
  469 + def _pixel_analysis(self, inx):
468 470 """
469 471 Performs analysis over the region of interest
470 472 corresponding to a single pixel of output map
... ... @@ -487,25 +489,25 @@ class csscs(ctools.csobservation):
487 489 dec = pixdir.dec_deg()
488 490  
489 491 # Write header for spatial bin
490   - msg = 'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx,ra,dec)
  492 + msg = 'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx, ra, dec)
491 493 self._log_header2(gammalib.TERSE, msg)
492 494  
493 495 # Initialize results
494 496 result = {}
495 497 for name in self._srcnames:
496   - result[name] = {'flux' : 0.0,
  498 + result[name] = {'flux': 0.0,
497 499 'flux_err': 0.0,
498   - 'TS' : 0.0,
499   - 'ulimit' : -1.0}
  500 + 'TS': 0.0,
  501 + 'ulimit': -1.0}
500 502  
501 503 # Mask observations
502 504 self._log_header3(gammalib.NORMAL, 'Masking observations')
503   - masked_obs = self._mask_observations(ra,dec,self['rad'].real())
  505 + masked_obs = self._mask_observations(ra, dec, self['rad'].real())
504 506  
505 507 # Set up likelihood analysis
506 508 self._log_header3(gammalib.NORMAL, 'Performing fit in region')
507 509 like = ctools.ctlike(masked_obs)
508   - like['edisp'] = self['edisp'].boolean()
  510 + like['edisp'] = self['edisp'].boolean()
509 511 like['nthreads'] = 1 # Avoids OpenMP conflict
510 512  
511 513 # If chatter level is verbose and debugging is requested then
... ... @@ -522,10 +524,10 @@ class csscs(ctools.csobservation):
522 524 # Prepare objects for flux extraction
523 525 # ROI
524 526 centre = gammalib.GSkyDir()
525   - centre.radec_deg(ra,dec)
526   - roi = gammalib.GSkyRegionCircle(centre,self['rad'].real())
  527 + centre.radec_deg(ra, dec)
  528 + roi = gammalib.GSkyRegionCircle(centre, self['rad'].real())
527 529 # Energy boundaries
528   - emin = gammalib.GEnergy(self['emin'].real(),'TeV')
  530 + emin = gammalib.GEnergy(self['emin'].real(), 'TeV')
529 531 emax = gammalib.GEnergy(self['emax'].real(), 'TeV')
530 532  
531 533 # Continue only if log-likelihood is non-zero
... ... @@ -539,11 +541,12 @@ class csscs(ctools.csobservation):
539 541  
540 542 # Get flux
541 543 # Integrate spectral model between emin and emax
542   - flux = source.spectral().flux(emin,emax)
  544 + flux = source.spectral().flux(emin, emax)
543 545 # Calculate correction factor
544 546 # Spatial model flux over ROI divided by ROI solid angle
545 547 corr_factor = source.spatial().flux(roi)
546   - corr_factor /= gammalib.twopi * (1 - math.cos(math.radians(self['rad'].real())))
  548 + corr_factor /= gammalib.twopi * (
  549 + 1 - math.cos(math.radians(self['rad'].real())))
547 550 # Multiply flux by correction factor
548 551 flux *= corr_factor
549 552 result[name]['flux'] = flux
... ... @@ -555,7 +558,7 @@ class csscs(ctools.csobservation):
555 558 # Relative error on flux is same as on normalization
556 559 # Avoid zero division error
557 560 if normpar.value() > 0.:
558   - flux_error = flux * normpar.error()/normpar.value()
  561 + flux_error = flux * normpar.error() / normpar.value()
559 562 else:
560 563 flux_error = 0.
561 564 result[name]['flux_error'] = flux_error
... ... @@ -591,58 +594,57 @@ class csscs(ctools.csobservation):
591 594 result[name]['ulimit'] = ulimit_value
592 595 except:
593 596 self._log_string(gammalib.NORMAL, 'Upper limit '
594   - 'calculation failed.')
  597 + 'calculation failed.')
595 598  
596 599  
597 600 else:
598 601 value = 'Likelihood is zero. Bin is skipped.'
599   - self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra,dec), value)
  602 + self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra, dec), value)
600 603  
601 604 # Return result
602 605 return result
603 606  
604   - def _write_hdu_keywords(self,hdu, is_flux = False):
  607 + def _write_hdu_keywords(self, hdu, is_flux=False):
605 608 """
606 609 Append cards to HDU
607 610  
608 611 Parameters
609 612 ----------
610   - hdu : `dict` HDU
  613 + hdu : `~gammalib.GFitsHDU` HDU
611 614 """
612 615  
613 616 # Flux units
614 617 if is_flux:
615   - hdu.card('BUNIT', 'ph/cm2/s/sr' , 'Photon flux')
  618 + hdu.card('BUNIT', 'ph/cm2/s/sr', 'Photon flux')
616 619  
617 620 # Minimum and maximum energy
618   - hdu.card("E_MIN", self['emin'].real(), "[TeV] Lower energy boundary")
619   - hdu.card("E_MAX", self['emax'].real(), "[TeV] Upper energy boundary")
620   - hdu.card("EUNIT", "TeV", "Units for E_MIN and E_MAX")
  621 + hdu.card("E_MIN", self['emin'].real(), "[TeV] Lower energy boundary")
  622 + hdu.card("E_MAX", self['emax'].real(), "[TeV] Upper energy boundary")
  623 + hdu.card("EUNIT", "TeV", "Units for E_MIN and E_MAX")
621 624  
622 625 # ROI size
623   - hdu.card('ROISZ', self['rad'].real() ,
  626 + hdu.card('ROISZ', self['rad'].real(),
624 627 '[deg] Size of ROI used for component separation')
625 628  
626 629 # Analysis method
627   - hdu.card('METHOD',self._method, 'Analysis method')
  630 + hdu.card('METHOD', self._method, 'Analysis method')
628 631  
629 632 # For On/Off log information
630 633 if self._method == 'ONOFF':
631 634  
632 635 # Exclusion map
633   - if (self._excl_reg_map is not None) and self['inexclusion'].is_valid():
634   - hdu.card("EXCLMAP", self['inexclusion'].filename().url(), "Exclusion map name")
  636 + if self._excl_reg_name is not None:
  637 + hdu.card("EXCLMAP", self._excl_reg_name.url(), "Exclusion map name")
635 638  
636 639 # Model background?
637   - hdu.card("BKGMODEL", self['use_model_bkg'].boolean(),'Use background model?')
  640 + hdu.card("BKGMODEL", self['use_model_bkg'].boolean(), 'Use background model?')
638 641  
639 642 # Fit statistic
640   - hdu.card("STAT", self['statistic'].string(),'Fit statistic')
  643 + hdu.card("STAT", self['statistic'].string(), 'Fit statistic')
641 644  
642 645 return
643 646  
644   -
645   - def _fill_fits(self,results):
  647 + def _fill_fits(self, results):
646 648 """
647 649 Fill FITS object to store the results
648 650  
... ... @@ -674,27 +676,29 @@ class csscs(ctools.csobservation):
674 676  
675 677 # Loop over pixels and fill maps
676 678 for inx in range(fmap.npix()):
677   -
678 679 # Fill maps
679   - fmap[inx,0] = results[inx][name]['flux']
  680 + fmap[inx, 0] = results[inx][name]['flux']
680 681 errmap[inx, 0] = results[inx][name]['flux_error']
681 682 tsmap[inx, 0] = results[inx][name]['TS']
682 683 ulmap[inx, 0] = results[inx][name]['ulimit']
683 684  
  685 + # Uppercase name to adhere to gammalib convention
  686 + Name = gammalib.toupper(name)
  687 +
684 688 # Write maps to Fits
685   - fhdu = fmap.write(fits, name + ' FLUX')
686   - self._write_hdu_keywords(fhdu,is_flux=True)
687   - errhdu = errmap.write(fits, name + ' FLUX ERROR')
  689 + fhdu = fmap.write(fits, Name + ' FLUX')
  690 + self._write_hdu_keywords(fhdu, is_flux=True)
  691 + errhdu = errmap.write(fits, Name + ' FLUX ERROR')
688 692 self._write_hdu_keywords(errhdu, is_flux=True)
689 693  
690 694 # If requested write TS map to fits
691 695 if self['calc_ts'].boolean():
692   - tshdu = tsmap.write(fits, name + ' TS')
  696 + tshdu = tsmap.write(fits, Name + ' TS')
693 697 self._write_hdu_keywords(tshdu)
694 698  
695 699 # If requested write upper limit map to fits
696 700 if self['calc_ulim'].boolean():
697   - ulhdu = ulmap.write(fits, name + ' FLUX UPPER LIMIT')
  701 + ulhdu = ulmap.write(fits, Name + ' FLUX UPPER LIMIT')
698 702 self._write_hdu_keywords(ulhdu, is_flux=True)
699 703  
700 704 # Return
... ... @@ -725,11 +729,13 @@ class csscs(ctools.csobservation):
725 729 skymap = None
726 730  
727 731 # Proceed only if fits has been filled
728   - if self._fits != None:
  732 + if self._fits != None:
729 733 # Check that the name is in the list of target sources
730 734 if name in self._srcnames:
  735 + # Uppercase name
  736 + Name = gammalib.toupper(name)
731 737 # Assemble HDU name
732   - hduname = name + ' ' + quantity
  738 + hduname = Name + ' ' + quantity
733 739 # Try to fetch HDU
734 740 try:
735 741 hdu = self._fits[hduname]
... ... @@ -747,7 +753,6 @@ class csscs(ctools.csobservation):
747 753 # Return
748 754 return skymap
749 755  
750   -
751 756 # Public methods
752 757 def run(self):
753 758 """
... ... @@ -773,8 +778,8 @@ class csscs(ctools.csobservation):
773 778 if self._nthreads > 1:
774 779  
775 780 # Compute values in pixels
776   - args = [(self, '_pixel_analysis', i)
777   - for i in range(npix)]
  781 + args = [(self, '_pixel_analysis', i)
  782 + for i in range(npix)]
778 783 poolresults = mputils.process(self._nthreads, mputils.mpfunc, args)
779 784  
780 785 # Construct results
... ... @@ -804,7 +809,6 @@ class csscs(ctools.csobservation):
804 809  
805 810 # Continue only if FITS file is valid
806 811 if self._fits != None:
807   -
808 812 # Get outmap parameter
809 813 outfile = self['outfile'].filename()
810 814  
... ... @@ -850,7 +854,7 @@ class csscs(ctools.csobservation):
850 854 # Return
851 855 return self._fits
852 856  
853   - def flux(self,name):
  857 + def flux(self, name):
854 858 """
855 859 Return flux skymap
856 860  
... ... @@ -865,11 +869,11 @@ class csscs(ctools.csobservation):
865 869 Flux skymap
866 870 """
867 871  
868   - skymap = self._get_skymap(name,'FLUX')
  872 + skymap = self._get_skymap(name, 'FLUX')
869 873  
870 874 return skymap
871 875  
872   - def flux_error(self,name):
  876 + def flux_error(self, name):
873 877 """
874 878 Return flux error skymap
875 879  
... ... @@ -884,11 +888,11 @@ class csscs(ctools.csobservation):
884 888 Flux error skymap
885 889 """
886 890  
887   - skymap = self._get_skymap(name,'FLUX ERROR')
  891 + skymap = self._get_skymap(name, 'FLUX ERROR')
888 892  
889 893 return skymap
890 894  
891   - def ts(self,name):
  895 + def ts(self, name):
892 896 """
893 897 Return TS skymap
894 898  
... ... @@ -905,7 +909,7 @@ class csscs(ctools.csobservation):
905 909  
906 910 # Check that TS computation is requested
907 911 if self['calc_ts'].boolean():
908   - skymap = self._get_skymap(name,'TS')
  912 + skymap = self._get_skymap(name, 'TS')
909 913 # Otherwise throw error
910 914 else:
911 915 msg = 'TS computation not requested. ' \
... ... @@ -915,7 +919,7 @@ class csscs(ctools.csobservation):
915 919  
916 920 return skymap
917 921  
918   - def ulimit(self,name):
  922 + def ulimit(self, name):
919 923 """
920 924 Return flux upper limit skymap
921 925  
... ... @@ -932,7 +936,7 @@ class csscs(ctools.csobservation):
932 936  
933 937 # Check that upper limit computation is requested
934 938 if self['calc_ulim'].boolean():
935   - skymap = self._get_skymap(name,'FLUX UPPER LIMIT')
  939 + skymap = self._get_skymap(name, 'FLUX UPPER LIMIT')
936 940 # Otherwise throw error
937 941 else:
938 942 msg = 'Upper limit computation not requested. ' \
... ... @@ -940,7 +944,6 @@ class csscs(ctools.csobservation):
940 944 'to obtain upper limit.'
941 945 raise RuntimeError(msg)
942 946  
943   -
944 947 return skymap
945 948  
946 949  
... ... @@ -948,7 +951,6 @@ class csscs(ctools.csobservation):
948 951 # Main routine entry point #
949 952 # ======================== #
950 953 if __name__ == '__main__':
951   -
952 954 # Create instance of application
953 955 app = csscs(sys.argv)
954 956  
... ...
cscripts/obsutils.py
  View file @9600793
... ... @@ -775,18 +775,14 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = &#39;&#39;):
775 775  
776 776 # Initialise inmodel, inexclusion and use_model_bkg
777 777 inmodel = 'NONE'
778   - inexclusion = 'NONE'
779 778 use_model_bkg = True
780 779  
781   - # Set inmodel, srcname, inexclusion and use_model_bkg if they are available
  780 + # Set inmodel, srcname, and use_model_bkg if they are available
782 781 if cls.has_par('inmodel') and obs.models().size() == 0:
783 782 if cls['inmodel'].is_valid():
784 783 inmodel = cls['inmodel'].value()
785 784 if srcname == '' and cls.has_par('srcname'):
786 785 srcname = cls['srcname'].value()
787   - if cls.has_par('inexclusion'):
788   - if cls['inexclusion'].is_valid():
789   - inexclusion = cls['inexclusion'].value()
790 786 if cls.has_par('use_model_bkg'):
791 787 use_model_bkg = cls['use_model_bkg'].boolean()
792 788  
... ... @@ -794,7 +790,6 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = &#39;&#39;):
794 790 phagen = cscripts.csphagen(obs)
795 791 phagen['inmodel'] = inmodel
796 792 phagen['srcname'] = srcname
797   - phagen['inexclusion'] = inexclusion
798 793 phagen['emin'] = cls['emin'].real()
799 794 phagen['emax'] = cls['emax'].real()
800 795 phagen['enumbins'] = cls['enumbins'].integer()
... ... @@ -830,13 +825,21 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = &#39;&#39;):
830 825 phagen['debug'] = cls['debug'].boolean()
831 826 phagen['nthreads'] = nthreads
832 827  
833   - # Pipe exclusion map
  828 + # Pipe exclusion map in memory
834 829 if hasattr(cls, 'exclusion_map'):
835 830 exclusion_map = cls.exclusion_map()
836 831 if exclusion_map is not None:
837 832 phagen.exclusion_map(exclusion_map)
  833 + # Otherwise try to read it as user parameter
  834 + elif cls.has_par('inexclusion'):
  835 + if cls['inexclusion'].is_valid():
  836 + inexclusion = cls['inexclusion'].value()
  837 + phagen['inexclusion'] = inexclusion
  838 + # Otherwise set to None
  839 + else:
  840 + phagen['inexclusion'] = 'NONE'
838 841  
839   - # Run csphagen
  842 + # Run csphagen
840 843 phagen.run()
841 844  
842 845 # Clone resulting observation container
... ...
test/test_csscs.py
  View file @9600793
... ... @@ -18,9 +18,9 @@
18 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 19 #
20 20 # ==========================================================================
21   -import os
22   -import gammalib
23 21 import cscripts
  22 +import gammalib
  23 +import os
24 24 from testing import test
25 25  
26 26  
... ... @@ -43,6 +43,15 @@ class Test(test):
43 43 # Call base class constructor
44 44 test.__init__(self)
45 45  
  46 + # Set test datasets and parameters
  47 + self._evs_offaxis = self._datadir + '/crab_offaxis1.fits'
  48 + self._exclusion = self._datadir + '/crab_exclusion.fits'
  49 + self._expcube = self._datadir + '/crab_expcube.fits'
  50 + self._psfcube = self._datadir + '/crab_psfcube.fits'
  51 + self._bkgcube = self._datadir + '/crab_bkgcube.fits'
  52 + self._bkgmodel = self._datadir + '/crab_bkgcube.xml'
  53 + self._edispcube = self._datadir + '/crab_edispcube.fits'
  54 +
46 55 # Return
47 56 return
48 57  
... ... @@ -71,11 +80,38 @@ class Test(test):
71 80 csscs = self._script('csscs')
72 81  
73 82 # Setup csscs command
74   - cmd = csscs+' logfile="csscs_cmd1.log" chatter=1'
  83 + cmd = csscs + ' inobs="' + self._cntcube + '"' + \
  84 + ' expcube="' + self._expcube + '"' + \
  85 + ' psfcube="' + self._psfcube + '"' + \
  86 + ' bkgcube="' + self._bkgcube + '"' + \
  87 + ' inmodel="' + self._bkgmodel + '"' + \
  88 + ' srcnames="Crab"' + \
  89 + ' emin=1. emax=100.0 nxpix=2 nypix=2' + \
  90 + ' binsz=0.1 coordsys="CEL" proj="CAR"' + \
  91 + ' xref=83.63 yref=22.01 rad=0.2' + \
  92 + ' outfile="csscs_cmd1.fits"' + \
  93 + ' logfile="csscs_cmd1.log" chatter=1'
75 94  
76 95 # Check if execution was successful
77 96 self.test_value(self._execute(cmd), 0,
78   - 'Check successful execution from command line')
  97 + 'Check successful execution from command line')
  98 +
  99 + # Setup csscs command
  100 + cmd = csscs + ' inobs="' + self._cntcube + '"' + \
  101 + ' expcube="' + self._expcube + '"' + \
  102 + ' psfcube="' + self._psfcube + '"' + \
  103 + ' bkgcube="' + self._bkgcube + '"' + \
  104 + ' inmodel="' + self._bkgmodel + '"' + \
  105 + ' srcnames="Crab;source_that_is_not_in_the_model"' + \
  106 + ' emin=1. emax=100.0 nxpix=2 nypix=2' + \
  107 + ' binsz=0.1 coordsys="CEL" proj="CAR"' + \
  108 + ' xref=83.63 yref=22.01 rad=0.2' + \
  109 + ' outfile="csscs_cmd1.fits"' + \
  110 + ' logfile="csscs_cmd1.log" chatter=1'
  111 +
  112 + # Check if execution failed
  113 + self.test_assert(self._execute(cmd, success=False) != 0,
  114 + 'Check invalid source name when executed from command line')
79 115  
80 116 # Check csscs --help
81 117 self._check_help(csscs)
... ... @@ -93,24 +129,266 @@ class Test(test):
93 129  
94 130 # Check that saving does not nothing
95 131 script['logfile'] = 'csscs_py0.log'
  132 + script['outfile'] = 'csscs_py0.fits'
96 133 script.logFileOpen()
97 134 script.save()
98 135 self.test_assert(not os.path.isfile('csscs_py0.fits'),
99   - 'Check that no FITS file has been created')
  136 + 'Check that no FITS file has been created')
100 137  
101 138 # Check that clearing does not lead to an exception or segfault
102 139 script.clear()
103 140  
104   - # Now set csscs parameters
  141 + # Set up csscs for simple binned analysis
105 142 script['logfile'] = 'csscs_py1.log'
  143 + script['inobs'] = self._cntcube
  144 + script['expcube'] = self._expcube
  145 + script['psfcube'] = self._psfcube
  146 + script['bkgcube'] = self._bkgcube
  147 + script['inmodel'] = self._bkgmodel
  148 + script['srcnames'] = 'Crab'
  149 + script['emin'] = 1.
  150 + script['emax'] = 100.
  151 + script['nxpix'] = 2
  152 + script['nypix'] = 2
  153 + script['binsz'] = 0.1
  154 + script['coordsys'] = 'CEL'
  155 + script['proj'] = 'CAR'
  156 + script['xref'] = 83.63
  157 + script['yref'] = 22.01
  158 + script['rad'] = 0.2
  159 + script['outfile'] = 'csscs_py1.fits'
  160 + script['chatter'] = 4
  161 +
  162 + # Run csscs script and save results
  163 + script.logFileOpen()
  164 + script.execute()
  165 +
  166 + # Check result
  167 + self._check_result_file('csscs_py1.fits')
  168 +
  169 + # Set up for same test without TS computation
  170 + script = cscripts.csscs()
  171 + script['logfile'] = 'csscs_py2.log'
  172 + script['inobs'] = self._cntcube
  173 + script['expcube'] = self._expcube
  174 + script['psfcube'] = self._psfcube
  175 + script['bkgcube'] = self._bkgcube
  176 + script['inmodel'] = self._bkgmodel
  177 + script['srcnames'] = 'Crab'
  178 + script['emin'] = 1.
  179 + script['emax'] = 100.
  180 + script['nxpix'] = 2
  181 + script['nypix'] = 2
  182 + script['binsz'] = 0.1
  183 + script['coordsys'] = 'CEL'
  184 + script['proj'] = 'CAR'
  185 + script['xref'] = 83.63
  186 + script['yref'] = 22.01
  187 + script['rad'] = 0.2
  188 + script['calc_ts'] = False
  189 + script['outfile'] = 'csscs_py2.fits'
  190 + script['chatter'] = 2
  191 +
  192 + # Run csscs script and save results
  193 + script.logFileOpen()
  194 + script.execute()
  195 +
  196 + # Check result
  197 + self._check_result_file('csscs_py2.fits', calc_ts=False)
  198 +
  199 + # # Set up for same test with upper limit computation
  200 + script = cscripts.csscs()
  201 + script['logfile'] = 'csscs_py3.log'
  202 + script['inobs'] = self._cntcube
  203 + script['expcube'] = self._expcube
  204 + script['psfcube'] = self._psfcube
  205 + script['bkgcube'] = self._bkgcube
  206 + script['inmodel'] = self._bkgmodel
  207 + script['srcnames'] = 'Crab'
  208 + script['emin'] = 1.
  209 + script['emax'] = 100.
  210 + script['nxpix'] = 2
  211 + script['nypix'] = 2
  212 + script['binsz'] = 0.1
  213 + script['coordsys'] = 'CEL'
  214 + script['proj'] = 'CAR'
  215 + script['xref'] = 83.63
  216 + script['yref'] = 22.01
  217 + script['rad'] = 0.2
  218 + script['calc_ts'] = True
  219 + script['calc_ulim'] = True
  220 + script['outfile'] = 'csscs_py3.fits'
  221 + script['chatter'] = 2
  222 +
  223 + # Run csscs script and save results
  224 + script.logFileOpen()
  225 + script.execute()
  226 +
  227 + # Check result
  228 + self._check_result_file('csscs_py3.fits', calc_ulim=True)
  229 +
  230 + # # Check that results can be piped in memory
  231 + self._check_results_mem(script, calc_ulim=True)
  232 +
  233 + # Set up for test without multiprocessing
  234 + script = cscripts.csscs()
  235 + script['logfile'] = 'csscs_py4.log'
  236 + script['inobs'] = self._cntcube
  237 + script['expcube'] = self._expcube
  238 + script['psfcube'] = self._psfcube
  239 + script['bkgcube'] = self._bkgcube
  240 + script['inmodel'] = self._bkgmodel
  241 + script['srcnames'] = 'Crab'
  242 + script['emin'] = 1.
  243 + script['emax'] = 100.
  244 + script['nxpix'] = 2
  245 + script['nypix'] = 2
  246 + script['binsz'] = 0.1
  247 + script['coordsys'] = 'CEL'
  248 + script['proj'] = 'CAR'
  249 + script['xref'] = 83.63
  250 + script['yref'] = 22.01
  251 + script['rad'] = 0.2
  252 + script['outfile'] = 'csscs_py4.fits'
  253 + script['nthreads'] = 1
  254 + script['chatter'] = 2
  255 +
  256 + # Run csscs script and save results
  257 + script.logFileOpen()
  258 + script.execute()
  259 +
  260 + # Check result
  261 + self._check_result_file('csscs_py4.fits')
  262 +
  263 + # Set up for test with energy dispersion
  264 + script = cscripts.csscs()
  265 + script['logfile'] = 'csscs_py5.log'
  266 + script['inobs'] = self._cntcube
  267 + script['expcube'] = self._expcube
  268 + script['psfcube'] = self._psfcube
  269 + script['bkgcube'] = self._bkgcube
  270 + script['edisp'] = True
  271 + script['edispcube'] = self._edispcube
  272 + script['inmodel'] = self._bkgmodel
  273 + script['srcnames'] = 'Crab'
  274 + script['emin'] = 1.
  275 + script['emax'] = 100.
  276 + script['nxpix'] = 2
  277 + script['nypix'] = 2
  278 + script['binsz'] = 0.1
  279 + script['coordsys'] = 'CEL'
  280 + script['proj'] = 'CAR'
  281 + script['xref'] = 83.63
  282 + script['yref'] = 22.01
  283 + script['rad'] = 0.2
  284 + script['outfile'] = 'csscs_py5.fits'
  285 + script['chatter'] = 2
  286 +
  287 + # Run csscs script and save results
  288 + script.logFileOpen()
  289 + script.execute()
  290 +
  291 + # Check result
  292 + self._check_result_file('csscs_py5.fits')
  293 +
  294 + # TODO: test unbinned analysis
  295 +
  296 + # Set up for test in On/Off mode with background model
  297 + script = cscripts.csscs()
  298 + script['logfile'] = 'csscs_py6.log'
  299 + script['inobs'] = self._evs_offaxis
  300 + script['caldb'] = self._caldb
  301 + script['irf'] = self._irf
  302 + script['inmodel'] = self._model
  303 + script['srcnames'] = 'Crab'
  304 + script['emin'] = 1.
  305 + script['emax'] = 100.
  306 + script['nxpix'] = 2
  307 + script['nypix'] = 2
  308 + script['binsz'] = 0.1
  309 + script['coordsys'] = 'CEL'
  310 + script['proj'] = 'CAR'
  311 + script['xref'] = 83.63
  312 + script['yref'] = 22.01
  313 + script['rad'] = 0.2
  314 + script['method'] = 'ONOFF'
  315 + script['inexclusion'] = self._exclusion
  316 + script['enumbins'] = 2
  317 + script['outfile'] = 'csscs_py6.fits'
  318 + script['chatter'] = 2
  319 +
  320 + # Run csscs script and save results
  321 + script.logFileOpen()
  322 + script.execute()
  323 +
  324 + # Check result
  325 + self._check_result_file('csscs_py6.fits')
  326 +
  327 + # Same test with exclusion map piped in memory
  328 + script = cscripts.csscs()
  329 + script['logfile'] = 'csscs_py7.log'
  330 + script['inobs'] = self._evs_offaxis
  331 + script['caldb'] = self._caldb
  332 + script['irf'] = self._irf
  333 + script['inmodel'] = self._model
  334 + script['srcnames'] = 'Crab'
  335 + script['emin'] = 1.
  336 + script['emax'] = 100.
  337 + script['nxpix'] = 2
  338 + script['nypix'] = 2
  339 + script['binsz'] = 0.1
  340 + script['coordsys'] = 'CEL'
  341 + script['proj'] = 'CAR'
  342 + script['xref'] = 83.63
  343 + script['yref'] = 22.01
  344 + script['rad'] = 0.2
  345 + script['method'] = 'ONOFF'
  346 + script['enumbins'] = 2
  347 + script['outfile'] = 'csscs_py7.fits'
  348 + script['chatter'] = 2
  349 +
  350 + # Set exclusion map
  351 + exclmap = gammalib.GSkyRegionMap(self._exclusion)
  352 + script.exclusion_map(exclmap)
  353 +
  354 + # Run csscs script and save results
  355 + script.logFileOpen()
  356 + script.execute()
  357 +
  358 + # Check result
  359 + self._check_result_file('csscs_py7.fits')
  360 +
  361 + # Set up for test in On/Off mode without background model
  362 + script = cscripts.csscs()
  363 + script['logfile'] = 'csscs_py8.log'
  364 + script['inobs'] = self._evs_offaxis
  365 + script['caldb'] = self._caldb
  366 + script['irf'] = self._irf
  367 + script['inmodel'] = self._model
  368 + script['srcnames'] = 'Crab'
  369 + script['emin'] = 1.
  370 + script['emax'] = 100.
  371 + script['nxpix'] = 2
  372 + script['nypix'] = 2
  373 + script['binsz'] = 0.1
  374 + script['coordsys'] = 'CEL'
  375 + script['proj'] = 'CAR'
  376 + script['xref'] = 83.63
  377 + script['yref'] = 22.01
  378 + script['rad'] = 0.2
  379 + script['method'] = 'ONOFF'
  380 + script['use_model_bkg'] = False
  381 + script['inexclusion'] = self._exclusion
  382 + script['enumbins'] = 2
  383 + script['outfile'] = 'csscs_py8.fits'
106 384 script['chatter'] = 2
107 385  
108   - # Run csscs script
  386 + # Run csscs script and save results
109 387 script.logFileOpen()
110   - script.run()
  388 + script.execute()
111 389  
112 390 # Check result
113   - self._check_result(script)
  391 + self._check_result_file('csscs_py8.fits')
114 392  
115 393 # Return
116 394 return
... ... @@ -126,32 +404,150 @@ class Test(test):
126 404 # Setup script for pickling text
127 405 script = cscripts.csscs()
128 406 script['logfile'] = 'csscs_py1_pickle.log'
  407 + script['inobs'] = self._cntcube
  408 + script['expcube'] = self._expcube
  409 + script['psfcube'] = self._psfcube
  410 + script['bkgcube'] = self._bkgcube
  411 + script['inmodel'] = self._bkgmodel
  412 + script['srcnames'] = 'Crab'
  413 + script['emin'] = 1.
  414 + script['emax'] = 100.
  415 + script['nxpix'] = 2
  416 + script['nypix'] = 2
  417 + script['binsz'] = 0.1
  418 + script['coordsys'] = 'CEL'
  419 + script['proj'] = 'CAR'
  420 + script['xref'] = 83.63
  421 + script['yref'] = 22.01
  422 + script['rad'] = 0.2
  423 + script['outfile'] = 'csscs_py1_pickle.fits'
129 424 script['chatter'] = 2
130 425  
131 426 # Perform pickeling tests of filled class
132 427 obj = self._pickeling(script)
133 428  
134 429 # Run csscs script and save result
135   - obj.logFileOpen() # Make sure we get a log file
136   - obj.run()
137   - obj.save()
  430 + obj.logFileOpen() # Make sure we get a log file
  431 + obj.execute()
138 432  
139   - # TODO: Check result file
  433 + # Check result file
  434 + self._check_result_file('csscs_py1_pickle.fits')
140 435  
141 436 # Return
142 437 return
143 438  
144   - # Check csscs result
145   - def _check_result(self, script):
  439 + # Check access to results in memory'
  440 + def _check_results_mem(self, script, srcnames=['Crab'], nx=2, ny=2, calc_ts=True,
  441 + calc_ulim=False):
146 442 """
147   - Check content of script
  443 + Check results in memory
148 444  
149 445 Parameters
150 446 ----------
151 447 script : `~cscripts.csscs`
152   - csscs instance
  448 + csscs object
  449 + srcnames : list, optional
  450 + Source names
  451 + nx : int, optional
  452 + Number of X pixels
  453 + ny : int, optional
  454 + Number of Y pixels
  455 + calc_ts : bool, optional
  456 + Check TS HDU
  457 + calc_ulim : bool, optional
  458 + Check upper limit HDU
  459 + """
  460 +
  461 + # Check that the entire fits container can be extracted
  462 + fits = script.fits()
  463 + # Dump to disk and verify
  464 + fits.saveto('csscs_py3_mem.fits', True)
  465 + self._check_result_file('csscs_py3_mem.fits', srcnames=srcnames, nx=nx, ny=ny,
  466 + calc_ts=calc_ts, calc_ulim=calc_ulim)
  467 +
  468 + # Now check direct access in memory for skymaps
  469 + # Loop over sources
  470 + for name in srcnames:
  471 + # Check flux map
  472 + skymap = script.flux(name)
  473 + self._check_map(skymap, nx=nx, ny=ny)
  474 + # Check flux error map
  475 + skymap = script.flux_error(name)
  476 + self._check_map(skymap, nx=nx, ny=ny)
  477 + # If requested check TS map
  478 + if calc_ts:
  479 + skymap = script.ts(name)
  480 + self._check_map(skymap, nx=nx, ny=ny)
  481 + # If requested check upper limit map
  482 + if calc_ulim:
  483 + skymap = script.ulimit(name)
  484 + self._check_map(skymap, nx=nx, ny=ny)
  485 +
  486 + return
  487 +
  488 + # Check result file
  489 + def _check_result_file(self, filename, srcnames=['Crab'], nx=2, ny=2,
  490 + calc_ts=True, calc_ulim=False):
  491 + """
  492 + Check result file
  493 +
  494 + Parameters
  495 + ----------
  496 + filename : str
  497 + Sky map file name
  498 + srcnames : list, optional
  499 + Source names
  500 + nx : int, optional
  501 + Number of X pixels
  502 + ny : int, optional
  503 + Number of Y pixels
  504 + calc_ts : bool, optional
  505 + Check TS HDU
  506 + calc_ulim : bool, optional
  507 + Check upper limit HDU
  508 + """
  509 +
  510 + print('ho', calc_ts)
  511 + # Loop over source names
  512 + for name in srcnames:
  513 + # Check flux map
  514 + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX]')
  515 + self._check_map(skymap, nx=nx, ny=ny)
  516 + # Check flux error map
  517 + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX ERROR]')
  518 + self._check_map(skymap, nx=nx, ny=ny)
  519 + # If requested check TS map
  520 + if calc_ts:
  521 + print('hoho')
  522 + skymap = gammalib.GSkyMap(filename + '[' + name + ' TS]')
  523 + self._check_map(skymap, nx=nx, ny=ny)
  524 + # If requested check upper limit map
  525 + if calc_ulim:
  526 + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX UPPER LIMIT]')
  527 + self._check_map(skymap, nx=nx, ny=ny)
  528 +
  529 + # Return
  530 + return
  531 +
  532 + # Check sky map
  533 + def _check_map(self, skymap, nx=2, ny=2):
153 534 """
154   - # TODO: Implement test on script result
  535 + Check sky map
  536 +
  537 + Parameters
  538 + ----------
  539 + skymap : `~gammalib.GSkyMap`
  540 + Sky map
  541 + nx : int, optional
  542 + Number of X pixels
  543 + ny : int, optional
  544 + Number of Y pixels
  545 + """
  546 +
  547 + # Check dimensions
  548 + self.test_value(skymap.nmaps(), 1, 'Check number of maps')
  549 + self.test_value(skymap.nx(), nx, 'Check for number of X pixels')
  550 + self.test_value(skymap.ny(), ny, 'Check for number of Y pixels')
155 551  
156 552 # Return
157 553 return
... ...