33 Generates an IACT observation definition XML file
35 This class implements the creation of a observation xml file for IACT
36 data analysis. This class is dedicated for use inside a IACT
37 Collaboration, i.e. it can only be used if you have access to IACT data
38 in FITS format. The FITS data has to be structured and in the format
40 http://gamma-astro-data-formats.readthedocs.org/en/latest/
49 self._init_cscript(self.__class__.__name__, ctools.__version__, argv)
52 self.
_obs = gammalib.GObservations()
56 self.
_xml = gammalib.GXml()
81 self._xml.append(gammalib.GXmlElement(
'observation_list '
82 'title="observation list"'))
85 self._xml.append(gammalib.GXmlElement(
'observation_list title="observation list"'))
95 Get parameters from parfile and setup the observation
99 datapath = self[
'datapath'].string()
100 if gammalib.toupper(datapath) !=
'NONE':
109 self[
'inmodel'].query()
112 self.
_prodname = self[
'prodname'].string()
121 runfile = open(self._runlistfile.url())
122 for line
in runfile.readlines():
127 if len(line.split()) > 0:
128 self._runlist.append(line.split()[0])
132 self.
_bkgpars = self[
'bkgpars'].integer()
135 if self._read_ahead():
136 self[
'outmodel'].filename()
137 self[
'outobs'].filename()
146 self.
_ev_hiera = self[
'ev_hiera'].string().split(
'|')
147 self.
_aeff_hiera = self[
'aeff_hiera'].string().split(
'|')
148 self.
_psf_hiera = self[
'psf_hiera'].string().split(
'|')
149 self.
_bkg_hiera = self[
'bkg_hiera'].string().split(
'|')
150 self.
_edisp_hiera = self[
'edisp_hiera'].string().split(
'|')
163 if not os.path.isfile(master_file):
164 raise RuntimeError(
'FITS data store not available. No master '
165 'index file found. Make sure the file is '
166 'copied from the server and your datapath '
170 json_data = open(master_file).read()
171 data = json.loads(json_data)
172 if not 'datasets' in data:
173 raise RuntimeError(
'Key "datasets" not available in master '
177 configs = data[
'datasets']
183 for config
in configs:
197 raise RuntimeError(
'*** ERROR: FITS data store "'+self.
_prodname+
198 '" not available. Run csiactdata to get a list '
199 'of available storage names')
202 filename = gammalib.GFilename(self.
_hdu_index+
'[HDU_INDEX]')
203 if not filename.is_fits():
204 raise RuntimeError(
'*** ERROR: HDU index file "'+self.
_hdu_index+
205 '[HDU_INDEX]" for FITS data store "'+
206 self.
_prodname+
'" not available. Check your '
207 'master index file or run csiactdata to get '
208 'a list of available storage names.')
215 filename = gammalib.GFilename(self.
_obs_index+
'[OBS_INDEX]')
218 if filename.is_fits():
225 if not fits[
'OBS_INDEX'].contains(
'BKG_SCALE'):
240 self._log_parameters(gammalib.TERSE)
247 Create a background spectrum model dependent on user parameters
252 Power law prefactor of spectral model
254 Power law index of spectral model
256 Minimum energy (in case a spectral node function is required)
258 Maximum energy (in case a spectral node function is required)
262 spec : `~gammalib.GModelSpectral()`
263 Spectral model for the background shape
266 if index == 0.0
and self.
_bkgpars <= 1:
267 spec = gammalib.GModelSpectralConst()
272 spec[
'Normalization'].value(prefactor)
276 spec[
'Normalization'].fix()
278 spec[
'Normalization'].free()
286 pivot = gammalib.GEnergy(1.0,
'TeV')
287 spec = gammalib.GModelSpectralPlaw(prefactor, index, pivot)
310 pivot = gammalib.GEnergy(1.0,
'TeV')
311 plaw = gammalib.GModelSpectralPlaw(prefactor, index, pivot)
314 spec = gammalib.GModelSpectralNodes()
317 bounds = gammalib.GEbounds(self.
_bkgpars,
318 gammalib.GEnergy(emin,
'TeV'),
319 gammalib.GEnergy(emax,
'TeV'))
322 for i
in range(bounds.size()):
323 energy = bounds.elogmean(i)
324 value = plaw.eval(energy)
327 spec.append(energy, value)
333 if 'Energy' in par.name():
337 elif 'Intensity' in par.name():
347 emin=0.01, emax=100):
349 Create an IACT background model
358 Background scaling factor
360 Type of background (irf,aeff, or gauss)
364 model : `~gammalib.GModelData()`
365 Background model for IACT observation
373 prefactor *= bkg_scale
379 bck = gammalib.GCTAModelIrfBackground(spec)
382 elif bkgtype ==
'aeff':
391 bck = gammalib.GCTAModelAeffBackground(spec)
394 elif bkgtype ==
'gauss':
404 bck = gammalib.GCTAModelRadialAcceptance(radial, spec)
407 msg =
'Background type "'+bkgtype+
'" unsupported'
408 raise RuntimeError(msg)
414 model.ids(str(obs_id))
417 model.instruments(telescope)
420 model.name(
'bkg_'+str(obs_id))
434 if self[
'inmodel'].is_valid():
437 filename = self[
'inmodel'].filename().url()
440 self._log_header1(gammalib.TERSE,
'Append input models')
441 self._log_value(gammalib.NORMAL,
'Input model file', filename)
444 models = gammalib.GModels(filename)
448 self._models.append(model)
449 self._log_value(gammalib.NORMAL,
'Append model',
450 '"'+model.name()+
'"')
457 Write observation summary
460 self._log_header1(gammalib.NORMAL,
'Observation summary')
464 if self._ebounds.size() > 0:
465 erange =
'%s - %s' % (str(self._ebounds.emin()),
466 str(self._ebounds.emax()))
468 erange =
'not available'
471 self._log_value(gammalib.NORMAL,
'Energy range', erange)
478 Retrieves a filename from a hdu index
482 hdu : `~gammalib.GFitsTable`
484 hierarchy : list of str
485 List of strings containing the hierarchy how to look for files
486 formats : list of str
487 File formats available for this observation
492 Filename of requested file
502 for version
in hierarchy:
505 n = formats.count(version)
509 index = formats.index(version)
514 filename = os.path.join(self.
_subdir,
515 hdu[
'FILE_DIR'][index],
516 hdu[
'FILE_NAME'][index])
517 filehdu = hdu[
'HDU_NAME'][index]
518 filename +=
'['+filehdu+
']'
525 Checks if a filename is present
539 True if filename is present, False otherwise
546 fname = gammalib.GFilename(filename)
550 msg =
'Skipping observation "%s": No %s found' % (obs_id, filetype)
551 self._log_string(gammalib.NORMAL, msg)
555 elif not fname.exists():
556 msg = (
'Skipping observation "%s": %s "%s" does not exist' %
557 (obs_id, filetype, filename))
558 self._log_string(gammalib.NORMAL, msg)
577 msg =
'Looping over %d runs' % len(self.
_runlist)
578 self._log_header1(gammalib.TERSE, msg)
581 self._ebounds.clear()
587 obs_selection =
'[OBS_ID=='+str(obs_id)+
']'
590 hduindx = gammalib.GFits(self.
_hdu_index+
'[HDU_INDEX]'+
592 hduindx_hdu = hduindx[
'HDU_INDEX']
596 for i
in range(hduindx_hdu.nrows()):
597 formats.append(hduindx_hdu[
'HDU_CLASS'][i])
607 if not self.
_is_present(eventfile, obs_id,
"event file"):
611 if not self.
_is_present(aefffile, obs_id,
"effective area file"):
615 if not self.
_is_present(psffile, obs_id,
"PSF file"):
619 if not gammalib.GFilename(edispfile).exists():
622 msg = (
'Warning: observation "%s" has no energy dispersion '
623 'information' % obs_id)
624 self._log_string(gammalib.NORMAL, msg)
634 if not gammalib.GFilename(bkgfile).exists():
640 if 'irf' in run_bkg_mod_hierarchy:
643 run_bkg_mod_hierarchy.remove(
'irf')
647 msg = (
'Warning: observation "%s" has no background '
648 'information (IRF background cannot be used)' %
650 self._log_string(gammalib.NORMAL, msg)
653 if len(run_bkg_mod_hierarchy) == 0:
656 msg = (
'Skipping observation "%s": No background can be '
658 self._log_string(gammalib.NORMAL, msg)
663 msg = (
'Observation "%s": Falling back to background "%s"' %
664 (obs_id, run_bkg_mod_hierarchy[0]))
665 self._log_string(gammalib.NORMAL, msg)
677 obsindx = gammalib.GFits(self.
_obs_index+
'[OBS_INDEX]'+
679 bkg_scale = obsindx[
'OBS_INDEX'][
'BKG_SCALE'][0]
683 fits = gammalib.GFits(eventfile)
686 events = fits[gammalib.GFilename(eventfile).extname()]
687 object_name = events.string(
'OBJECT')
688 telescope = events.string(
'TELESCOP')
694 aeff_fits = gammalib.GFits(aefffile)
695 aeff_table = aeff_fits[gammalib.GFilename(aefffile).extname()]
698 if aeff_table.has_card(
'LO_THRES')
and aeff_table.has_card(
'HI_THRES'):
701 run_emin = gammalib.GEnergy(aeff_table.real(
'LO_THRES'),
'TeV')
702 run_emax = gammalib.GEnergy(aeff_table.real(
'HI_THRES'),
'TeV')
705 self._ebounds.append(run_emin, run_emax)
709 run_emin = gammalib.GEnergy(10,
'GeV')
710 run_emax = gammalib.GEnergy(100,
'TeV')
719 run_bkg_mod_hierarchy[0],
724 obs = self._xml_obslist.append(
'observation')
725 obs.attribute(
'name', object_name)
726 obs.attribute(
'id', str(obs_id))
727 obs.attribute(
'instrument', telescope)
730 ev = gammalib.GXmlElement(
'parameter name="EventList"')
731 ev.attribute(
'file', eventfile)
734 aeff = gammalib.GXmlElement(
'parameter name="EffectiveArea"')
735 aeff.attribute(
'file', aefffile)
738 psf = gammalib.GXmlElement(
'parameter name="PointSpreadFunction"')
739 psf.attribute(
'file', psffile)
742 edisp = gammalib.GXmlElement(
'parameter name="EnergyDispersion"')
743 edisp.attribute(
'file', edispfile)
746 bck = gammalib.GXmlElement(
'parameter name="Background"')
747 bck.attribute(
'file', bkgfile)
757 self._log_value(gammalib.NORMAL,
'Append observation',
'%s ("%s")' %
758 (obs_id, object_name))
761 self._log_value(gammalib.EXPLICIT,
' Event file', eventfile)
762 self._log_value(gammalib.EXPLICIT,
' Effective area file', aefffile)
763 self._log_value(gammalib.EXPLICIT,
' Point spread function file',
765 self._log_value(gammalib.EXPLICIT,
' Energy dispersion file',
767 self._log_value(gammalib.EXPLICIT,
' Background file', bkgfile)
769 self._log_value(gammalib.EXPLICIT,
' Background scale',
779 if self._xml_obslist.size() == 0:
780 self._log_string(gammalib.NORMAL,
'WARNING: No observation was '
781 'appended from specified runlist')
788 Save observation definition and model definition XML file
791 self._log_header1(gammalib.TERSE,
'Save XML files')
794 outobs = self[
'outobs'].filename()
795 outmodel = self[
'outmodel'].filename()
798 self._log_value(gammalib.NORMAL,
'Obs. definition XML file',
802 self._xml.save(outobs.url())
805 self._log_value(gammalib.NORMAL,
'Model definiton XML file',
809 self._models.save(outmodel.url())
821 List of observation IDs
826 Input runlist is not a Python list
835 raise RuntimeError(
'Argument is not a Python list')
842 Return runlist energy boundaries
849 Return observations container
856 self._obs.read(self.
_xml)
868 if __name__ ==
'__main__':