diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 1173aae..5a6c64c 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -152,7 +152,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/models_spectral.xml \ $(top_srcdir)/test/data/model_filefunction.txt \ $(top_srcdir)/test/data/model_map.fits.gz \ - $(top_srcdir)/test/data/model_mapcube.fits.gz + $(top_srcdir)/test/data/model_mapcube.fits.gz \ + $(top_srcdir)/test/data/iact_runlist.dat # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \ diff --git a/cscripts/csiactcopy.py b/cscripts/csiactcopy.py index 8abf1ce..48d891d 100755 --- a/cscripts/csiactcopy.py +++ b/cscripts/csiactcopy.py @@ -91,7 +91,18 @@ class csiactcopy(ctools.cscript): def _copy(self, source, clobber): """ Copy file to outpath + + Parameters + ---------- + source : string + Path of file to be copied + + Returns + ------- + filesize : float + Size of the file that was copied """ + # Get file destination destination = os.path.join(self._outpath, os.path.relpath(source, self._remote_base)) @@ -147,10 +158,21 @@ class csiactcopy(ctools.cscript): def _merge(self, localfits, remotefits, hduname, clobber): """ Merge remote and local fits files - If the local fits file is not present, a new one is created. + + Parameters + ---------- + localfits : string + Path of local index FITS file + remotefits : string + Path of remote index FITS file + hduname : string + Name of HDU extension to be merged + clobber : bool + Flag if remote content should overwrite local content """ - # ... + + # Logging if self._logTerse(): self._log_value('Local file', localfits) self._log_value('Remote file', remotefits) @@ -214,7 +236,7 @@ class csiactcopy(ctools.cscript): tmp_hdu = remote[hduname] # Get old size of remote HDU - size = remote_hdu.size() + size = remote_hdu.nrows() # Extend remote HDU remote_hdu.append_rows(tmp_hdu.nrows()) @@ -229,8 +251,8 @@ class csiactcopy(ctools.cscript): robs_col = remote_hdu['OBS_ID'] for obs_id in robs_col: remote_obs.append(obs_id) - - # Dump the number of remote and local entries + + # Log entries of index files if self._logTerse(): self._log_value('Remote entries', len(remote_obs)) self._log_value('Local entries', len(local_obs)) @@ -250,12 +272,18 @@ class csiactcopy(ctools.cscript): # Remove local obs_id entries table_has_obsid = True + + # Loop over obs ids and overwrite with remote content while table_has_obsid: + + # First remove old rows for i in range(local_hdu.nrows()): if remote_obs_id == local_hdu['OBS_ID'][i]: local_hdu.remove_rows(i,1) removed_rows += 1 break + + # Replace with remote content table_has_obsid = False for i in range(local_hdu.nrows()): if remote_obs_id == local_hdu['OBS_ID'][i]: @@ -290,8 +318,8 @@ class csiactcopy(ctools.cscript): # Add remote HDUs local_hdu.insert_rows(old_local_rows, remote_hdu.nrows()) - - # Dump ... + + # Log actions if self._logTerse(): if clobber: what = 'Removed local rows' @@ -344,6 +372,8 @@ class csiactcopy(ctools.cscript): # Check if runlist file is available if self._runlist.exists(): runfile = open(self._runlist.url()) + + # Read runlist for line in runfile.readlines(): if len(line) == 0: continue @@ -351,6 +381,8 @@ class csiactcopy(ctools.cscript): continue if len(line.split()) > 0: self._runs.append(int(line.split()[0])) + + # Close runlist file runfile.close() # Logging @@ -428,8 +460,8 @@ class csiactcopy(ctools.cscript): # Initialise flag if SIZE column is present has_size = table.contains('SIZE') + # Initialise file size if has_size: - # Initialise file size cp_size = 0 # Initialise remote observation IDs @@ -439,6 +471,8 @@ class csiactcopy(ctools.cscript): # Log runs that are not available remotely for run in self._runs: + + # Check for run not in remote data store if not run in remote_ids: if self._logNormal(): self._log('Skip observation '+str(run)+': ') @@ -529,16 +563,22 @@ class csiactcopy(ctools.cscript): self._log('\n') self._log.header1('Copying files') - # intialise counter + # Intialise counter k = 0 - # initialise values for logging + # Initialise values for logging last_fraction = 0.0 fraction_increment = 20.0 + + # Use 10% step increase if self._logNormal(): fraction_increment = 10.0 + + # Use 5% step increase if self._logTerse(): fraction_increment = 5.0 + + # Use 2% step increase if self._logExplicit(): fraction_increment = 2.0 @@ -552,6 +592,8 @@ class csiactcopy(ctools.cscript): # Log progress fraction = float(k) / float(len(files)) * 100.0 while fraction > last_fraction: + + # Print status of copying procedure if self._logNormal() and not self._logVerbose(): self._log_value('Status', '%d %%' % int(last_fraction)) last_fraction += fraction_increment @@ -565,7 +607,7 @@ class csiactcopy(ctools.cscript): total_size += filesize n_copied += 1 - # Increment counter + # Increment file counter k += 1 # Logging @@ -573,14 +615,17 @@ class csiactcopy(ctools.cscript): self._log.parformat('Status') self._log('Finished') self._log('\n') + + # Logging about index files if self._logTerse(): self._log('\n') self._log.header1('Updating index files') - # Build local index file names + # Build local hdu index file name local_hdu = os.path.join(self._outpath, os.path.relpath(remote_hdu, self._remote_base)) + # Build local obs index file name local_obs = os.path.join(self._outpath, os.path.relpath(remote_obs, self._remote_base)) @@ -596,6 +641,7 @@ class csiactcopy(ctools.cscript): # Merge remote index files with local files self._merge(local_hdu, remote_hdu, 'HDU_INDEX', self._clobber()) + # Logging if self._logTerse(): self._log('\n') self._log.header3('OBS index') @@ -607,7 +653,8 @@ class csiactcopy(ctools.cscript): # If all files were copied, just copy index files too self._copy(remote_hdu, self._clobber()) self._copy(remote_obs, self._clobber()) - + + # Logging if self._logTerse(): self._log('\n') self._log.header3('Master index file') diff --git a/cscripts/csiactdata.py b/cscripts/csiactdata.py index 61c6fa3..1c54d56 100755 --- a/cscripts/csiactdata.py +++ b/cscripts/csiactdata.py @@ -95,6 +95,12 @@ class csiactdata(ctools.cscript): def _log_configuration(self, config): """ Log configuration + + Parameters + ---------- + config : dict + Dictionary containing data configuration + """ # Log only if the level is at least "Terse" if self._logTerse(): @@ -107,6 +113,8 @@ class csiactdata(ctools.cscript): # Print additional information for key in config: + + # Skip keys that were treated above if key in ['name','hduindx','obsindx']: continue self._log_value(key, config[key]) diff --git a/cscripts/csiactobs.py b/cscripts/csiactobs.py index 278548d..37d9418 100755 --- a/cscripts/csiactobs.py +++ b/cscripts/csiactobs.py @@ -47,9 +47,10 @@ class csiactobs(ctools.cscript): """ # Set name and version self._name = 'csiactobs' - self._version = '1.1.0' + self._version = '1.2.0' # Initialise some members + self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() self._datapath = os.getenv('VHEFITS','') self._inmodels = gammalib.GModels() @@ -59,8 +60,8 @@ class csiactobs(ctools.cscript): self._runlist = [] self._runlistfile = gammalib.GFilename() self._bkgpars = 0 - self._outmodel = gammalib.GFilename() - self._outobs = gammalib.GFilename() + #self._outmodel = gammalib.GFilename() + #self._outobs = gammalib.GFilename() self._master_indx = '' self._use_bkg_scale = False self._ev_hiera = [''] @@ -80,8 +81,7 @@ class csiactobs(ctools.cscript): self._subdir = '' self._debug = False - # Initialise application by calling the appropriate class - # constructor. + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Return @@ -91,7 +91,7 @@ class csiactobs(ctools.cscript): # Private methods def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile and setup the observation """ # Set data path if self._datapath == '': @@ -105,7 +105,7 @@ class csiactobs(ctools.cscript): self._datapath = gammalib.expand_env(self._datapath) # Read FITS production - self._prodname = self['prodname'].string() + self._prodname = self['prodname'].string() # Read runlist file if list not already filled if len(self._runlist) == 0: @@ -125,13 +125,12 @@ class csiactobs(ctools.cscript): runfile.close() # Read number of background parameters - self._bkgpars = self['bkgpars'].integer() - - # Output model file - self._outmodel = self['outmodel'].filename() - - # Observation outfile - self._outobs = self['outobs'].filename() + self._bkgpars = self['bkgpars'].integer() + + # Query ahead output parameters + if self._read_ahead(): + self['outmodel'].filename() + self['outobs'].filename() # Master index file name self._master_indx = self['master_indx'].string() @@ -174,27 +173,45 @@ class csiactobs(ctools.cscript): # Get HDUs for config in configs: + + # Check if prodname is present if self._prodname == config['name']: self._hdu_index = str(os.path.join(self._datapath, config['hduindx'])) self._obs_index = str(os.path.join(self._datapath, config['obsindx'])) + + # Leave loop if index file names were found break - # Check HDUs + # Check index files if self._hdu_index == '' or self._obs_index == '': raise RuntimeError('*** ERROR: FITS data store "'+self._prodname+'" not available. Run csiactdata to get a list of available storage names') + + # Check HDU names filename = gammalib.GFilename(self._hdu_index+'[HDU_INDEX]') if not filename.is_fits(): raise RuntimeError('*** ERROR: HDU index file "'+self._hdu_index+'[HDU_INDEX]" for FITS data store "'+self._prodname+'" not available. Check your master index file or run csiactdata to get a list of available storage names.') # Check for existence of 'BKG_SCALE' in the observation index file if required if self._use_bkg_scale: + + # Create filename filename = gammalib.GFilename(self._obs_index+'[OBS_INDEX]') + + # Check if it is a FITS file if filename.is_fits(): + + # Open FITS file fits = gammalib.GFits(self._obs_index) + + # Check if column "BKG_SCALE" is found and signal its possible usage if not fits['OBS_INDEX'].contains('BKG_SCALE'): self._use_bkg_scale = False + + # Close FITS file fits.close() + else: + # Signal that there is no background scale self._use_bkg_scale = False # Create base data directory from hdu index file location @@ -209,33 +226,37 @@ class csiactobs(ctools.cscript): # Return return - def runlist(self, runlist): - """ - Set observation list + def _background_spectrum(self, prefactor, index, emin = 0.01, emax = 100.0): """ + Create a background spectrum model dependent on user parameters - # Check if convertable to a python list - try: - - # Store a copy of input argument - self._runlist = list(runlist) - - except: - - # raise error if wrong type detected - raise RuntimeError('Argument of csiactobs.runlist() must be of type "list"') + Parameters + ---------- + prefactor : float + Prefactor parameter of spectral model + index : float + Index parameter of spectral model + emin : float + Minimum energy (in case a spectral node function is required) + emax : float + Maximum energy (in case a spectral node function is required) - # Return - return - - def _background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0): + Returns + ------- + spec : `~gammalib.GModelSpectral()` + Spectral model for the background shape + """ # Handle constant spectral model if index == 0.0 and self._bkgpars <= 1: spec = gammalib.GModelSpectralConst() + + # Set parameter range spec['Normalization'].min(prefactor / self._bkg_range_factor) spec['Normalization'].max(prefactor * self._bkg_range_factor) spec['Normalization'].value(prefactor) + + # Freeze or release normalisation parameter if self._bkgpars == 0: spec['Normalization'].fix() else: @@ -245,8 +266,10 @@ class csiactobs(ctools.cscript): # Create power law model if self._bkgpars <= 2: - e = gammalib.GEnergy(1.0,'TeV') - spec = gammalib.GModelSpectralPlaw(prefactor, index, e) + + # Set Power Law model with arbitrary pivot energy + pivot = gammalib.GEnergy(1.0,'TeV') + spec = gammalib.GModelSpectralPlaw(prefactor, index, pivot) # Set parameter ranges spec[0].min(prefactor / self._bkg_range_factor) @@ -269,19 +292,31 @@ class csiactobs(ctools.cscript): else: # Create reference powerlaw - plaw = gammalib.GModelSpectralPlaw(prefactor, index, gammalib.GEnergy(1.0,'TeV')) + pivot = gammalib.GEnergy(1.0,'TeV') + plaw = gammalib.GModelSpectralPlaw(prefactor, index, pivot) # Create spectral model and energy values spec = gammalib.GModelSpectralNodes() + + # Create logarithmic energy nodes bounds = gammalib.GEbounds(self._bkgpars,gammalib.GEnergy(emin,'TeV'),gammalib.GEnergy(emax,'TeV'), True) + + # Loop over bounds and set intensity value for i in range(bounds.size()): energy = bounds.elogmean(i) - value = plaw.eval(energy, gammalib.GTime()) + value = plaw.eval(energy) + + # Append energy, value - tuple to node function spec.append(energy, value) + + # Loop over parameters for par in spec: + # Fix energy nodes if 'Energy' in par.name(): par.fix() + + # Release intensity nodes elif 'Intensity' in par.name(): value = par.value() par.scale(value) @@ -292,13 +327,36 @@ class csiactobs(ctools.cscript): return spec def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100): + """ + Create an IACT background model + + Parameters + ---------- + telescope : string + Name of telescope + obs_id : string + Observation ID + bkg_scale : float + Background scaling factor + bkgtype : string + Type of background (irf,aeff, or gauss) + Returns + ------- + model : `~gammalib.GModelData()` + Background model for IACT observation + """ + # handle IrfBackground if bkgtype == 'irf': + + # Set parameters to have a constant spectral model prefactor = 1.0 index = 0.0 prefactor *= bkg_scale - spec = self._background_spectrum(obs_id, prefactor, index, emin, emax) + + # Create background spectrum + spec = self._background_spectrum(prefactor, index, emin, emax) # Create background model instance bck = gammalib.GCTAModelIrfBackground(spec) @@ -306,7 +364,7 @@ class csiactobs(ctools.cscript): # Set AeffBackground elif bkgtype == 'aeff': prefactor = bkg_scale * self._bkg_aeff_norm - spec = self._background_spectrum(obs_id, prefactor, self._bkg_aeff_index, emin, emax) + spec = self._background_spectrum(prefactor, self._bkg_aeff_index, emin, emax) # Create background model instance bck = gammalib.GCTAModelAeffBackground(spec) @@ -314,12 +372,14 @@ class csiactobs(ctools.cscript): # Set Gaussian Background elif bkgtype == 'gauss': prefactor = bkg_scale * self._bkg_gauss_norm - spec = self._background_spectrum(obs_id, prefactor, self._bkg_gauss_index, emin, emax) + + # Create background spectrum + spec = self._background_spectrum(prefactor, self._bkg_gauss_index, emin, emax) radial = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma) bck = gammalib.GCTAModelRadialAcceptance(radial, spec) else: - msg = 'Background type "'+self._bkgtype+'" unsupported' + msg = 'Background type "'+bkgtype+'" unsupported' raise RuntimeError(msg) # Copy model @@ -383,6 +443,7 @@ class csiactobs(ctools.cscript): eventfile = aefffile = psffile = edispfile = bkgfile = '' eventhdu = aeffhdu = '' + # Initialise types and formats types = [] formats = [] for i in range(hduindx_hdu.nrows()): @@ -402,7 +463,11 @@ class csiactobs(ctools.cscript): hduindx_hdu['FILE_NAME'][index]) eventhdu = hduindx_hdu['HDU_NAME'][index] eventfile += '['+eventhdu+']' + + # Check if psf file is available if not gammalib.GFilename(eventfile).is_fits(): + + # Print warning that observation will be skipped since no events are present if self._logTerse(): self._log('Skipping observation '+str(obs_id)+ ': eventfile "'+eventfile+'" not found\n') @@ -421,7 +486,11 @@ class csiactobs(ctools.cscript): hduindx_hdu['FILE_NAME'][index]) aeffhdu = hduindx_hdu['HDU_NAME'][index] aefffile += '['+aeffhdu+']' + + # Check if psf file is available if not gammalib.GFilename(aefffile).is_fits(): + + # Print warning that observation will be skipped since no effective area is present if self._logTerse(): self._log('Skipping observation '+str(obs_id)+ ': effective area "'+aefffile+'" not found\n') @@ -439,7 +508,11 @@ class csiactobs(ctools.cscript): hduindx_hdu['FILE_DIR'][index], hduindx_hdu['FILE_NAME'][index]) psffile += '['+hduindx_hdu['HDU_NAME'][index]+']' + + # Check if psf file is available if not gammalib.GFilename(psffile).is_fits(): + + # Print warning that observation will be skipped since no PSF is present if self._logTerse(): self._log('Skipping observation '+str(obs_id)+ ': point spread function "'+psffile+ @@ -458,7 +531,11 @@ class csiactobs(ctools.cscript): hduindx_hdu['FILE_DIR'][index], hduindx_hdu['FILE_NAME'][index]) edispfile += '['+hduindx_hdu['HDU_NAME'][index]+']' + + # Check if edisp file is available if not gammalib.GFilename(edispfile).is_fits(): + + # Print warning that edisp cannot be used for this observation if self._logTerse(): self._log('Warning: observation '+str(obs_id)+ ' has no energy dispersion "'+edispfile+ @@ -480,17 +557,22 @@ class csiactobs(ctools.cscript): bkgfile += '['+hduindx_hdu['HDU_NAME'][index]+']' if 'irf' in bkg_mod_hierarchy and not gammalib.GFilename(bkgfile).is_fits(): bkg_mod_hierarchy.remove('irf') + + # Throw warning if IRF background cannot be used if self._logTerse(): self._log('Warning: observation '+str(obs_id)+ ' has no background information (file="'+ bkgfile+'"). IRF background cannot be used\n') bkgfile = '' if len(bkg_mod_hierarchy) == 0: + + # Skip observation if no background can be used if self._logTerse(): self._log('Skipping observation '+str(obs_id)+ ': No background can be used\n') continue else: + # Log if we fall back to next background approach if self._logTerse(): self._log('Observation '+str(obs_id)+ ': Falling back to background "'+ @@ -499,15 +581,21 @@ class csiactobs(ctools.cscript): # Close hdu index file hduindx.close() - # Handle background scale information if available + # Initialise background scale bkg_scale = 1.0 + + # Check if background scale should be used if self._use_bkg_scale: + + # Read background scale from fits file obsindx = gammalib.GFits(self._obs_index+'[OBS_INDEX]'+obs_selection) bkg_scale = obsindx['OBS_INDEX']['BKG_SCALE'][0] obsindx.close() - # Open fits file to determine the observation name + # Open event FITS file fits = gammalib.GFits(eventfile) + + # Get object and telescope strings from event hdu events = fits[eventhdu] object_name = events.string('OBJECT') telescope = events.string('TELESCOP') @@ -550,6 +638,7 @@ class csiactobs(ctools.cscript): self._log('Adding observation '+str(obs_id)+ ' ("'+object_name+'")\n') + # log more details if self._logExplicit(): self._log(' Event file: '+eventfile+'\n') self._log(' Effective area: '+aefffile+'\n') @@ -615,10 +704,13 @@ class csiactobs(ctools.cscript): self._log('\n') else: + # Print warning if no observation could be be used from runlist self._log.header2('WARNING: No observation from given runlist available') # Append models provided by 'inmodels' if necessary - if not self._inmodels == None: + if self._inmodels.size() > 0: + + # Logging if self._logTerse(): self._log('\n') self._log.header1('Appending models') @@ -650,17 +742,64 @@ class csiactobs(ctools.cscript): def save(self): """ - Save model definition XML file + Save observation definition and model definition XML file """ - # Save observation XML file - self._xml.save(self._outobs) + # Write header + if self._logTerse(): + self._log('\n') + self._log.header1('Save XML files') + + # Get filenames + outobs = self['outobs'].filename() + outmodel = self['outmodel'].filename() + + # Log file name + if self._logTerse(): + self._log(gammalib.parformat('Obs. definition XML file')) + self._log(outobs.url()) + self._log('\n') + + # Save observation definition XML file + self._xml.save(outobs.url()) + + # Log file name + if self._logTerse(): + self._log(gammalib.parformat('Model definiton XML file')) + self._log(outmodel.url()) + self._log('\n') # Save model XML file - self._models.save(self._outmodel) + self._models.save(outmodel.url()) # Return return + def runlist(self, runlist): + """ + Set runlist + + Parameters + ---------- + runlist : list + List of observation IDs + + Raises + ------ + RuntimeError + Input runlist is not a Python list + """ + # If the runlist is convertable to a Python list then store the + # runlist in the class + try: + self._runlist = list(runlist) + + # ... otherwise raise an exception + except: + raise RuntimeError('Argument is not a Python list') + + # Return + return + def ebounds(self): """ Return runlist energy boundaries @@ -670,20 +809,20 @@ class csiactobs(ctools.cscript): def obs(self): """ - Returns GObservations object + Return observations container """ - # Initialise observations - obs = gammalib.GObservations() + # Clear observations + self._obs.clear() # Read observations if XML is filled if self._xml.size(): - obs.read(self._xml) + self._obs.read(self._xml) # Assign models - obs.models(self._models) + self._obs.models(self._models) # Return observations - return obs + return self._obs # ======================== # diff --git a/pyext/ctool.i b/pyext/ctool.i index 5aec3a7..123c69b 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -198,7 +198,7 @@ cscript._init_cscript = _init_cscript # Log the value of a parameter def _log_value(self, name, value): - self._log.parformat(name) + self._log.parformat(str(name)) self._log(str(value)) self._log('\n') cscript._log_value = _log_value diff --git a/test/data/iact_runlist.dat b/test/data/iact_runlist.dat new file mode 100644 index 0000000..ff95acd --- /dev/null +++ b/test/data/iact_runlist.dat @@ -0,0 +1 @@ +0 diff --git a/test/data/iactdata/master.json b/test/data/iactdata/master.json index 51bb14c..77df9ec 100644 --- a/test/data/iactdata/master.json +++ b/test/data/iactdata/master.json @@ -3,7 +3,14 @@ { "name": "unit-test", "obsindx": "obs-index.fits.gz", - "hduindx": "hdu-index.fits.gz" + "hduindx": "hdu-index.fits.gz", + "comment": "this is a comment" + }, + { + "name": "unit-test-that-doesnt-exist", + "obsindx": "badfile.fits.gz", + "hduindx": "unavailable_file.fits.gz", + "comment": "this is a comment" } ] } diff --git a/test/test_csiactdata.py b/test/test_csiactdata.py index 47f1db6..0333ef5 100755 --- a/test/test_csiactdata.py +++ b/test/test_csiactdata.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the csiactdata script. +# This scripts performs unit tests for the csiactdata script # # Copyright (C) 2016 Juergen Knoedlseder # @@ -108,6 +108,13 @@ class Test(test): # Run csiactdata script and save run list iactdata.logFileOpen() # Make sure we get a log file iactdata.run() + + # Get available configs + names = iactdata.names() + + # Check if there is at least one config available + self.test_assert(len(names) > 0, + 'Check successful execution from Python') # Return return diff --git a/test/test_csiactobs.py b/test/test_csiactobs.py index 0c961f9..08cf4af 100755 --- a/test/test_csiactobs.py +++ b/test/test_csiactobs.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the csiactobs script. +# This scripts performs unit tests for the csiactobs script # # Copyright (C) 2016 Juergen Knoedlseder # @@ -44,6 +44,7 @@ class Test(test): # Set data members self._datapath = self._datadir + '/iactdata' + self._runlist = self._datadir + '/iact_runlist.dat' # Return return @@ -74,10 +75,10 @@ class Test(test): # Setup csiactobs command cmd = csiactobs+' datapath="'+self._datapath+'"'+ \ ' prodname="unit-test"'+ \ - ' infile="runlist_cmd1.dat"'+ \ + ' infile="'+self._runlist+'"'+ \ ' bkgpars=1'+\ - ' outobs="obs_iactobs_cmd1.xml"'+ \ - ' outmodel="bgd_iactobs_cmd1.xml"'+ \ + ' outobs="csiactobs_obs_cmd1.xml"'+ \ + ' outmodel="csiactobs_bgd_cmd1.xml"'+ \ ' logfile="csiactobs_cmd1.log" chatter=1' # Check if execution of wrong command fails @@ -89,23 +90,36 @@ class Test(test): 'Check successful execution from command line') # Check observation definition XML file - self._check_obsdef('obs_iactobs_cmd1.xml') + self._check_obsdef('csiactobs_obs_cmd1.xml') # Check model definition XML file - self._check_moddef('bgd_iactobs_cmd1.xml') + self._check_moddef('csiactobs_bgd_cmd1.xml') # Setup csiactobs command cmd = csiactobs+' datapath="data_path_that_does_not_exist"'+ \ ' prodname="unit-test"'+ \ - ' infile="runlist_cmd1.dat"'+ \ + ' infile="'+self._runlist+'"'+ \ ' bkgpars=1'+\ - ' outobs="obs_iactobs_cmd2.xml"'+ \ - ' outmodel="bgd_iactobs_cmd2.xml"'+ \ + ' outobs="csiactobs_obs_cmd2.xml"'+ \ + ' outmodel="csiactobs_bgd_cmd2.xml"'+ \ ' logfile="csiactobs_cmd2.log" chatter=1' # Check if execution failed self.test_assert(self._execute(cmd) != 0, - 'Check invalid input file when executed from command line') + 'Check invalid input datapath when executed from command line') + + # Setup csiactobs command + cmd = csiactobs+' datapath="'+self._datapath+'"'+ \ + ' prodname="unit-test-doesnt-exist"'+ \ + ' infile="'+self._runlist+'"'+ \ + ' bkgpars=1'+\ + ' outobs="csiactobs_obs_cmd3.xml"'+ \ + ' outmodel="csiactobs_bgd_cmd3.xml"'+ \ + ' logfile="csiactobs_cmd3.log" chatter=1' + + # Check if execution failed + self.test_assert(self._execute(cmd) != 0, + 'Check invalid input prodname when executed from command line') # Return return @@ -119,10 +133,10 @@ class Test(test): iactobs = cscripts.csiactobs() iactobs['datapath'] = self._datapath iactobs['prodname'] = 'unit-test' - iactobs['infile'] = 'runlist_py1.dat' + iactobs['infile'] = self._runlist iactobs['bkgpars'] = 1 - iactobs['outobs'] = 'obs_iactobs_py1.xml' - iactobs['outmodel'] = 'bgd_iactobs_py1.xml' + iactobs['outobs'] = 'csiactobs_obs_py1.xml' + iactobs['outmodel'] = 'csiactobs_bgd_py1.xml' iactobs['logfile'] = 'csiactobs_py1.log' iactobs['chatter'] = 2 @@ -132,10 +146,88 @@ class Test(test): iactobs.save() # Check observation definition XML file - self._check_obsdef('obs_iactobs_py1.xml') + self._check_obsdef('csiactobs_obs_py1.xml') # Check model definition XML file - self._check_moddef('bgd_iactobs_py1.xml') + self._check_moddef('csiactobs_bgd_py1.xml') + + # Create test runlist + runlist = ['0'] + + # Set-up csiactobs using a runlist with 2 background parameters + iactobs = cscripts.csiactobs() + iactobs['datapath'] = self._datapath + iactobs['prodname'] = 'unit-test' + iactobs['bkgpars'] = 2 + iactobs['outobs'] = 'csiactobs_obs_py2.xml' + iactobs['outmodel'] = 'csiactobs_bgd_py2.xml' + iactobs['logfile'] = 'csiactobs_py2.log' + iactobs['chatter'] = 3 + iactobs.runlist(runlist) + + # Run csiactobs script and save run list + iactobs.logFileOpen() # Make sure we get a log file + iactobs.run() + iactobs.save() + + # Test return functions + self.test_value(iactobs.obs().size(), 1, + 'Check number of observations in container') + self.test_value(iactobs.ebounds().size(), 0, + 'Check number of energy boundaries') + + # Check observation definition XML file + self._check_obsdef('csiactobs_obs_py2.xml') + + # Check model definition XML file + self._check_moddef('csiactobs_bgd_py2.xml') + + # Set-up csiactobs with a large number of free parameters and "aeff" + # background + iactobs = cscripts.csiactobs() + iactobs['datapath'] = self._datapath + iactobs['prodname'] = 'unit-test' + iactobs['infile'] = self._runlist + iactobs['bkgpars'] = 8 + iactobs['bkg_mod_hiera'] = 'aeff' + iactobs['outobs'] = 'csiactobs_obs_py3.xml' + iactobs['outmodel'] = 'csiactobs_bgd_py3.xml' + iactobs['logfile'] = 'csiactobs_py3.log' + iactobs['chatter'] = 4 + + # Execute csiactobs script + iactobs.execute() + + # Check observation definition XML file + self._check_obsdef('csiactobs_obs_py3.xml') + + # Check model definition XML file + self._check_moddef('csiactobs_bgd_py3.xml') + + # Set-up csiactobs with a "gauss" background and "inmodel" parameter + iactobs = cscripts.csiactobs() + iactobs['datapath'] = self._datapath + iactobs['inmodel'] = self._model + iactobs['prodname'] = 'unit-test' + iactobs['infile'] = self._runlist + iactobs['bkgpars'] = 1 + iactobs['bkg_mod_hiera'] = 'gauss' + iactobs['outobs'] = 'NONE' + iactobs['outmodel'] = 'NONE' + iactobs['logfile'] = 'csiactobs_py4.log' + iactobs['chatter'] = 4 + + # Run csiactobs script + iactobs.logFileOpen() # Make sure we get a log file + iactobs.run() + + # Check number of observations + self.test_value(iactobs.obs().size(), 1, + 'Check number of observations in container') + + # Check number of models + self.test_value(iactobs.obs().models().size(), 3, + 'Check number of models in container') # Return return @@ -154,34 +246,26 @@ class Test(test): # Check number of observations self.test_value(obs.size(), 1, 'Check for single observation in XML file') - self.test_assert(obs[0].eventfile().file() == 'events.fits.gz', - 'Check event file name') - self.test_assert(obs[0].eventfile().extname() == 'EVENTS', - 'Check event extension name') - self.test_assert(rsp.aeff().filename().file() == - 'irf_file.fits.gz', - 'Check effective area name') - self.test_assert(rsp.aeff().filename().extname() == - 'EFFECTIVE AREA', - 'Check effective area extension name') - self.test_assert(rsp.psf().filename().file() == - 'irf_file.fits.gz', - 'Check point spread function name') - self.test_assert(rsp.psf().filename().extname() == - 'POINT SPREAD FUNCTION', - 'Check point spread function extension name') - self.test_assert(rsp.edisp().filename().file() == - 'irf_file.fits.gz', - 'Check energy dispersion name') - self.test_assert(rsp.edisp().filename().extname() == - 'ENERGY DISPERSION', - 'Check energy dispersion extension name') - self.test_assert(rsp.background().filename().file() == - 'irf_file.fits.gz', - 'Check background name') - self.test_assert(rsp.background().filename().extname() == - 'BACKGROUND', - 'Check background extension name') + self.test_value(obs[0].eventfile().file(), 'events.fits.gz', + 'Check event file name') + self.test_value(obs[0].eventfile().extname(), 'EVENTS', + 'Check event extension name') + self.test_value(rsp.aeff().filename().file(), 'irf_file.fits.gz', + 'Check effective area file name') + self.test_value(rsp.aeff().filename().extname(), 'EFFECTIVE AREA', + 'Check effective area extension name') + self.test_value(rsp.psf().filename().file(), 'irf_file.fits.gz', + 'Check point spread function file name') + self.test_value(rsp.psf().filename().extname(), 'POINT SPREAD FUNCTION', + 'Check point spread function extension name') + self.test_value(rsp.edisp().filename().file(), 'irf_file.fits.gz', + 'Check energy dispersion file name') + self.test_value(rsp.edisp().filename().extname(), 'ENERGY DISPERSION', + 'Check energy dispersion extension name') + self.test_value(rsp.background().filename().file(), 'irf_file.fits.gz', + 'Check background file name') + self.test_value(rsp.background().filename().extname(), 'BACKGROUND', + 'Check background extension name') # Return return