From b17eaf02157e070096055cd99bb7bfe95174f98d Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 16 Oct 2017 18:16:56 +0200 Subject: [PATCH] Integrate csphagen script (#2208) --- AUTHORS | 1 + ChangeLog | 3 ++- NEWS | 5 ++++- cscripts/csobsselect.py | 2 +- cscripts/csphagen.par | 19 ++++++++++--------- cscripts/csphagen.py | 310 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------- doc/source/admin/release_history/1.5.rst | 3 +++ doc/source/users/reference_manual/csphagen.rst | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------- test/test_csphagen.py | 223 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------- 9 files changed, 387 insertions(+), 256 deletions(-) diff --git a/AUTHORS b/AUTHORS index daee14b..42793a0 100644 --- a/AUTHORS +++ b/AUTHORS @@ -15,4 +15,5 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de) Stefan Ohm (stefan.ohm81@gmail.com) David Sanchez (david.sanchez@lapp.in2p3.fr) Hubert Siejkowski (h.siejkowski@cyfronet.pl) +Luigi Tibaldo (luigi.tibaldo@irap.omp.eu) Alexander Ziegler (alexander.ziegler@fau.de) diff --git a/ChangeLog b/ChangeLog index 65bf1cd..223b363 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-08-08 +2017-10-16 * Version 1.5.0 released ======================== + Add csphagen script (#2208) Use t-type parameters for all "tmin" and "tmax" parameters (#1864) diff --git a/NEWS b/NEWS index 2940dcf..1a26f1d 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.5.0 -8 August 2017 +16 October 2017 Introduction @@ -122,6 +122,9 @@ cscripts Use t-type parameters for "tmin" and "tmax" to all specifying times in various formats (#1864). +The script csphagen was added to allow the generation of PHA, ARF and RMF datasets +for classical On/Off IACT analysis (#2208). + Examples -------- diff --git a/cscripts/csobsselect.py b/cscripts/csobsselect.py index 54f66e2..523e539 100755 --- a/cscripts/csobsselect.py +++ b/cscripts/csobsselect.py @@ -320,7 +320,7 @@ class csobsselect(ctools.cscript): # Otherwise log filename and save file else: # Log filename - self._log_value(gammalib.NORMAL, 'Observation definition XML file', + self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs.url()) # Save observations diff --git a/cscripts/csphagen.par b/cscripts/csphagen.par index 396b032..4cb633e 100644 --- a/cscripts/csphagen.par +++ b/cscripts/csphagen.par @@ -16,11 +16,12 @@ # # Input and output parameters #============================ -inobs, f, a, obs.xml,,, "Input event list or observation definition XML file" -caldb, s, a, prod2,,, "Calibration database" -irf, s, a, South_0.5h,,, "Instrument response function" -exclusion, f, h, NONE,,, "FITS with WCS map of exclusion regions" -outroot, s, a, spectrum,,, "Root of the name for output PHA, ARF, RMF, XML, and ds9 reg files" +inobs, f, a, obs.xml,,, "Input event list or observation definition XML file" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" +inexclusion, f, h, NONE,,, "Sky map of exclusion regions" +outobs, s, a, onoff_obs.xml,,, "Output observation definition XML file" +prefix, s, h, onoff,,, "Prefix for output files" # # Script parameters @@ -30,7 +31,6 @@ emax, r, a, 100.0,,, "Upper energy limit (TeV)" enumbins, i, a, 120,,, "Number of energy bins" ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm" ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" -bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation" srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region" coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" ra, r, a, 83.63,0,360, "Right Ascension of source region centre (deg)" @@ -38,9 +38,10 @@ dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)" glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)" glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)" rad, r, a, 0.2,0,180, "Radius of source region circle (deg)" -bkgregmin, i, h, 2,1,1000, "Minimum number of background region for REFLECTED method" -maxoffset, r, h, 4.,,, "Maximum offset of source from camera center to accept observation" -stack, b, a, no,,, "Stack multiple observations into single On/Off PHA (if no produce run-by-run files)" +bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation" +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" +stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?" # # Standard parameters diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py index 60cd3c7..b35dc98 100755 --- a/cscripts/csphagen.py +++ b/cscripts/csphagen.py @@ -3,7 +3,7 @@ # Computes the PHA spectra for source/background and ARF/RMF files using the # reflected region method # -# Copyright (C) 2017- Luigi Tibaldo +# Copyright (C) 2017 Luigi Tibaldo # # 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 @@ -30,7 +30,7 @@ import sys # =============== # class csphagen(ctools.cscript): """ - Generate pha files for classical IACT spectral analysis + Generate PHA, ARF and RMF files for classical IACT spectral analysis """ # Constructor @@ -48,6 +48,15 @@ class csphagen(ctools.cscript): # Initialise application by calling the appropriate class constructor self._init_cscript(argv) + # Initialise other variables + self._outobs = gammalib.GObservations() + self._ebounds = gammalib.GEbounds() + self._src_dir = gammalib.GSkyDir() + self._src_reg = gammalib.GSkyRegions() + self._bkg_regs = [] + self._excl_reg = None + self._has_exclusion = False + # Return return @@ -63,49 +72,43 @@ class csphagen(ctools.cscript): # Set energy bounds self._ebounds = self._create_ebounds() - # Initialise background estimation method - self._bkgmethod = self["bkgmethod"].string() - if self._bkgmethod == "REFLECTED": - self._bkgregmin = self["bkgregmin"].integer() - self._maxoffset = self["maxoffset"].real() - # Initialise source position/region querying relevant parameters self._src_dir = gammalib.GSkyDir() self._src_reg = gammalib.GSkyRegions() coordsys = self['coordsys'].string() - if coordsys == "CEL": - ra = self['ra'].real() + if coordsys == 'CEL': + ra = self['ra'].real() dec = self['dec'].real() self._src_dir.radec_deg(ra, dec) - elif coordsys == "GAL": + elif coordsys == 'GAL': glon = self['glon'].real() glat = self['glat'].real() self._src_dir.lb_deg(glon, glat) self._srcshape = self['srcshape'].string() - if self._srcshape == "CIRCLE": + if self._srcshape == 'CIRCLE': self._rad = self['rad'].real() - self._src_reg.append( - gammalib.GSkyRegionCircle(self._src_dir, self._rad)) + self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad)) + + # Query background estimation method and parameters + bkgmethod = self['bkgmethod'].string() + if bkgmethod == 'REFLECTED': + self['bkgregmin'].integer() + self['maxoffset'].real() - # exclusion map - if self["exclusion"].filename().is_fits(): - self._excl_reg = gammalib.GSkyRegionMap( - self["exclusion"].filename()) + # Exclusion map + if self['inexclusion'].filename().is_fits(): + self._excl_reg = gammalib.GSkyRegionMap(self['inexclusion'].filename()) self._has_exclusion = True else: self._has_exclusion = False - # Stacking - self._stack = self['stack'].boolean() + # Query remaining parametersStacking + self['stack'].boolean() # Query ahead output parameters if (self._read_ahead()): - self._outroot = self['outroot'].string() - - # Set some fixed parameters - self._chatter = self["chatter"].integer() - self._clobber = self["clobber"].boolean() - self._debug = self["debug"].boolean() + self['outobs'].string() + self['prefix'].string() # If there are no observations in container then get them from the # parameter file @@ -121,50 +124,66 @@ class csphagen(ctools.cscript): def _reflected_regions(self, obs): """ Calculate list of reflected regions for a single observation (pointing) - :param obs: observation - :return: list of reflected regions + + Parameters + ---------- + obs : `~gammalib.GObservation()` + CTA observation + + Returns + ------- + regions : `~gammalib.GSkyRegions` + List of reflected regions """ - outregions = [] + # Initialise list of reflected regions + regions = gammalib.GSkyRegions() + + # Get offset angle of source pnt_dir = obs.pointing().dir() - offset = pnt_dir.dist_deg(self._src_dir) - if offset <= self._rad or offset >= self["maxoffset"].real(): - if self._logExplicit(): - msg = 'Observation {} pointed at {} deg from source\n'.format( - obs.id(), offset) - self._log(msg) + offset = pnt_dir.dist_deg(self._src_dir) + + # If ... + if offset <= self._rad or offset >= self['maxoffset'].real(): + self._log_string(gammalib.EXPLICIT, 'Observation {} pointed at {} ' + 'deg from source'.format(obs.id(), offset)) + + # ... otherwise else: posang = pnt_dir.posang_deg(self._src_dir) - if self._srcshape == "CIRCLE": + if self._srcshape == 'CIRCLE': # angular separation of reflected regions wrt camera center # and number alpha = 1.05 * 2 * self._rad / offset - # 1.05 ensures background regions do not overlap due to numerical precision issues - N = int(2 * math.pi / alpha) - if N < self._bkgregmin + 3: - if self._logExplicit(): - msg = 'Observation {}: insufficient regions for background estimation\n'.format( - obs.id(), offset) - self._log = (msg) + # 1.05 ensures background regions do not overlap due to + # numerical precision issues + N = int(2.0 * math.pi / alpha) + if N < self['bkgregmin'].integer() + 3: + self._log_string(gammalib.EXPLICIT, 'Observation {}: ' + 'insufficient regions for background ' + 'estimation'.format(obs.id(), offset)) + + # Otherwise loop over position angle to create reflected + # regions else: - alpha = 360. / N - # loop to create reflected regions + alpha = 360.0 / N for s in range(2, N - 1): - dphi = s * alpha + dphi = s * alpha ctr_dir = pnt_dir.clone() ctr_dir.rotate_deg(posang + dphi, offset) region = gammalib.GSkyRegionCircle(ctr_dir, self._rad) if self._has_exclusion: if self._excl_reg.overlaps(region): - if self._logVerbose(): - msg = 'Observation {}: reflected region overlaps with exclusion region\n'.format( - obs.id(), offset) - self._log(msg) + self._log_string(gammalib.VERBOSE, 'Observation ' + '{}: reflected region overlaps ' + 'with exclusion region'. + format(obs.id(), offset)) else: - outregions.append(region) + regions.append(region) else: - outregions.append(region) + regions.append(region) - return outregions + # Return reflected regions + return regions # Public methods def run(self): @@ -181,74 +200,135 @@ class csphagen(ctools.cscript): # Write observation into logger self._log_observations(gammalib.NORMAL, self._obs, 'Observation') + # Write header self._log_header1(gammalib.NORMAL, 'Generation of source and background spectra') + + # Initialise run variables + self._outobs = gammalib.GObservations() + self._bkg_regs = [] + # Loop through observations and generate pha, arf, rmf files - # Save ON/OFF regions in ds9 reg files - outobs = gammalib.GObservations() - etrue = self._ebounds - ereco = self._ebounds for obs in self._obs: + + # Initialise background regions for this observation bkg_reg = gammalib.GSkyRegions() - if self._bkgmethod == "REFLECTED": - regions = self._reflected_regions(obs) - for region in regions: - bkg_reg.append(region) - if bkg_reg.size() >= 1 or ( - self._bkgmethod == "REFLECTED" and bkg_reg.size() >= self._bkgregmin): - onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, etrue, - ereco, self._src_reg, + + # If reflected background is requested then created reflected + # background regions + if self['bkgmethod'].string() == 'REFLECTED': + bkg_reg = self._reflected_regions(obs) + + # If there are reflected regions then create On/Off observation + # and append it to the output container + if bkg_reg.size() >= self['bkgregmin'].integer(): + onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, + self._ebounds, + self._ebounds, + self._src_reg, bkg_reg) onoff.id(obs.id()) - outobs.append(onoff) - self._src_reg.save( - self._outroot + '_{}_on.reg'.format(obs.id())) - bkg_reg.save(self._outroot + '_{}_off.reg'.format(obs.id())) + self._outobs.append(onoff) + self._bkg_regs.append({'regions': bkg_reg, 'id': obs.id()}) else: - if self._logNormal(): - msg = 'Observation {} not included in spectra generation\n'.format( - obs.id()) - self._log(msg) - - # Save PHA, ARF and RMFs - if outobs.size() < 2 or self._stack == False: - for obs in outobs: - obs.on_spec().save( - self._outroot + '_{}_pha_on.fits'.format(obs.id()), - True) - obs.off_spec().save( - self._outroot + '_{}_pha_off.fits'.format(obs.id()), - True) - obs.arf().save(self._outroot + '_{}_arf.fits'.format(obs.id()), - True) - obs.rmf().save(self._outroot + '_{}_rmf.fits'.format(obs.id()), - True) - else: - if self._logVerbose(): - msg = 'Stacking {} observations\n'.format( - outobs.size()) - self._log(msg) - stackobs = gammalib.GCTAOnOffObservation(outobs) - stackobs.on_spec().save( - self._outroot + '_stacked_pha_on.fits', - True) - stackobs.off_spec().save( - self._outroot + '_stacked_pha_off.fits', - True) - stackobs.arf().save(self._outroot + '_stacked_arf.fits', - True) - stackobs.rmf().save(self._outroot + '_stacked_rmf.fits', - True) - outobs = gammalib.GObservations() - outobs.append(stackobs) - - # Save On/Off observations - outname = self._outroot + '.xml' - outobs.save(outname) - # Log filename - if self._logNormal(): - self._log( - 'Output observation definition XML file: {}\n'.format(outname)) + self._log_string(gammalib.NORMAL, 'Observation {} not included ' + 'in spectra generation'.format(obs.id())) + + # Stack observations + if self._outobs.size() > 1 and self['stack'].boolean() == True: + + # Write header + self._log_header1(gammalib.NORMAL, + 'Stacking {} observations'.format(self._outobs.size())) + + # Stack observations + stacked_obs = gammalib.GCTAOnOffObservation(self._outobs) + + # Put stacked observations in output container + self._outobs = gammalib.GObservations() + self._outobs.append(stacked_obs) + + # Return + return + + def save(self): + """ + Save data + """ + # Write header + self._log_header1(gammalib.TERSE, 'Save data') + + # Get XMK output filename, prefix and clobber + outobs = self['outobs'].string() + prefix = self['prefix'].string() + clobber = self['clobber'].boolean() + + # Loop over all observation in container + for obs in self._outobs: + + # Set filenames + if self['stack'].boolean(): + onname = prefix + '_stacked_pha_on.fits' + offname = prefix + '_stacked_pha_off.fits' + arfname = prefix + '_stacked_arf.fits' + rmfname = prefix + '_stacked_rmf.fits' + elif self._outobs.size() > 1: + onname = prefix + '_{}_pha_on.fits'.format(obs.id()) + offname = prefix + '_{}_pha_off.fits'.format(obs.id()) + arfname = prefix + '_{}_arf.fits'.format(obs.id()) + rmfname = prefix + '_{}_rmf.fits'.format(obs.id()) + else: + onname = prefix + '_pha_on.fits' + offname = prefix + '_pha_off.fits' + arfname = prefix + '_arf.fits' + rmfname = prefix + '_rmf.fits' + + # Save files + obs.on_spec().save(onname, clobber) + obs.off_spec().save(offname, clobber) + obs.arf().save(arfname, clobber) + obs.rmf().save(rmfname, clobber) + + # Log file names + self._log_value(gammalib.NORMAL, 'PHA on file', onname) + self._log_value(gammalib.NORMAL, 'PHA off file', offname) + self._log_value(gammalib.NORMAL, 'ARF file', arfname) + self._log_value(gammalib.NORMAL, 'RMF file', arfname) + + # Save observation definition XML file + self._outobs.save(outobs) + + # Log file name + self._log_value(gammalib.NORMAL, 'Obs. definition XML file', outobs) + + # Save ds9 On region file + regname = prefix + '_on.reg' + self._src_reg.save(regname) + self._log_value(gammalib.NORMAL, 'On region file', regname) + + # Save ds9 Off region files + for bkg_reg in self._bkg_regs: + + # Set filename + if len(self._bkg_regs) > 1: + regname = prefix + '_{}_off.reg'.format(bkg_reg['id']) + else: + regname = prefix + '_off.reg' + + # Save ds9 region file + bkg_reg['regions'].save(regname) + + # Log file name + self._log_value(gammalib.NORMAL, 'Off region file', regname) + + # Return + return + + def obs(self): + """ + Return observation container + """ + return self._outobs def execute(self): """ @@ -263,6 +343,9 @@ class csphagen(ctools.cscript): # Run the script self.run() + # Save result + self.save() + # Return return @@ -271,6 +354,7 @@ class csphagen(ctools.cscript): # Main routine entry point # # ======================== # if __name__ == '__main__': + # Create instance of application app = csphagen(sys.argv) diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index 9e4ab6b..07cd29e 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -10,6 +10,7 @@ ctools 1.5 is a major release that adds significant functionality. In particular, this release provides: +* Support of classical IACT analysis * Specifcation of times in various formats (UTC, MJD, JD, MET) @@ -28,5 +29,7 @@ Improvements New features ------------ +* [`2208 `_] - + Add :ref:`csphagen` script * [`1864 `_] - Use t-type parameters for all ``tmin`` and ``tmax`` parameters diff --git a/doc/source/users/reference_manual/csphagen.rst b/doc/source/users/reference_manual/csphagen.rst index 9117888..9e9fb94 100644 --- a/doc/source/users/reference_manual/csphagen.rst +++ b/doc/source/users/reference_manual/csphagen.rst @@ -1,57 +1,73 @@ .. _csphagen: csphagen -========= +======== -Creates files necessary to perform a classical region-based spectral analysis on IACT data, i.e., source and background PHA files, ARF and RMF files. +Creates files necessary to perform a classical region-based spectral analysis +on IACT data, i.e., source and background PHA files, ARF and RMF files. Synopsis -------- -This script can be used to derive from an observation or a set of observations the count spectra in a source region and in background regions, as well as the detector response (effective area, energy redistribution matrix), that can be used to perform a classical 1D spectral analysis. The output files are saved in the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF). `See here `__. +This script can be used to derive from an observation or a set of observations +the count spectra in a source region and in background regions, as well as the +detector response (effective area, energy redistribution matrix), that can be +used to perform a classical 1D spectral analysis. The output files are saved in +the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF). +`See here `__. -The outputs are: 1) the PHA, ARF, RMF files, either separately for each observation, or stacked for a set of observations, 2) DS9 regions files listing the source and background regions for each observation; 3) a new observation definition XML file. +The outputs are: +1) the PHA, ARF, RMF files, either separately for each observation, or stacked + for a set of observations, +2) DS9 regions files listing the source and background regions for each + observation; +3) a new observation definition XML file. General parameters ------------------ ``inobs [file]`` - Input event list or observation definition XML file + Input event list or observation definition XML file. ``caldb [string]`` Calibration database. - + ``irf [string]`` Instrument response function. -``(exclusion = NONE) [file]`` - Optional FITS file containing a WCS map in the first hdu that defines sky regions not to be used for background estimation (where map value > 0) +``(inexclusion = NONE) [file]`` + Optional FITS file containing a WCS map in the first hdu that defines sky + regions not to be used for background estimation (where map value != 0). -``outroot [string]`` - Root of the file name for output PHA, ARF, RMF, XML, and ds9 reg files +``outobs [string]`` + Output observation definition XML file. + +``(prefix = onoff) [string]`` + Prefix of the file name for output PHA, ARF, RMF, XML, and DS9 region files. ``emin [real]`` - Lower energy limit (TeV) if LIN or LOG binning algorithms are used + Lower energy limit (TeV) if LIN or LOG binning algorithms are used. ``emax [real]`` - Upper energy limit (TeV) if LIN or LOG binning algorithms are used + Upper energy limit (TeV) if LIN or LOG binning algorithms are used. ``enumbins [integer]`` - Number of energy bins. At least 30 bins per decade are recommended for proper evaluation of the instrument response + Number of energy bins. At least 30 bins per decade are recommended for + proper evaluation of the instrument response. ``ebinalg [string]’’ - Energy binning algorithm: LOG - logarithmically spaced energy bins; LIN - linearly spaced energy bins; FILE - energy bounds retrieved from file + Energy binning algorithm: LOG - logarithmically spaced energy bins; + LIN - linearly spaced energy bins; FILE - energy bounds retrieved from file. ``ebinfile [file]`` - Name of the file containing the energy bin definition if FILE algorithm is used - -``(bkgmethod = REFLECTED) [string]’’ - Method for background estimation. REFLECTED: background evaluated in regions with the same shape as source region reflected w.r.t. pointing direction for each observation + Name of the file containing the energy bin definition if FILE algorithm + is used. ``(srcshape = CIRCLE) [string]`` - Shape of the source region. CIRCLE: circular region around given position + Shape of the source region. + CIRCLE: circular region around given position. ``coordsys [string]`` Coordinate system (CEL - celestial, GAL - galactic) @@ -71,14 +87,23 @@ General parameters ``rad [real]`` Radius of source region circle (deg) +``(bkgmethod = REFLECTED) [string]’’ + Method for background estimation. + REFLECTED: background evaluated in regions with the same shape as source + region reflected w.r.t. pointing direction for each observation. + ``(bkgregmin = 2) [integer]`` - Minimum number of background regions for REFLECTED method. If not available, observation is skipped. (Default: 2) + Minimum number of background regions that are required for an observation. + If this number of background regions is not available the observation is + skipped. -``(maxoffset = 4.) [real]`` - Maximum offset of source from camera center to accept observation +``(maxoffset = 4.0) [real]`` + Maximum offset in degrees of source from camera center to accept the + observation. ``stack [boolean]`` - Yes: Stack multiple observations into single PHA. No: produce run-by-run PHA files + Specifies whether multiple observations should be stacked (``yes``) or + whether run-wise PHA, ARF and RMF files should be produced (``no``). Standard parameters @@ -100,13 +125,13 @@ Standard parameters Specifies whether an existing output runlist should be overwritten. ``(debug = no) [boolean]`` - Enables debug mode. In debug mode the executable will dump any log file output to the console. + Enables debug mode. In debug mode the executable will dump any log file + output to the console. ``(mode = ql) [string]`` Mode of automatic parameters (default is "ql", i.e. "query and learn"). -``(logfile = csfindobs.log) [filename]`` +``(logfile = csphagen.log) [filename]`` Log filename. - diff --git a/test/test_csphagen.py b/test/test_csphagen.py index f06f045..bf61270 100644 --- a/test/test_csphagen.py +++ b/test/test_csphagen.py @@ -41,12 +41,13 @@ class Test(test): """ # Call base class constructor test.__init__(self) - self._myevents1 = self._datadir + '/crab_offaxis1.fits' - self._myevents2 = self._datadir + '/crab_offaxis2.fits' - self._exclusion = self._datadir + '/crab_exclusion.fits' + self._myevents1 = self._datadir + '/crab_offaxis1.fits' + self._myevents2 = self._datadir + '/crab_offaxis2.fits' + self._exclusion = self._datadir + '/crab_exclusion.fits' self._nreg_with_excl = 5 - self._nreg_wo_excl = 8 - self._nreg_mul = [self._nreg_with_excl, 6] + self._nreg_wo_excl = 8 + self._nreg_mul = [self._nreg_with_excl, 6] + self._nbins = 20 # number of expected background regions with/wo exclusion, # and for two different runs @@ -75,14 +76,17 @@ class Test(test): """ # Set script name csphagen = self._script('csphagen') - nbins = 120 - cmd = csphagen + ' inobs="' + self._myevents1 + '" caldb="' + \ - self._caldb + '" irf="' + self._irf + \ - '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \ - str(nbins) + ' coordsys="CEL" ra=83.633 dec=22.0145' + \ - ' rad=0.2 stack="no" exclusion="' + self._exclusion + \ - '" outroot="phagen_cmd1" logfile="csphagen_cmd1.log" chatter=1' + # Setup csphagen command + cmd = csphagen + ' inobs="' + self._myevents1 + \ + '" caldb="' + self._caldb + '" irf="' + self._irf + \ + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \ + str(self._nbins) + \ + ' coordsys="CEL" ra=83.633 dec=22.0145' + \ + ' rad=0.2 stack="no" inexclusion="' + \ + self._exclusion + \ + '" outobs="csphagen_cmd1.xml" prefix="csphagen_cmd1" ' + \ + 'logfile="csphagen_cmd1.log" chatter=2' # Check if execution of wrong command fails self.test_assert(self._execute('command_that_does_not_exist') != 0, @@ -93,15 +97,19 @@ class Test(test): 'Check successful execution from command line') # Check output files - self._check_output('phagen_cmd1_', nbins, self._nreg_with_excl) - self._check_outobs('phagen_cmd1', 1) - - cmd = csphagen + ' inobs="events_that_do_not_exist.fits" caldb="' + \ - self._caldb + '" irf="' + self._irf + \ - '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \ - str(nbins) + ' coordsys="CEL" ra=83.633 dec=22.0145' + \ - ' rad=0.2 stack="no" exclusion="' + self._exclusion + \ - '" outroot="phagen_cmd1" logfile="csphagen_cmd1.log" chatter=1' + self._check_output('csphagen_cmd1', self._nbins, self._nreg_with_excl) + self._check_outobs('csphagen_cmd1', 1) + + # Setup csphagen command + cmd = csphagen + ' inobs="events_that_do_not_exist.fits" ' + \ + 'caldb="' + self._caldb + '" irf="' + self._irf + \ + '" ebinalg="LOG" emin=0.1 emax=100. enumbins=' + \ + str(self._nbins) + \ + ' coordsys="CEL" ra=83.633 dec=22.0145' + \ + ' rad=0.2 stack="no" inexclusion="' + \ + self._exclusion + \ + '" outobs="csphagen_cmd2.xml" prefix="csphagen_cmd2" ' + \ + 'logfile="csphagen_cmd2.log" chatter=1' # Check if execution failed self.test_assert(self._execute(cmd) != 0, @@ -117,58 +125,58 @@ class Test(test): """ Test csphagen from Python """ - nbins = 120 - # Same test as from command line phagen = cscripts.csphagen() - phagen['inobs'] = self._myevents1 - phagen['caldb'] = self._caldb - phagen['irf'] = self._irf - phagen['ebinalg'] = 'LOG' - phagen['emin'] = 0.1 - phagen['emax'] = 100. - phagen['enumbins'] = nbins - phagen['coordsys'] = 'CEL' - phagen['ra'] = 83.633 - phagen['dec'] = 22.0145 - phagen['rad'] = 0.2 - phagen['stack'] = False - phagen['exclusion'] = self._exclusion - phagen['outroot'] = 'phagen_py1' - phagen['logfile'] = 'csphagen_py1.log' - phagen['chatter'] = 1 + phagen['inobs'] = self._myevents1 + phagen['caldb'] = self._caldb + phagen['irf'] = self._irf + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins + phagen['coordsys'] = 'CEL' + phagen['ra'] = 83.633 + phagen['dec'] = 22.0145 + phagen['rad'] = 0.2 + phagen['stack'] = False + phagen['inexclusion'] = self._exclusion + phagen['outobs'] = 'csphagen_py1.xml' + phagen['prefix'] = 'csphagen_py1' + phagen['logfile'] = 'csphagen_py1.log' + phagen['chatter'] = 1 # Run script phagen.execute() # Check output - self._check_output('phagen_py1_', nbins, self._nreg_with_excl) - self._check_outobs('phagen_py1', 1) + self._check_output('csphagen_py1', self._nbins, self._nreg_with_excl) + self._check_outobs('csphagen_py1', 1) # Now test without exclusion region phagen = cscripts.csphagen() - phagen['inobs'] = self._myevents1 - phagen['caldb'] = self._caldb - phagen['irf'] = self._irf - phagen['ebinalg'] = 'LOG' - phagen['emin'] = 0.1 - phagen['emax'] = 100. - phagen['enumbins'] = nbins + phagen['inobs'] = self._myevents1 + phagen['caldb'] = self._caldb + phagen['irf'] = self._irf + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins phagen['coordsys'] = 'CEL' - phagen['ra'] = 83.633 - phagen['dec'] = 22.0145 - phagen['rad'] = 0.2 - phagen['stack'] = False - phagen['outroot'] = 'phagen_py2' - phagen['logfile'] = 'csphagen_py2.log' - phagen['chatter'] = 1 + phagen['ra'] = 83.633 + phagen['dec'] = 22.0145 + phagen['rad'] = 0.2 + phagen['stack'] = False + phagen['outobs'] = 'csphagen_py2.xml' + phagen['prefix'] = 'csphagen_py2' + phagen['logfile'] = 'csphagen_py2.log' + phagen['chatter'] = 2 # Run script phagen.execute() # Check output - self._check_output('phagen_py2_', nbins, self._nreg_wo_excl) - self._check_outobs('phagen_py2', 1) + self._check_output('csphagen_py2', self._nbins, self._nreg_wo_excl) + self._check_outobs('csphagen_py2', 1) # Test with multiple input observations, no stacking # Create observation container @@ -181,55 +189,57 @@ class Test(test): # Setup csphagen phagen = cscripts.csphagen(obs) - phagen['ebinalg'] = 'LOG' - phagen['emin'] = 0.1 - phagen['emax'] = 100. - phagen['enumbins'] = nbins - phagen['coordsys'] = 'CEL' - phagen['ra'] = 83.633 - phagen['dec'] = 22.0145 - phagen['rad'] = 0.2 - phagen['stack'] = False - phagen['exclusion'] = self._exclusion - phagen['outroot'] = 'phagen_py3' - phagen['logfile'] = 'csphagen_py3.log' - phagen['chatter'] = 1 + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins + phagen['coordsys'] = 'CEL' + phagen['ra'] = 83.633 + phagen['dec'] = 22.0145 + phagen['rad'] = 0.2 + phagen['stack'] = False + phagen['inexclusion'] = self._exclusion + phagen['outobs'] = 'csphagen_py3.xml' + phagen['prefix'] = 'csphagen_py3' + phagen['logfile'] = 'csphagen_py3.log' + phagen['chatter'] = 3 # Run script phagen.execute() # Check output for s in range(2): - self._check_output('phagen_py3_' + str(s + 1), nbins, + self._check_output('csphagen_py3_' + str(s + 1), self._nbins, self._nreg_mul[s]) - self._check_outobs('phagen_py3', 2) + self._check_outobs('csphagen_py3', 2) # Test with multiple input observations and stacking # Setup csphagen phagen = cscripts.csphagen(obs) - phagen['ebinalg'] = 'LOG' - phagen['emin'] = 0.1 - phagen['emax'] = 100. - phagen['enumbins'] = nbins - phagen['coordsys'] = 'CEL' - phagen['ra'] = 83.633 - phagen['dec'] = 22.0145 - phagen['rad'] = 0.2 - phagen['stack'] = True - phagen['exclusion'] = self._exclusion - phagen['outroot'] = 'phagen_py4' - phagen['logfile'] = 'csphagen_py4.log' - phagen['chatter'] = 1 + phagen['ebinalg'] = 'LOG' + phagen['emin'] = 0.1 + phagen['emax'] = 100.0 + phagen['enumbins'] = self._nbins + phagen['coordsys'] = 'CEL' + phagen['ra'] = 83.633 + phagen['dec'] = 22.0145 + phagen['rad'] = 0.2 + phagen['stack'] = True + phagen['inexclusion'] = self._exclusion + phagen['outobs'] = 'csphagen_py4.xml' + phagen['prefix'] = 'csphagen_py4' + phagen['logfile'] = 'csphagen_py4.log' + phagen['chatter'] = 4 # Run script phagen.execute() # Check output for s in range(2): - self._check_output('phagen_py4_stacked', nbins, + self._check_output('csphagen_py4_stacked', self._nbins, 0, check_regions=False) - self._check_outobs('phagen_py4', 1) + self._check_outobs('csphagen_py4', 1) return @@ -274,8 +284,7 @@ class Test(test): # Check FITS table structure self.test_value(table.ncols(), len(cols), - 'Check for %d columns in PHA table' % len( - cols)) + 'Check for %d columns in PHA table' % len(cols)) self.test_value(table.nrows(), bins, 'Check for %d rows in PHA table' % bins) for col in cols: @@ -313,8 +322,7 @@ class Test(test): # Check FITS table structure self.test_value(table.ncols(), len(cols), - 'Check for %d columns in ARF table' % len( - cols)) + 'Check for %d columns in ARF table' % len(cols)) self.test_value(table.nrows(), bins, 'Check for %d rows in ARF table' % bins) for col in cols: @@ -350,8 +358,7 @@ class Test(test): # Check FITS table structure self.test_value(table.ncols(), len(cols), - 'Check for %d columns in RMF table' % len( - cols)) + 'Check for %d columns in RMF table' % len(cols)) self.test_value(table.nrows(), bins, 'Check for %d rows in RMF table' % bins) for col in cols: @@ -372,33 +379,39 @@ class Test(test): """ Check the output from a csphagen run """ - # OGIP files self._check_pha(filenameroot + '_pha_on.fits', bins) self._check_pha(filenameroot + '_pha_off.fits', bins) self._check_arf(filenameroot + '_arf.fits', bins) self._check_rmf(filenameroot + '_rmf.fits', bins) - # Observations + # Optionally check for regions + if check_regions: + # Check On region + onregion = (filenameroot+'_on.reg').replace('_1_', '_').replace('_2_', '_') + reg = gammalib.GSkyRegions(onregion) + self.test_value(reg.size(), 1, 'Check for 1 region in source region file') - # Regions - if check_regions: - reg = gammalib.GSkyRegions(filenameroot + "_on.reg") - self.test_value(reg.size(), 1, 'Check for ' + str( - 1) + ' region in source region file') - reg = gammalib.GSkyRegions(filenameroot + "_off.reg") - self.test_value(reg.size(), nreg, 'Check for ' + str( - nreg) + ' region in background region file') + # Check for Off region + offregion = filenameroot + '_off.reg' + reg = gammalib.GSkyRegions(offregion) + self.test_value(reg.size(), nreg, 'Check for ' + str(nreg) + + ' region in background region file') + # Return return def _check_outobs(self, filenameroot, nout): """ Check the output XML file containing ON/OFF observations """ + # Load observation container obs = gammalib.GObservations(filenameroot + '.xml') - self.test_value(obs.size(), nout, 'Check for ' + str( - nout) + ' observations in XML file') + # Check container size + self.test_value(obs.size(), nout, 'Check for ' + str(nout) + + ' observations in XML file') + + # Return return -- libgit2 0.21.4