#!/usr/bin/env python # ========================================================================== # Generates an IACT observation definition XML file. # # Copyright (C) 2015-2017 Michael Mayer # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import sys import os import json import gammalib import ctools # =============== # # csiactobs class # # =============== # class csiactobs(ctools.cscript): """ Generates an IACT observation definition XML file This class implements the creation of a observation xml file for IACT data analysis. This class is dedicated for use inside a IACT Collaboration, i.e. it can only be used if you have access to IACT data in FITS format. The FITS data has to be structured and in the format described here: http://gamma-astro-data-formats.readthedocs.org/en/latest/ """ # Constructor def __init__(self, *argv): """ Constructor """ # Set name and version self._name = 'csiactobs' self._version = ctools.__version__ # Initialise some members self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() self._datapath = os.getenv('VHEFITS','') self._prodname = '' self._xml = gammalib.GXml() self._models = gammalib.GModels() self._runlist = [] self._runlistfile = gammalib.GFilename() self._bkgpars = 0 self._master_indx = '' self._use_bkg_scale = False self._ev_hiera = [''] self._aeff_hiera = [''] self._psf_hiera = [''] self._bkg_hiera = [''] self._edisp_hiera = [''] self._bkg_mod_hiera = [''] self._bkg_gauss_norm = 1.0 self._bkg_gauss_index = 0.0 self._bkg_gauss_sigma = 1.0 self._bkg_aeff_index = 0.0 self._bkg_aeff_norm = 1.0 self._bkg_range_factor = 1.0 self._hdu_index = '' self._obs_index = '' self._subdir = '' self._debug = False # Initialise empty observation definition XML file self._xml.append(gammalib.GXmlElement('observation_list ' 'title="observation list"')) # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Append an observation list to XML instance self._xml.append(gammalib.GXmlElement('observation_list title="observation list"')) self._xml_obslist = self._xml.element('observation_list', 0) # Return return # Private methods def _get_parameters(self): """ Get parameters from parfile and setup the observation """ # Set data path if self._datapath == '': self._datapath = self['datapath'].string() # Query input parameters self['inmodel'].filename() # Expand environment self._datapath = gammalib.expand_env(self._datapath) # Read FITS production self._prodname = self['prodname'].string() # Read runlist file if list not already filled if len(self._runlist) == 0: # Get file name self._runlistfile = self['infile'].filename() # Read runlist from file runfile = open(self._runlistfile.url()) for line in runfile.readlines(): if len(line) == 0: continue if line[0] == '#': continue if len(line.split()) > 0: self._runlist.append(line.split()[0]) runfile.close() # Read number of background parameters 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() # Read flag for background scaling factor self._use_bkg_scale = self['bkg_scale'].boolean() # Read hierarchy of file loading self._ev_hiera = self['ev_hiera'].string().split('|') self._aeff_hiera = self['aeff_hiera'].string().split('|') self._psf_hiera = self['psf_hiera'].string().split('|') self._bkg_hiera = self['bkg_hiera'].string().split('|') self._edisp_hiera = self['edisp_hiera'].string().split('|') self._bkg_mod_hiera = self['bkg_mod_hiera'].string().split('|') # Read hidden background parameters self._bkg_gauss_norm = self['bkg_gauss_norm'].real() self._bkg_gauss_index = self['bkg_gauss_index'].real() self._bkg_gauss_sigma = self['bkg_gauss_sigma'].real() self._bkg_aeff_index = self['bkg_aeff_index'].real() self._bkg_aeff_norm = self['bkg_aeff_norm'].real() self._bkg_range_factor = self['bkg_range_factor'].real() # Open master index file and look for prodname master_file = os.path.join(self._datapath, self._master_indx) if not os.path.isfile(master_file): raise RuntimeError('FITS data store not available. No master index file found. Make sure the file is copied from the server and your datapath is set correctly.') # Open and load JSON file json_data = open(master_file).read() data = json.loads(json_data) if not 'datasets' in data: raise RuntimeError('Key "datasets" not available in master index file.') # Get configurations configs = data['datasets'] # Initialise HDUs self._hdu_index = self._obs_index = '' # 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 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 self._subdir = os.path.dirname(self._hdu_index) self._debug = False # Debugging in client tools # Write input parameters into logger self._log_parameters(gammalib.TERSE) # Return return def _background_spectrum(self, prefactor, index, emin = 0.01, emax = 100.0): """ Create a background spectrum model dependent on user parameters Parameters ---------- prefactor : float Power law prefactor of spectral model index : float Power law index 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) 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: spec['Normalization'].free() else: # Create power law model if self._bkgpars <= 2: # 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) spec[0].max(prefactor * self._bkg_range_factor) spec[1].scale(1) spec[1].min(-5.0) spec[1].max(5.0) # Set number of free parameters if self._bkgpars == 0: spec[0].fix() spec[1].fix() elif self._bkgpars == 1: spec[0].free() spec[1].fix() else: spec[0].free() spec[1].free() else: # Create reference powerlaw 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) # 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) par.min(value / self._bkg_range_factor) par.max(value * self._bkg_range_factor) # Return spectrum return spec def _iact_background(self, telescope, obs_id, bkg_scale, bkgtype, emin=0.01, emax=100): """ Create an IACT background model Parameters ---------- telescope : str Name of telescope obs_id : str Observation ID bkg_scale : float Background scaling factor bkgtype : str 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 # Create background spectrum spec = self._background_spectrum(prefactor, index, emin, emax) # Create background model bck = gammalib.GCTAModelIrfBackground(spec) # Set AeffBackground elif bkgtype == 'aeff': # Set background components prefactor = bkg_scale * self._bkg_aeff_norm spec = self._background_spectrum(prefactor, self._bkg_aeff_index, emin, emax) # Create background model bck = gammalib.GCTAModelAeffBackground(spec) # Set Gaussian Background elif bkgtype == 'gauss': # Set background components prefactor = bkg_scale * self._bkg_gauss_norm spec = self._background_spectrum(prefactor, self._bkg_gauss_index, emin, emax) radial = gammalib.GCTAModelRadialGauss(self._bkg_gauss_sigma) # Create background model bck = gammalib.GCTAModelRadialAcceptance(radial, spec) else: msg = 'Background type "'+bkgtype+'" unsupported' raise RuntimeError(msg) # Copy model model = bck.clone() # Assign specific run id model.ids(str(obs_id)) # Assign instrument model.instruments(telescope) # Set name (arbitrary) model.name('bkg_'+str(obs_id)) # Turn off TS calculation for background model model.tscalc(False) # Return model return model def _append_inmodels(self): """ Append input models """ # If there are models provided by "inmodels" then append them now to # the model container filename = self['inmodel'].filename().url() if filename != '' and filename != 'NONE': # Log header and input model file self._log_header1(gammalib.TERSE, 'Append input models') self._log_value(gammalib.NORMAL, 'Input model file', filename) # Load input models models = gammalib.GModels(filename) # Loop over input models and append them to the model container for model in models: self._models.append(model) self._log_value(gammalib.NORMAL, 'Append model', '"'+model.name()+'"') # Return return def _write_summary(self): """ Write observation summary """ # Log header self._log_header1(gammalib.NORMAL, 'Observation summary') # Set energy range dependent on whether boundaries exist or # not if self._ebounds.size() > 0: erange = '%s - %s' % (str(self._ebounds.emin()), str(self._ebounds.emax())) else: erange = 'not available' # Log energy range self._log_value(gammalib.NORMAL, 'Energy range', erange) # Return return def _get_filename(self, hdu, hierarchy, formats): """ Retrieves a filename from a hdu index Parameters ---------- hdu : `~gammalib.GFitsTable` HDU FITS table hierarchy : list of str List of strings containing the hierarchy how to look for files formats : list of str File formats available for this observation Returns ------- filename : str Filename of requested file """ # Initialise filename filename = '' # Handle hierarchies index = -1 # Loop over possible formats for version in hierarchy: # Count number of same available versions n = formats.count(version) # Leave loop and store index of first occurence if n > 0: index = formats.index(version) break # Build filename from information stored in hduindx_hdu if index >= 0: filename = os.path.join(self._subdir, hdu['FILE_DIR'][index], hdu['FILE_NAME'][index]) filehdu = hdu['HDU_NAME'][index] filename += '['+filehdu+']' # Return filename return filename def _is_present(self, filename, obs_id, filetype): """ Checks if a filename is present Parameters ---------- filename : str Filename obs_id : int Observation ID filetype : str Type of file Returns ------- present : bool True if filename is present, False otherwise """ # Initialise return value present = True # Create filename instance fname = gammalib.GFilename(filename) # Check if file is empty if fname.is_empty(): msg = 'Skipping observation "%s": No %s found' % (obs_id, filetype) self._log_string(gammalib.NORMAL, msg) present = False # Check if file exists elif not fname.exists(): msg = ('Skipping observation "%s": %s "%s" does not exist' % (obs_id, filetype, filename)) self._log_string(gammalib.NORMAL, msg) present = False # Return present flag return present # Public methods def run(self): """ Run the script """ # Switch screen logging on in debug mode if self._logDebug(): self._log.cout(True) # Get parameters self._get_parameters() # Clear models self._models.clear() # Log output msg = 'Looping over %d runs' % len(self._runlist) self._log_header1(gammalib.TERSE, msg) # Initialise energy range values for logging self._ebounds.clear() # Loop over runs for obs_id in self._runlist: # Create selection string obs_selection = '[OBS_ID=='+str(obs_id)+']' # Open HDU index file hduindx = gammalib.GFits(self._hdu_index+'[HDU_INDEX]'+ obs_selection) hduindx_hdu = hduindx['HDU_INDEX'] # Initialise types and formats formats = [] for i in range(hduindx_hdu.nrows()): formats.append(hduindx_hdu['HDU_CLASS'][i]) # Retrieve filenames eventfile = self._get_filename(hduindx_hdu, self._ev_hiera, formats) aefffile = self._get_filename(hduindx_hdu, self._aeff_hiera, formats) psffile = self._get_filename(hduindx_hdu, self._psf_hiera, formats) edispfile = self._get_filename(hduindx_hdu, self._edisp_hiera, formats) bkgfile = self._get_filename(hduindx_hdu, self._bkg_hiera, formats) # Check for presence of required event file if not self._is_present(eventfile, obs_id, "event file"): continue # Check for presence of required effective area file if not self._is_present(aefffile, obs_id, "effective area file"): continue # Check for presence of required psf file if not self._is_present(psffile, obs_id, "PSF file"): continue # Check for optional energy dispersion file if not gammalib.GFilename(edispfile).exists(): # Print warning that edisp cannot be used for this observation msg = ('Warning: observation "%s" has no energy dispersion ' 'information' % obs_id) self._log_string(gammalib.NORMAL, msg) # Set energy dispersion file as empty edispfile = '' # Create a copy of background model hierarchy for this run run_bkg_mod_hierarchy = list(self._bkg_mod_hiera) # Check if background file is available # Remove IRF background from hierarchy if file doesnt exist if not gammalib.GFilename(bkgfile).exists(): # Set background file as empty bkgfile = '' # Check for IRF background in background model hierarchy if 'irf' in run_bkg_mod_hierarchy: # Remove IRF background if file doesnt exist run_bkg_mod_hierarchy.remove('irf') # Print warning that edisp cannot be used for this # observation msg = ('Warning: observation "%s" has no background ' 'information (IRF background cannot be used)' % obs_id) self._log_string(gammalib.NORMAL, msg) # Check if background model hierarchy is empty if len(run_bkg_mod_hierarchy) == 0: # Skip observation if no background can be used msg = ('Skipping observation "%s": No background can be ' 'used' % obs_id) self._log_string(gammalib.NORMAL, msg) continue else: # Log if we fall back to next background approach msg = ('Observation "%s": Falling back to background "%s"' % (obs_id, run_bkg_mod_hierarchy[0])) self._log_string(gammalib.NORMAL, msg) # Close hdu index file hduindx.close() # 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 event FITS file fits = gammalib.GFits(eventfile) # Get object and telescope strings from event hdu events = fits[gammalib.GFilename(eventfile).extname()] object_name = events.string('OBJECT') telescope = events.string('TELESCOP') # Close FITS file fits.close() # Open effective area file to look for threshold aeff_fits = gammalib.GFits(aefffile) aeff_table = aeff_fits[gammalib.GFilename(aefffile).extname()] # Set energy range from header keyword if present if aeff_table.has_card('LO_THRES') and aeff_table.has_card('HI_THRES'): # Get safe energy range run_emin = gammalib.GEnergy(aeff_table.real('LO_THRES'), 'TeV') run_emax = gammalib.GEnergy(aeff_table.real('HI_THRES'), 'TeV') # Append to ebounds self._ebounds.append(run_emin, run_emax) else: # Set default values for energy range run_emin = gammalib.GEnergy(10,'GeV') run_emax = gammalib.GEnergy(100,'TeV') # Close Aeff fits file aeff_fits.close() # Append instrumental background model self._models.append(self._iact_background(telescope, obs_id, bkg_scale, run_bkg_mod_hierarchy[0], run_emin.TeV(), run_emax.TeV())) # Append observation to XML and set attributes obs = self._xml_obslist.append('observation') obs.attribute('name', object_name) obs.attribute('id', str(obs_id)) obs.attribute('instrument', telescope) # Append event file ev = gammalib.GXmlElement('parameter name="EventList"') ev.attribute('file', eventfile) # Append effective area aeff = gammalib.GXmlElement('parameter name="EffectiveArea"') aeff.attribute('file', aefffile) # Append PSF psf = gammalib.GXmlElement('parameter name="PointSpreadFunction"') psf.attribute('file', psffile) # Append energy dispersion edisp = gammalib.GXmlElement('parameter name="EnergyDispersion"') edisp.attribute('file', edispfile) # Append background bck = gammalib.GXmlElement('parameter name="Background"') bck.attribute('file', bkgfile) # Assign events and IRFs to observation obs.append(ev) obs.append(aeff) obs.append(psf) obs.append(edisp) obs.append(bck) # Log the observation ID and object name that has been appended self._log_value(gammalib.NORMAL, 'Append observation', '%s ("%s")' % (obs_id, object_name)) # Log the file names of the event and response files self._log_value(gammalib.EXPLICIT, ' Event file', eventfile) self._log_value(gammalib.EXPLICIT, ' Effective area file', aefffile) self._log_value(gammalib.EXPLICIT, ' Point spread function file', psffile) self._log_value(gammalib.EXPLICIT, ' Energy dispersion file', edispfile) self._log_value(gammalib.EXPLICIT, ' Background file', bkgfile) if self._use_bkg_scale: self._log_value(gammalib.EXPLICIT, ' Background scale', bkg_scale) # Append models provided by "inmodel" to the model container self._append_inmodels() # Write observation summary self._write_summary() # Write warning in no observation could be used from runlist if self._xml_obslist.size() == 0: self._log_string(gammalib.NORMAL, 'WARNING: No observation was ' 'appended from specified runlist') # Return return def execute(self): """ Execute the script """ # Open logfile self.logFileOpen() # Read ahead output parameters self._read_ahead(True) # Run the script self.run() # Save residual map self.save() # Return return def save(self): """ Save observation definition and model definition XML file """ # Write header self._log_header1(gammalib.TERSE, 'Save XML files') # Get filenames outobs = self['outobs'].filename() outmodel = self['outmodel'].filename() # Log observation definition file name self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs.url()) # Save observation definition XML file self._xml.save(outobs.url()) # Log model definition file name self._log_value(gammalib.NORMAL, 'Model definiton XML file', outmodel.url()) # Save model XML file 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 """ # Return energy boundaries return self._ebounds def obs(self): """ Return observations container """ # Clear observations self._obs.clear() # Read observations if XML is filled if self._xml.size(): self._obs.read(self._xml) # Assign models self._obs.models(self._models) # Return observations return self._obs # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': # Create instance of application app = csiactobs(sys.argv) # Execute application app.execute()