From 9600793a6504cb003846a40731fec91b2d33b0df Mon Sep 17 00:00:00 2001 From: Luigi Tibaldo Date: Fri, 14 Feb 2020 10:44:34 +0100 Subject: [PATCH] added test unit --- cscripts/csscs.py | 162 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------- cscripts/obsutils.py | 19 +++++++++++-------- test/test_csscs.py | 432 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------ 3 files changed, 507 insertions(+), 106 deletions(-) diff --git a/cscripts/csscs.py b/cscripts/csscs.py index 59b09fb..647521f 100755 --- a/cscripts/csscs.py +++ b/cscripts/csscs.py @@ -18,12 +18,12 @@ # along with this program. If not, see . # # ========================================================================== -import sys -import gammalib import ctools +import gammalib +import math +import sys from cscripts import mputils from cscripts import obsutils -import math # ================= # @@ -53,7 +53,8 @@ class csscs(ctools.csobservation): self._srcnames = [] self._method = None self._fits = None - self._excl_reg_map = None # Exclusion region map for on/off analysis + self._excl_reg_map = None # Exclusion region map for on/off analysis + self._excl_reg_name = None # Return return @@ -70,12 +71,13 @@ class csscs(ctools.csobservation): """ # Set pickled dictionary # Set pickled dictionary - state = {'base' : ctools.csobservation.__getstate__(self), - 'nthreads' : self._nthreads, - 'method' : self._method, - 'srcnames' : self._srcnames, - 'fits' : self._fits, - 'excl_reg_map' : self._excl_reg_map} + state = {'base': ctools.csobservation.__getstate__(self), + 'nthreads': self._nthreads, + 'method': self._method, + 'srcnames': self._srcnames, + 'fits': self._fits, + 'excl_reg_map': self._excl_reg_map, + 'excl_reg_name' : self._excl_reg_name} # Return pickled dictionary return state @@ -96,11 +98,11 @@ class csscs(ctools.csobservation): self._method = state['method'] self._fits = state['fits'] self._excl_reg_map = state['excl_reg_map'] + self._excl_reg_name = state['excl_reg_name'] # Return return - # Private methods def _get_onoff_parameters(self): """ @@ -117,8 +119,8 @@ class csscs(ctools.csobservation): # Exclusion map set and is not empty pass elif self['inexclusion'].is_valid(): - inexclusion = self['inexclusion'].filename() - self._excl_reg_map = gammalib.GSkyRegionMap(inexclusion) + self._excl_reg_name = self['inexclusion'].filename() + self._excl_reg_map = gammalib.GSkyRegionMap(self._excl_reg_name) else: msg = 'csscs in On/Off mode requires input exclusion region.' raise RuntimeError(msg) @@ -150,8 +152,8 @@ class csscs(ctools.csobservation): # Collect number of unbinned, binned and On/Off observations in # observation container n_unbinned = 0 - n_binned = 0 - n_onoff = 0 + n_binned = 0 + n_onoff = 0 for obs in self.obs(): if obs.classname() == 'GCTAObservation': if obs.eventtype() == 'CountsCube': @@ -160,16 +162,16 @@ class csscs(ctools.csobservation): n_unbinned += 1 elif obs.classname() == 'GCTAOnOffObservation': n_onoff += 1 - n_cta = n_unbinned + n_binned + n_onoff + n_cta = n_unbinned + n_binned + n_onoff n_other = self.obs().size() - n_cta # Check if there are On/Off or non-IACT observations if n_onoff > 0 or n_other > 0: - msg = 'On/Off or non-CTA observations found. '\ + msg = 'On/Off or non-CTA observations found. ' \ 'csscs does not support this type of observations.' raise RuntimeError(msg) # Otherwise check if there is a mix of binned and unbinned - elif n_unbinned > 0 and n_binned>0: + elif n_unbinned > 0 and n_binned > 0: msg = 'Mix of unbinned and binned CTA observations ' \ 'found in observation container. csscs does not ' \ 'support this mix.' @@ -228,7 +230,7 @@ class csscs(ctools.csobservation): # Loop over models and verify if model is of type GModelSky for model in self.obs().models(): if model.classname() == 'GModelSky': - nsources +=1 + nsources += 1 # If there are background gamma-ray sources in the model # throw runtime error if nsources > len(self._srcnames): @@ -327,7 +329,7 @@ class csscs(ctools.csobservation): # Free the normalisation parameter which is assumed to be # the first spectral parameter if normpar.is_fixed(): - self._log_string(gammalib.EXPLICIT, ' Freeing "'+normpar.name()+'"') + self._log_string(gammalib.EXPLICIT, ' Freeing "' + normpar.name() + '"') normpar.free() # Optionally compute Test Statistic value @@ -338,7 +340,7 @@ class csscs(ctools.csobservation): elif self['fix_bkg'].boolean() and not model.classname() == 'GModelSky': if normpar.is_free(): - self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"') + self._log_string(gammalib.EXPLICIT, ' Fixing "' + normpar.name() + '"') normpar.fix() @@ -346,13 +348,13 @@ class csscs(ctools.csobservation): elif self['fix_srcs'].boolean() and model.classname() == 'GModelSky': if normpar.is_free(): - self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"') + self._log_string(gammalib.EXPLICIT, ' Fixing "' + normpar.name() + '"') normpar.fix() # Return return - def _mask_cube(self,obs,ra,dec,rad): + def _mask_cube(self, obs, ra, dec, rad): """ Mask cube observation @@ -392,7 +394,7 @@ class csscs(ctools.csobservation): # Return return new_obs - def _mask_evlist(self,obs,ra,dec,rad): + def _mask_evlist(self, obs, ra, dec, rad): """ Mask event list @@ -434,7 +436,7 @@ class csscs(ctools.csobservation): # Return return new_obs - def _mask_observations(self,ra,dec,rad): + def _mask_observations(self, ra, dec, rad): """ Create observations restricted to circular ROI @@ -457,14 +459,14 @@ class csscs(ctools.csobservation): elif self._method == 'BINNED': new_obs = self._mask_cube(self.obs(), ra, dec, rad) elif self._method == 'ONOFF': - new_obs = obsutils.get_onoff_obs(self,self.obs(),nthreads=1, - ra = ra, dec = dec, - srcname = self._srcnames[0]) + new_obs = obsutils.get_onoff_obs(self, self.obs(), nthreads=1, + ra=ra, dec=dec, + srcname=self._srcnames[0]) # Return return new_obs - def _pixel_analysis(self,inx): + def _pixel_analysis(self, inx): """ Performs analysis over the region of interest corresponding to a single pixel of output map @@ -487,25 +489,25 @@ class csscs(ctools.csobservation): dec = pixdir.dec_deg() # Write header for spatial bin - msg = 'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx,ra,dec) + msg = 'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx, ra, dec) self._log_header2(gammalib.TERSE, msg) # Initialize results result = {} for name in self._srcnames: - result[name] = {'flux' : 0.0, + result[name] = {'flux': 0.0, 'flux_err': 0.0, - 'TS' : 0.0, - 'ulimit' : -1.0} + 'TS': 0.0, + 'ulimit': -1.0} # Mask observations self._log_header3(gammalib.NORMAL, 'Masking observations') - masked_obs = self._mask_observations(ra,dec,self['rad'].real()) + masked_obs = self._mask_observations(ra, dec, self['rad'].real()) # Set up likelihood analysis self._log_header3(gammalib.NORMAL, 'Performing fit in region') like = ctools.ctlike(masked_obs) - like['edisp'] = self['edisp'].boolean() + like['edisp'] = self['edisp'].boolean() like['nthreads'] = 1 # Avoids OpenMP conflict # If chatter level is verbose and debugging is requested then @@ -522,10 +524,10 @@ class csscs(ctools.csobservation): # Prepare objects for flux extraction # ROI centre = gammalib.GSkyDir() - centre.radec_deg(ra,dec) - roi = gammalib.GSkyRegionCircle(centre,self['rad'].real()) + centre.radec_deg(ra, dec) + roi = gammalib.GSkyRegionCircle(centre, self['rad'].real()) # Energy boundaries - emin = gammalib.GEnergy(self['emin'].real(),'TeV') + emin = gammalib.GEnergy(self['emin'].real(), 'TeV') emax = gammalib.GEnergy(self['emax'].real(), 'TeV') # Continue only if log-likelihood is non-zero @@ -539,11 +541,12 @@ class csscs(ctools.csobservation): # Get flux # Integrate spectral model between emin and emax - flux = source.spectral().flux(emin,emax) + flux = source.spectral().flux(emin, emax) # Calculate correction factor # Spatial model flux over ROI divided by ROI solid angle corr_factor = source.spatial().flux(roi) - corr_factor /= gammalib.twopi * (1 - math.cos(math.radians(self['rad'].real()))) + corr_factor /= gammalib.twopi * ( + 1 - math.cos(math.radians(self['rad'].real()))) # Multiply flux by correction factor flux *= corr_factor result[name]['flux'] = flux @@ -555,7 +558,7 @@ class csscs(ctools.csobservation): # Relative error on flux is same as on normalization # Avoid zero division error if normpar.value() > 0.: - flux_error = flux * normpar.error()/normpar.value() + flux_error = flux * normpar.error() / normpar.value() else: flux_error = 0. result[name]['flux_error'] = flux_error @@ -591,58 +594,57 @@ class csscs(ctools.csobservation): result[name]['ulimit'] = ulimit_value except: self._log_string(gammalib.NORMAL, 'Upper limit ' - 'calculation failed.') + 'calculation failed.') else: value = 'Likelihood is zero. Bin is skipped.' - self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra,dec), value) + self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra, dec), value) # Return result return result - def _write_hdu_keywords(self,hdu, is_flux = False): + def _write_hdu_keywords(self, hdu, is_flux=False): """ Append cards to HDU Parameters ---------- - hdu : `dict` HDU + hdu : `~gammalib.GFitsHDU` HDU """ # Flux units if is_flux: - hdu.card('BUNIT', 'ph/cm2/s/sr' , 'Photon flux') + hdu.card('BUNIT', 'ph/cm2/s/sr', 'Photon flux') # Minimum and maximum energy - hdu.card("E_MIN", self['emin'].real(), "[TeV] Lower energy boundary") - hdu.card("E_MAX", self['emax'].real(), "[TeV] Upper energy boundary") - hdu.card("EUNIT", "TeV", "Units for E_MIN and E_MAX") + hdu.card("E_MIN", self['emin'].real(), "[TeV] Lower energy boundary") + hdu.card("E_MAX", self['emax'].real(), "[TeV] Upper energy boundary") + hdu.card("EUNIT", "TeV", "Units for E_MIN and E_MAX") # ROI size - hdu.card('ROISZ', self['rad'].real() , + hdu.card('ROISZ', self['rad'].real(), '[deg] Size of ROI used for component separation') # Analysis method - hdu.card('METHOD',self._method, 'Analysis method') + hdu.card('METHOD', self._method, 'Analysis method') # For On/Off log information if self._method == 'ONOFF': # Exclusion map - if (self._excl_reg_map is not None) and self['inexclusion'].is_valid(): - hdu.card("EXCLMAP", self['inexclusion'].filename().url(), "Exclusion map name") + if self._excl_reg_name is not None: + hdu.card("EXCLMAP", self._excl_reg_name.url(), "Exclusion map name") # Model background? - hdu.card("BKGMODEL", self['use_model_bkg'].boolean(),'Use background model?') + hdu.card("BKGMODEL", self['use_model_bkg'].boolean(), 'Use background model?') # Fit statistic - hdu.card("STAT", self['statistic'].string(),'Fit statistic') + hdu.card("STAT", self['statistic'].string(), 'Fit statistic') return - - def _fill_fits(self,results): + def _fill_fits(self, results): """ Fill FITS object to store the results @@ -674,27 +676,29 @@ class csscs(ctools.csobservation): # Loop over pixels and fill maps for inx in range(fmap.npix()): - # Fill maps - fmap[inx,0] = results[inx][name]['flux'] + fmap[inx, 0] = results[inx][name]['flux'] errmap[inx, 0] = results[inx][name]['flux_error'] tsmap[inx, 0] = results[inx][name]['TS'] ulmap[inx, 0] = results[inx][name]['ulimit'] + # Uppercase name to adhere to gammalib convention + Name = gammalib.toupper(name) + # Write maps to Fits - fhdu = fmap.write(fits, name + ' FLUX') - self._write_hdu_keywords(fhdu,is_flux=True) - errhdu = errmap.write(fits, name + ' FLUX ERROR') + fhdu = fmap.write(fits, Name + ' FLUX') + self._write_hdu_keywords(fhdu, is_flux=True) + errhdu = errmap.write(fits, Name + ' FLUX ERROR') self._write_hdu_keywords(errhdu, is_flux=True) # If requested write TS map to fits if self['calc_ts'].boolean(): - tshdu = tsmap.write(fits, name + ' TS') + tshdu = tsmap.write(fits, Name + ' TS') self._write_hdu_keywords(tshdu) # If requested write upper limit map to fits if self['calc_ulim'].boolean(): - ulhdu = ulmap.write(fits, name + ' FLUX UPPER LIMIT') + ulhdu = ulmap.write(fits, Name + ' FLUX UPPER LIMIT') self._write_hdu_keywords(ulhdu, is_flux=True) # Return @@ -725,11 +729,13 @@ class csscs(ctools.csobservation): skymap = None # Proceed only if fits has been filled - if self._fits != None: + if self._fits != None: # Check that the name is in the list of target sources if name in self._srcnames: + # Uppercase name + Name = gammalib.toupper(name) # Assemble HDU name - hduname = name + ' ' + quantity + hduname = Name + ' ' + quantity # Try to fetch HDU try: hdu = self._fits[hduname] @@ -747,7 +753,6 @@ class csscs(ctools.csobservation): # Return return skymap - # Public methods def run(self): """ @@ -773,8 +778,8 @@ class csscs(ctools.csobservation): if self._nthreads > 1: # Compute values in pixels - args = [(self, '_pixel_analysis', i) - for i in range(npix)] + args = [(self, '_pixel_analysis', i) + for i in range(npix)] poolresults = mputils.process(self._nthreads, mputils.mpfunc, args) # Construct results @@ -804,7 +809,6 @@ class csscs(ctools.csobservation): # Continue only if FITS file is valid if self._fits != None: - # Get outmap parameter outfile = self['outfile'].filename() @@ -850,7 +854,7 @@ class csscs(ctools.csobservation): # Return return self._fits - def flux(self,name): + def flux(self, name): """ Return flux skymap @@ -865,11 +869,11 @@ class csscs(ctools.csobservation): Flux skymap """ - skymap = self._get_skymap(name,'FLUX') + skymap = self._get_skymap(name, 'FLUX') return skymap - def flux_error(self,name): + def flux_error(self, name): """ Return flux error skymap @@ -884,11 +888,11 @@ class csscs(ctools.csobservation): Flux error skymap """ - skymap = self._get_skymap(name,'FLUX ERROR') + skymap = self._get_skymap(name, 'FLUX ERROR') return skymap - def ts(self,name): + def ts(self, name): """ Return TS skymap @@ -905,7 +909,7 @@ class csscs(ctools.csobservation): # Check that TS computation is requested if self['calc_ts'].boolean(): - skymap = self._get_skymap(name,'TS') + skymap = self._get_skymap(name, 'TS') # Otherwise throw error else: msg = 'TS computation not requested. ' \ @@ -915,7 +919,7 @@ class csscs(ctools.csobservation): return skymap - def ulimit(self,name): + def ulimit(self, name): """ Return flux upper limit skymap @@ -932,7 +936,7 @@ class csscs(ctools.csobservation): # Check that upper limit computation is requested if self['calc_ulim'].boolean(): - skymap = self._get_skymap(name,'FLUX UPPER LIMIT') + skymap = self._get_skymap(name, 'FLUX UPPER LIMIT') # Otherwise throw error else: msg = 'Upper limit computation not requested. ' \ @@ -940,7 +944,6 @@ class csscs(ctools.csobservation): 'to obtain upper limit.' raise RuntimeError(msg) - return skymap @@ -948,7 +951,6 @@ class csscs(ctools.csobservation): # Main routine entry point # # ======================== # if __name__ == '__main__': - # Create instance of application app = csscs(sys.argv) diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index f30feff..0a1a8b9 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -775,18 +775,14 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = ''): # Initialise inmodel, inexclusion and use_model_bkg inmodel = 'NONE' - inexclusion = 'NONE' use_model_bkg = True - # Set inmodel, srcname, inexclusion and use_model_bkg if they are available + # Set inmodel, srcname, and use_model_bkg if they are available if cls.has_par('inmodel') and obs.models().size() == 0: if cls['inmodel'].is_valid(): inmodel = cls['inmodel'].value() if srcname == '' and cls.has_par('srcname'): srcname = cls['srcname'].value() - if cls.has_par('inexclusion'): - if cls['inexclusion'].is_valid(): - inexclusion = cls['inexclusion'].value() if cls.has_par('use_model_bkg'): use_model_bkg = cls['use_model_bkg'].boolean() @@ -794,7 +790,6 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = ''): phagen = cscripts.csphagen(obs) phagen['inmodel'] = inmodel phagen['srcname'] = srcname - phagen['inexclusion'] = inexclusion phagen['emin'] = cls['emin'].real() phagen['emax'] = cls['emax'].real() phagen['enumbins'] = cls['enumbins'].integer() @@ -830,13 +825,21 @@ def get_onoff_obs(cls, obs, nthreads=0, ra = None, dec = None, srcname = ''): phagen['debug'] = cls['debug'].boolean() phagen['nthreads'] = nthreads - # Pipe exclusion map + # Pipe exclusion map in memory if hasattr(cls, 'exclusion_map'): exclusion_map = cls.exclusion_map() if exclusion_map is not None: phagen.exclusion_map(exclusion_map) + # Otherwise try to read it as user parameter + elif cls.has_par('inexclusion'): + if cls['inexclusion'].is_valid(): + inexclusion = cls['inexclusion'].value() + phagen['inexclusion'] = inexclusion + # Otherwise set to None + else: + phagen['inexclusion'] = 'NONE' - # Run csphagen + # Run csphagen phagen.run() # Clone resulting observation container diff --git a/test/test_csscs.py b/test/test_csscs.py index fb47c42..52bde60 100644 --- a/test/test_csscs.py +++ b/test/test_csscs.py @@ -18,9 +18,9 @@ # along with this program. If not, see . # # ========================================================================== -import os -import gammalib import cscripts +import gammalib +import os from testing import test @@ -43,6 +43,15 @@ class Test(test): # Call base class constructor test.__init__(self) + # Set test datasets and parameters + self._evs_offaxis = self._datadir + '/crab_offaxis1.fits' + self._exclusion = self._datadir + '/crab_exclusion.fits' + self._expcube = self._datadir + '/crab_expcube.fits' + self._psfcube = self._datadir + '/crab_psfcube.fits' + self._bkgcube = self._datadir + '/crab_bkgcube.fits' + self._bkgmodel = self._datadir + '/crab_bkgcube.xml' + self._edispcube = self._datadir + '/crab_edispcube.fits' + # Return return @@ -71,11 +80,38 @@ class Test(test): csscs = self._script('csscs') # Setup csscs command - cmd = csscs+' logfile="csscs_cmd1.log" chatter=1' + cmd = csscs + ' inobs="' + self._cntcube + '"' + \ + ' expcube="' + self._expcube + '"' + \ + ' psfcube="' + self._psfcube + '"' + \ + ' bkgcube="' + self._bkgcube + '"' + \ + ' inmodel="' + self._bkgmodel + '"' + \ + ' srcnames="Crab"' + \ + ' emin=1. emax=100.0 nxpix=2 nypix=2' + \ + ' binsz=0.1 coordsys="CEL" proj="CAR"' + \ + ' xref=83.63 yref=22.01 rad=0.2' + \ + ' outfile="csscs_cmd1.fits"' + \ + ' logfile="csscs_cmd1.log" chatter=1' # Check if execution was successful self.test_value(self._execute(cmd), 0, - 'Check successful execution from command line') + 'Check successful execution from command line') + + # Setup csscs command + cmd = csscs + ' inobs="' + self._cntcube + '"' + \ + ' expcube="' + self._expcube + '"' + \ + ' psfcube="' + self._psfcube + '"' + \ + ' bkgcube="' + self._bkgcube + '"' + \ + ' inmodel="' + self._bkgmodel + '"' + \ + ' srcnames="Crab;source_that_is_not_in_the_model"' + \ + ' emin=1. emax=100.0 nxpix=2 nypix=2' + \ + ' binsz=0.1 coordsys="CEL" proj="CAR"' + \ + ' xref=83.63 yref=22.01 rad=0.2' + \ + ' outfile="csscs_cmd1.fits"' + \ + ' logfile="csscs_cmd1.log" chatter=1' + + # Check if execution failed + self.test_assert(self._execute(cmd, success=False) != 0, + 'Check invalid source name when executed from command line') # Check csscs --help self._check_help(csscs) @@ -93,24 +129,266 @@ class Test(test): # Check that saving does not nothing script['logfile'] = 'csscs_py0.log' + script['outfile'] = 'csscs_py0.fits' script.logFileOpen() script.save() self.test_assert(not os.path.isfile('csscs_py0.fits'), - 'Check that no FITS file has been created') + 'Check that no FITS file has been created') # Check that clearing does not lead to an exception or segfault script.clear() - # Now set csscs parameters + # Set up csscs for simple binned analysis script['logfile'] = 'csscs_py1.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['outfile'] = 'csscs_py1.fits' + script['chatter'] = 4 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py1.fits') + + # Set up for same test without TS computation + script = cscripts.csscs() + script['logfile'] = 'csscs_py2.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['calc_ts'] = False + script['outfile'] = 'csscs_py2.fits' + script['chatter'] = 2 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py2.fits', calc_ts=False) + + # # Set up for same test with upper limit computation + script = cscripts.csscs() + script['logfile'] = 'csscs_py3.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['calc_ts'] = True + script['calc_ulim'] = True + script['outfile'] = 'csscs_py3.fits' + script['chatter'] = 2 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py3.fits', calc_ulim=True) + + # # Check that results can be piped in memory + self._check_results_mem(script, calc_ulim=True) + + # Set up for test without multiprocessing + script = cscripts.csscs() + script['logfile'] = 'csscs_py4.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['outfile'] = 'csscs_py4.fits' + script['nthreads'] = 1 + script['chatter'] = 2 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py4.fits') + + # Set up for test with energy dispersion + script = cscripts.csscs() + script['logfile'] = 'csscs_py5.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['edisp'] = True + script['edispcube'] = self._edispcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['outfile'] = 'csscs_py5.fits' + script['chatter'] = 2 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py5.fits') + + # TODO: test unbinned analysis + + # Set up for test in On/Off mode with background model + script = cscripts.csscs() + script['logfile'] = 'csscs_py6.log' + script['inobs'] = self._evs_offaxis + script['caldb'] = self._caldb + script['irf'] = self._irf + script['inmodel'] = self._model + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['method'] = 'ONOFF' + script['inexclusion'] = self._exclusion + script['enumbins'] = 2 + script['outfile'] = 'csscs_py6.fits' + script['chatter'] = 2 + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py6.fits') + + # Same test with exclusion map piped in memory + script = cscripts.csscs() + script['logfile'] = 'csscs_py7.log' + script['inobs'] = self._evs_offaxis + script['caldb'] = self._caldb + script['irf'] = self._irf + script['inmodel'] = self._model + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['method'] = 'ONOFF' + script['enumbins'] = 2 + script['outfile'] = 'csscs_py7.fits' + script['chatter'] = 2 + + # Set exclusion map + exclmap = gammalib.GSkyRegionMap(self._exclusion) + script.exclusion_map(exclmap) + + # Run csscs script and save results + script.logFileOpen() + script.execute() + + # Check result + self._check_result_file('csscs_py7.fits') + + # Set up for test in On/Off mode without background model + script = cscripts.csscs() + script['logfile'] = 'csscs_py8.log' + script['inobs'] = self._evs_offaxis + script['caldb'] = self._caldb + script['irf'] = self._irf + script['inmodel'] = self._model + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['method'] = 'ONOFF' + script['use_model_bkg'] = False + script['inexclusion'] = self._exclusion + script['enumbins'] = 2 + script['outfile'] = 'csscs_py8.fits' script['chatter'] = 2 - # Run csscs script + # Run csscs script and save results script.logFileOpen() - script.run() + script.execute() # Check result - self._check_result(script) + self._check_result_file('csscs_py8.fits') # Return return @@ -126,32 +404,150 @@ class Test(test): # Setup script for pickling text script = cscripts.csscs() script['logfile'] = 'csscs_py1_pickle.log' + script['inobs'] = self._cntcube + script['expcube'] = self._expcube + script['psfcube'] = self._psfcube + script['bkgcube'] = self._bkgcube + script['inmodel'] = self._bkgmodel + script['srcnames'] = 'Crab' + script['emin'] = 1. + script['emax'] = 100. + script['nxpix'] = 2 + script['nypix'] = 2 + script['binsz'] = 0.1 + script['coordsys'] = 'CEL' + script['proj'] = 'CAR' + script['xref'] = 83.63 + script['yref'] = 22.01 + script['rad'] = 0.2 + script['outfile'] = 'csscs_py1_pickle.fits' script['chatter'] = 2 # Perform pickeling tests of filled class obj = self._pickeling(script) # Run csscs script and save result - obj.logFileOpen() # Make sure we get a log file - obj.run() - obj.save() + obj.logFileOpen() # Make sure we get a log file + obj.execute() - # TODO: Check result file + # Check result file + self._check_result_file('csscs_py1_pickle.fits') # Return return - # Check csscs result - def _check_result(self, script): + # Check access to results in memory' + def _check_results_mem(self, script, srcnames=['Crab'], nx=2, ny=2, calc_ts=True, + calc_ulim=False): """ - Check content of script + Check results in memory Parameters ---------- script : `~cscripts.csscs` - csscs instance + csscs object + srcnames : list, optional + Source names + nx : int, optional + Number of X pixels + ny : int, optional + Number of Y pixels + calc_ts : bool, optional + Check TS HDU + calc_ulim : bool, optional + Check upper limit HDU + """ + + # Check that the entire fits container can be extracted + fits = script.fits() + # Dump to disk and verify + fits.saveto('csscs_py3_mem.fits', True) + self._check_result_file('csscs_py3_mem.fits', srcnames=srcnames, nx=nx, ny=ny, + calc_ts=calc_ts, calc_ulim=calc_ulim) + + # Now check direct access in memory for skymaps + # Loop over sources + for name in srcnames: + # Check flux map + skymap = script.flux(name) + self._check_map(skymap, nx=nx, ny=ny) + # Check flux error map + skymap = script.flux_error(name) + self._check_map(skymap, nx=nx, ny=ny) + # If requested check TS map + if calc_ts: + skymap = script.ts(name) + self._check_map(skymap, nx=nx, ny=ny) + # If requested check upper limit map + if calc_ulim: + skymap = script.ulimit(name) + self._check_map(skymap, nx=nx, ny=ny) + + return + + # Check result file + def _check_result_file(self, filename, srcnames=['Crab'], nx=2, ny=2, + calc_ts=True, calc_ulim=False): + """ + Check result file + + Parameters + ---------- + filename : str + Sky map file name + srcnames : list, optional + Source names + nx : int, optional + Number of X pixels + ny : int, optional + Number of Y pixels + calc_ts : bool, optional + Check TS HDU + calc_ulim : bool, optional + Check upper limit HDU + """ + + print('ho', calc_ts) + # Loop over source names + for name in srcnames: + # Check flux map + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX]') + self._check_map(skymap, nx=nx, ny=ny) + # Check flux error map + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX ERROR]') + self._check_map(skymap, nx=nx, ny=ny) + # If requested check TS map + if calc_ts: + print('hoho') + skymap = gammalib.GSkyMap(filename + '[' + name + ' TS]') + self._check_map(skymap, nx=nx, ny=ny) + # If requested check upper limit map + if calc_ulim: + skymap = gammalib.GSkyMap(filename + '[' + name + ' FLUX UPPER LIMIT]') + self._check_map(skymap, nx=nx, ny=ny) + + # Return + return + + # Check sky map + def _check_map(self, skymap, nx=2, ny=2): """ - # TODO: Implement test on script result + Check sky map + + Parameters + ---------- + skymap : `~gammalib.GSkyMap` + Sky map + nx : int, optional + Number of X pixels + ny : int, optional + Number of Y pixels + """ + + # Check dimensions + self.test_value(skymap.nmaps(), 1, 'Check number of maps') + self.test_value(skymap.nx(), nx, 'Check for number of X pixels') + self.test_value(skymap.ny(), ny, 'Check for number of Y pixels') # Return return -- libgit2 0.21.4